PS2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1> EEE485 HW3 </h1></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Kazim Sanlav**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T22:03:34.241982Z",
     "start_time": "2018-12-12T22:03:33.379496Z"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.io as sio\n",
    "import numpy as np\n",
    "from numpy.linalg import inv\n",
    "from sklearn import preprocessing\n",
    "from numpy import linalg as LA \n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import Latex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T22:03:35.713878Z",
     "start_time": "2018-12-12T22:03:35.653644Z"
    }
   },
   "outputs": [],
   "source": [
    "q2 = pd.read_csv('student.csv', header=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T22:03:36.553016Z",
     "start_time": "2018-12-12T22:03:36.530566Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1426.0</td>\n",
       "      <td>3.61</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1277.5</td>\n",
       "      <td>3.67</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>2104.7</td>\n",
       "      <td>4.00</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1534.0</td>\n",
       "      <td>3.19</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>2634.8</td>\n",
       "      <td>2.93</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   0       1     2  3\n",
       "0  1  1426.0  3.61  3\n",
       "1  1  1277.5  3.67  3\n",
       "2  1  2104.7  4.00  1\n",
       "3  1  1534.0  3.19  4\n",
       "4  1  2634.8  2.93  4"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T22:03:38.860388Z",
     "start_time": "2018-12-12T22:03:38.844133Z"
    }
   },
   "outputs": [],
   "source": [
    "def logit(x):\n",
    "    '''\n",
    "    x: vector of input\n",
    "    '''\n",
    "    return x/(1+x) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T22:03:39.621138Z",
     "start_time": "2018-12-12T22:03:39.606482Z"
    }
   },
   "outputs": [],
   "source": [
    "def gradient_decent(X,y):\n",
    "    '''\n",
    "    Gradient decent algorithm\n",
    "    X:  Design matrix [1 x_1;...;1 x_n]\n",
    "    y:  True response [y_1,...,y_n]\n",
    "    pi: [pi_1...pi_n] where pi_i = Pr(Y_i=1|X_i=x_i)\n",
    "    '''\n",
    "\n",
    "    treshhold = 0.01\n",
    "    lrate = 10^-1 #learning rate\n",
    "    converged = False\n",
    "    beta = np.zeros((2,1)) # inital weights\n",
    "    i = 0 # step number\n",
    "    \n",
    "    \n",
    "    while not converged and i < 100:\n",
    "        i+=1 \n",
    "        #beta_n+1 = beta_n - lamda*gradient_NLL_beta_n\n",
    "        #gradient = X.T(pi-y)\n",
    "        pi = logit(np.exp(X@beta)) #pi\n",
    "        gradient_NLL = (X.T@(pi-y)) #gradient\n",
    "        beta -= lrate*gradient_NLL #parameter update\n",
    "        \n",
    "        if (LA.norm(gradient_NLL)<treshhold):\n",
    "            converged = True\n",
    "        if (i%10 == 0):\n",
    "            print(gradient_NLL)\n",
    "#             print('step:',i,' beta:',beta)\n",
    "    return beta #maximum likelihood estimate of beta\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T22:03:40.470602Z",
     "start_time": "2018-12-12T22:03:40.445539Z"
    }
   },
   "outputs": [],
   "source": [
    "# add ones to first column of X\n",
    "q2_X = np.asarray(q2.loc[:,1])\n",
    "q2_X/=1000\n",
    "a = np.ones((400,2))\n",
    "a[:,1] = q2_X\n",
    "q2_X = a\n",
    "\n",
    "q2_Y = np.asarray(q2.loc[:,0])\n",
    "q2_Y = np.reshape(q2_Y,(400,1)) #response variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T22:03:41.393855Z",
     "start_time": "2018-12-12T22:03:41.284371Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n",
      "[[-200.     ]\n",
      " [-334.67436]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[-217712.00000188],\n",
       "       [-365521.61664147]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gradient_decent(q2_X,q2_Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = sio.loadmat('PS2Q3train.mat')\n",
    "df_train = np.concatenate((df_train['x_tr'],df_train['y_tr']), axis = 1)\n",
    "df_train = pd.DataFrame(df_train, columns=['X','Y'])\n",
    "df_test = sio.loadmat('PS2Q3test.mat')\n",
    "df_test = np.concatenate((df_test['x_te'],df_test['y_te']), axis = 1)\n",
    "df_test = pd.DataFrame(df_test, columns=['X','Y'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.627295</td>\n",
       "      <td>0.304737</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-1.859568</td>\n",
       "      <td>-0.951721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-1.761680</td>\n",
       "      <td>-1.059035</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.210733</td>\n",
       "      <td>0.049572</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.679922</td>\n",
       "      <td>-0.688042</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          X         Y\n",
       "0 -0.627295  0.304737\n",
       "1 -1.859568 -0.951721\n",
       "2 -1.761680 -1.059035\n",
       "3  0.210733  0.049572\n",
       "4 -1.679922 -0.688042"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Train Data')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXuUXdd52Pf79jnnvueBGQwGb4AgCAoiTMoSFZq1raq2ZCNWLC+7TuW0dp0sZyl1Y8d145UuZyVWm9Sp47jxiuXEDmt5RX7Kcey0dC1T1sMqTZmkRFISBQIECALEczDv+77nuXf/2Ode3BnM4EFiMIC0f2vNmvs495zvnnPP/vb33GKMweFwOByOm0VttgAOh8PhuLdwisPhcDgct4RTHA6Hw+G4JZzicDgcDsct4RSHw+FwOG4JpzgcDofDcUs4xeFwvAVExBORtojs3WxZHI47hVMcjm8o8kG+/6dFpDf0/L+71f0ZYzJjTM0Yc/5NyHJQRMzQ8a+IyJ+IyHfewj7+roh8/laP7XC8FZzicHxDkQ/yNWNMDTgPfO/Qa7+7ensR8e+UTMA3A58DnhSRH97o4zocbxanOByOIUTkfxeRPxCR3xeRFvDDIvK4iDwnInURmRGRXxGRIN/ez62G/fnz38nf/zMRaYnIsyJy380c2xgzY4z5ZeCfA78oIpLv85+IyJl8f6+IyAfz178J+FXg23OLZSF//YMi8pV8+/Mi8k9v82lyfIPjFIfDcS3fD/weMAb8AZACPwVsBb4VOAr8vet8/r8F/ikwgbVq/vktHv+PgR3Awfz5qfy4Y8DPA78nItPGmK8BPwH8ZW61bM23bwM/nG//vcBPicjfuEUZHI51cYrD4biWZ4wxf2KM0caYnjHmS8aY540xqTHmDPAE8F9e5/P/yRjzgjEmAX4XeMctHv9y/n8CwBjzH3NrRBtjfg94A3h0vQ8bYz5njDmWb/9V4BM3kNfhuCWc4nA4ruXC8BMReZuI/GkevG4C/wxrfazHlaHHXaB2i8fflf9fyo//t0Xkq7mrrA687XrHz11rnxeReRFpAH/3BvI6HLeEUxwOx7Wsbhn974FjwEFjzCjwc4Bs4PG/H6t8TovIAeDXgB8HJo0x48CrQ8dfq731J4A/AvYYY8aA39hgeR3fYDjF4XDcmBGgAXRE5DDXj2+8aURkWkT+AfBPgP/F2DUPaljlMG83kb+LtTj6zAK7+8H6IXmXjDGhiHwL8EMbIa/jGxenOByOG/MPgR8FWljr4w9u5877dRzAy8B3Az9gjPktAGPMy8CvAF8EZrBK4/mhj38aeA2YFZG+i+zHgf8jzwr7x8B/vJ3yOhziFnJyOBwOx63gLA6Hw+Fw3BJOcTgcDofjlnCKw+FwOBy3xKYpDhHZIyJ/ISIn8jYKP7XGNpK3bzgtIi+LyDs3Q1aHw+FwXGXDG7hdhxT4h8aYl0RkBHhRRD5tjDk+tM1fBx7I/x7D5rM/dqMdb9261ezfv38DRHY4HI6vT1588cUFY8zUzWy7aYrDGDODTS/EGNMSkRPYitlhxfF9wG/l+ezPici4iOzIP7su+/fv54UXXtgo0R0Oh+PrDhE5d7Pb3hUxjryz6DezMj8drCIZbv9wkavtGBwOh8OxCWy64hCRGrY9wv9kjGmufnuNj6xZeCIiHxaRF0Tkhfn5+dstpsPhcDhyNlVx5G0S/gj4XWPMH6+xyUVgz9Dz3VztHLoCY8wTxphHjTGPTk3dlJvO4XA4HG+CzcyqEuBjwAljzL9eZ7Mngf8+z676FqBxo/iGw+FwODaWzcyq+lbgR4CvichX8tf+MbAXwBjz68Ange8BTmPbU/+dTZDT4XA4HENsZlbVM9yg1XOeTfX374xEDofDce9xYqbBU8dmuVTvsWu8zNEj0xzeMbahx9z04LjD4XA43hwnZho88fRZGr2EHWMlGr2EJ54+y4mZxoYe1ykOh8PhuEd56tgsY+WAsXKAEhk8furY7IYe1ykOh8PhuEe5VO8xUloZcRgp+Vyq9zb0uJsZHHc4HA7Hm6Af13jx3BKdKKVW9JkaKXFwW5XA89g1Xt7Q4zvF4XA4HPcQf/ryJT762dfpxCmL7YgkMyx3YurdmIvLXd62fZQPffehDZXBKQ6Hw+G4Rzgx0+Cjn3sdBLpxQpIZUm0QoBGmJBkUPHFZVQ6Hw+GwPHVsliTTjJZ82lFGpg2eAiWAgTBJefbMosuqcjgcDoflUr3HZLVAlGoyDZIrjNRY5aFECBO94Sm5TnE4HA7HPULRE2abIa9eaZJqQ5YrDcA+1obAkw1PyXWKw+FwOO4BTsw0eG2uzWInpuR717Td0MYqj06U8dpsY0NTcp3icDgcjnuAp47NEqWaHWMlaiWfsXKw4n3J/3xP+OLZZbpRsmGyuKwqh8PhuAe4VO8RpRlbKgUmqtbeeHWmQSPMAPAUBEoQJSQpnFvaOIvDKQ6Hw+G4yzkx0+D8UpfL9ZALuovvKapFH08pRDIkj3Mk2g7qoyWPdpRumDxOcTgcDsddTL+RYSVQpJkmyTRRognjjEQbjLExh1LgkWlDpg3FwGN6tLRhMjnF4XA4HHcpJ2YafOTJ4yy2I3pJxtRIkWYvoRGmGG0IlGAwJBlESUbBV/i+oh1l/OjjezdMLqc4HA6H4y6kb2kstWNKvnBpOUYbg6eE8XKAUgIGDIZqwWepG+MpRbXgsWdLmQ88vGvDZNvsNcd/U0TmROTYOu+/V0QaIvKV/O/n7rSMDofDsRn0W6aXC4pL9QhPCZ4SkszQjlIE8JRgtG074inFg9tHeHTfFr7l/q0bKttmp+P+B+DoDbb5S2PMO/K/f3YHZHI4HI5Np98y3RjQWpNpTZhoUm0wxsYylECUGbpxRiVQNHsJX75Q59B0dUNl21TFYYx5GljaTBkcDofjbmTXeJlWmLLcjeklGUlmFQVAqqEZJvSSDAV0k5RekhEo4Z17xzk129lQ2Tbb4rgZHheRr4rIn4nIQ+ttJCIfFpEXROSF+fn5Oymfw+Fw3HaOHpnmwlKXhXaEiFAMFJ4SfIGSLwSeDYJHmWbPljL7Jqsk2lApehu+kNPdrjheAvYZYx4BPgr83+ttaIx5whjzqDHm0ampqTsmoMPhcGwEh3eMMT1apOApRMS2GDEGDcSZQRvDaMlntBQQJppS4FH0Fccvt76xF3IyxjSHHn9SRP6diGw1xixsplwOh8OxEfRX9rtU77FrvMxcK+LwjlGWu3ahpsVEo8i74mKbGooYstjGPTCG5V7K0SPTGyrnXW1xiMh2EXuKROSvYeVd3FypHA6H4/bTT79t9BJ2jJVo9BIuLveoBApPKZQoAk8w2IaGnrJWiAdUCh7tKAMRvu3g5IYv5LSpFoeI/D7wXmCriFwEPgIEAMaYXwd+EPhxEUmBHvBDxhizSeI6HA7HhtFPv+03LxwrBzw4XePklTbbRwucXWgT5T3UqwVF0fdohgkFz+O7v2krRd+n0Uv4kcf3bbism6o4jDF/6wbv/yrwq3dIHIfD4dg0LtV77Bhb2SZk72SV2WbIlVYMCKUAip5HnGm0sYV/RV+RZLBtJOBD79694dYG3OWuKofD4fhGoZ9+O0wrTNEIh7bVKAcKY4RUa0q+ougrtlQLVIt3fv7vFIfD4XDcBRw9Mk2jl9DoJWhjBo97UcIXXl8kTDSCIU41jTBloRWz0IoIPDWIiWz0krF9nOJwOByOu4DDO8b48HvuI04zPnNili+eXSKKU84vh2TaUCt6aGNIdb7aH9BNMs4vdfjcq3PEabbhS8b2uavTcR0Oh+MbjV6iObStxrFLDZ46Pos2dmW/kZKPiO2G28dXYICZ5S6ZNrxjz9iGrsPRx1kcDofDcZfw1LFZtNZ88Y0lzi/3BuuKG6AdpsSpwRPsH+QxD0MzSsm05sTMxhf/gbM4HA6H467hUr3HTCOk2UvxlZAZwBgybNGfyes34sxQUEJm7EJOCMy3YorBxhf/gbM4HA6H465h13iZxU5MqjWeCJ4Ck5sdmSEv/rNKQ3lQ8BUi4HuCr4RttaJLx3U4HI5vJI4emSbwFEps2m2/Slzl7qmSrxgt+QS+IstsrKNW8tk3UWWyVmD7+MYtFzuMUxwOh8Nxl3B4xxg/+R33M14JSLXtRRUo8JVQ8BW7xktUCj4iMFEtsKVaGPx/2/YR3n4HrA1wMQ6Hw+G4q/jAw7s4MFXjVz7zGp87OQdAwVNoY1jqJpQDRSXwKPge79w7zt7JKq0wpdFL7kh8A5zF4XA4HHcdh3eM8Ws/8ig/+M7d7J2oECaaKNF0opTlbkKUGQ5srTDTjJhphIyVAz78nvvuSHwDnMXhcDgcdwWrW6ofmq5yuRFyablHovPajczgG4OvhIVOzOEdY/zS33zkjsvqFIfD4XBsMv2W6mPlgB1jJc7Ot/md584ReEIv0Su2TTXEqebCUo/3Hd6+KfI6xeFwOBybzHBL9ddmmzx/dplmLybVsNY6ElFmSE3KoenqHZcVXIzD4XA4Np1L9R4jJZ+FdsjzZ5cRoOh7ayqNPpmG33vu/J0ScQXO4nA4HI5Noh/XOH65yWuzLTJt1xIvBwptFMTZup/1gS++scyfvnyJDzy8684JzSZbHCLymyIyJyLH1nlfRORXROS0iLwsIu+80zI6HA7HRjC8VOwje0ZphSlvLHbwxBCnBhGh4K3/eVtRbvj4s3fe6thsV9V/AI5e5/2/DjyQ/30Y+LU7IJPD4XBsOMNxjW0jZR6/f4JKwSNKbVxjslpgz5YKnlz7WcFu4ylhthneYck3f+nYp0Vk/3U2+T7gt/J1xp8TkXER2WGMmbkjAjocDscGcaneI/DguTNNmmHCaCngXXvHeWWmzTt2j3GlGbLUSdg2WqTeiemlVyMenrJKw1PC9OidaTMyzN0e49gFXBh6fjF/7RrFISIfxlol7N27944I53A4HG+Woic8e2aJkZLPSNEnTDJOzkYc2TnCfVM1CoHH4/eXOXpkmqeOzfLMqVm+drkFQOAJSgQR4Ucfv/Pj3d2uONYw0tZONDDGPAE8AfDoo49eLxnB4XA4Np28G/oK4lRzbqlHpRiwa9wqjX41+PmlLrVSwGtzHVphSsFXfPg9++94YBzufsVxEdgz9Hw3cHmTZHE4HI7bRpwZ3n3fFs4sdGmHKUrsin5hkq1YQ7zfSuTD77mPp47NMjVavkap3GnudsXxJPATIvIJ4DGg4eIbDofj64Fd42UavYTHD0wC8LlXZ2lHGYlO+OOXLlL0FSOlgN9+9hz/4gce5vCOsU1TFKvZVMUhIr8PvBfYKiIXgY8AAYAx5teBTwLfA5wGusDf2RxJHQ6H49ZY3Xvq6JFpzsy3+fiz55lthgQCnSSjWggoBcLrcx27JKyviBJNL87wlfDM6UVOzDTuGqUBm59V9bdu8L4B/v4dEsfhcDhuC6t7TzV6CR/5f17h9GwLUYpUW8UAhsoWj0v1BA0ESvBECNOMNDNcrPfYM17hqWOzd5Xi2Ow6DofD4fi6Y7hGQ4nkPahatOMMX0GWaTwFIDTChO1jJfZOlOklmnaUorVBKUgyQ5hmHJ9pbPZXWoFTHA6Hw3Gb6feeGqYTZ/n64Mqu7pcZ0syw1EnwRCgHHiqvzTCA1rZ9er2X8MrlJifuIuXhFIfD4XDcZnaNl2mF6YrXlAgCJJkerCduAGPg9fkOr860SFLbQl1h1xov+gpPBCXwxNNn7xrl4RSHw+Fw3GaOHpmm0Uto9BLmWj0+f3IOEet6Wu7EGMOgZboSSLMMg8H3BGM0UaYJPKFS9NlaK7JrvMJYOeCpY7Ob/dUApzgcDofjttOvu0jSjC+cXgLgPQ9sZawckGjIhkqUs1yJAIwUfYqBT7XoM1oKmKoVUUo4uK3KSMnnUr23Cd/mWpzicDgcjtvMiZkGv/PsOT51fJZWLyHVhslakR3jZYqedT15wqD7rQECT2EwdOOMXpyx2InR2vCufeNsrZVohSm7xsub+r363O0FgA6Hw3FP0K/beOVyg9NzbdJMk2aagq+4VO+x3I7oJhptrLlR8BVJdnVZ2CjVhKlGCZR8hed7dOIMbczA7fWhd+/erK+3AmdxOBwOx1tkeG2NZi+hHaa0o4yCpxAET2Cxm5BkGSKCr4Qo1QMXFeTuK2NjHmFqeGj7CCMln69eaDJWDgatR+4GnMXhcDgcb5Hhuo12ZAPdBV9sHUaSEiUaEfBEYQwYYyh4QpIaNFc7twa+4CtF0VckBt5zaIqZRshPv//QZn69a3CKw+FwONZhrbYhq2f9J2Ya/PnxK2BgpBxYBaGENDN0opRy4NFLbApVnGlqBUUv7ffFNQOlIdjK8WLgsX20SDtM76q4xjDOVeVwOBxrMOx+Gu5WO1xL0d+m4CkCT4iSjE6U4onQjlIwhl6SDTrfVgs+vRR2jRUpeIIROwgrrNXRiTVhnHK5HrLYjji32OHokenNOgXr4iwOh8PhWINh9xMw+D/cN+qjn3mNF8/X6cQpmTZM1YoUA0Wzp8m0tSU8gbGST5YX9EWpphFmJNowWvIJPEUzTBEgzTRhpvE8YbRcQMlaSxJtPk5xOByOW+ZmXDh3wz7fCpfqPXaMrVyWdbiW4k9fvsTnX1ugHChGij69JGOmEaLEUCkGHN4xwvmlHoEn7MzdTZeWumRas9yN81iHRklAKc+w8pSACH/j4Z1MjVgr525rcAjOVeVwOG6Rm3Hh3A37fKus1TZkOObw8WfPUwk8fKVQSqgWfXxPUErxwPQID+8eJ/AUmYbFdkSYZESZpuB7jJd8PAVRCu0oYWqkSKXoUyv5bK0VmBqxCutuKvobxikOh8NxS6zV+fWttsPYiH2+VYbbhmhjeGOhzXNnFnnlcoNf/vQpLi512TZaINW2WaExtjFhqg0Hp6qIQNkXmmHClWbExeUeCiHwFLsmKuyfrOKJrRq/uNylE6V044xHdl+1Llxw3OFwfF2wVufXtzoz3oh9vlX6bUPGygGvXmly8kqbQ9tqHN4xSqOXEKWaXqLZMVbC84Q40ygljBR9ROCvTi/SjjW1okc5UMSZoR1n9OKUV6+0eGOxQ5rXbmQGyoHdLvDUiqI/FxxfhYgcBf4N4AG/YYz5hVXv/23gXwGX8pd+1RjzG3dUSIfjFtksX/2dOm5/ydN+sBje+sx4I/Z5O+gv1/rLnz41aDQINlB+ZOcIX77QoDjusWusRCvK8JSwb6LCsUtNOlGKrwSDx7YRnwvLIbE2RJlBuFq7oYEDExU+8PBOzs63mWlGBL7HrvEyH3r37rsuvgGbqDhExAP+LfB+4CLwJRF50hhzfNWmf2CM+Yk7LqDD8SZYa+W3J54+uyFVv8OKouAJs82IPROVDT/u0SPTPPH0WcBaBa0wXbcdxs0qs1vZ51rc7HHW2+5Gn18rUP7wni2EicaIMNsMmR4t8ZPfcT8AH3nyOI1eQtFXTI8UWWzHRENl4kM9DhFgoR0CsG9rlULg8Ut/85Gb+t6bxWZaHH8NOG2MOQMgIp8Avg9YrTgcjnuGm0nhvB2sVlBPn5qnFaZsHyuiJNiw48JVF87wQLvWzHgtJfqLT51k51iJKDMrBuib3efNnIv1lOZ6273v8BSfOTF/3c+vtojmWyGvXG6SZIb3v316hQJ64umz7BovowQyDbPNkHas15QdbKPDXmLfvxusrJthMxXHLuDC0POLwGNrbPdfi8h7gFPATxtjLqyxDSLyYeDDAHv37r3NojocN8eNUjhvF6sVVJIZakWP03MdttZuf0bOWjPyG7XBWC1jnGacX+yy1Il5z6Gpawbo/t+tcjPK+sRMg488eZyldsxErcDBqeogc+njz57n7TtGr/v5YYsoTFK+dHYZAzx2YMuK79GX5aGdozR7CcvdZKAU1iPJbKuRu62R4fXYzOD4WpUtZtXzPwH2G2MeBj4DfHy9nRljnjDGPGqMeXRqauo2iulw3DzDKZzzrZBnzyzyZ1+7wvml7m1NLV0dTK6VfARohsngtds1e32zqbKrZTw936FW9GwQ+TZmTt0osN6Xf7EdsaXiEyUZL52vM98KGSn5zDbDG37+qWOztKOE4zNNa91FCQbDmfkuSZYNvkdflqmREo/fP8meifI1g9ow/ViHiCFOs7uqkeH12EyL4yKwZ+j5buDy8AbGmMWhp/8X8C/vgFzfEGxkILW/FsGXLzQwGL55zzg/8vi+u+aG2Mjv3p+ZLrUjTs22QcBTsGO0eFtjDqtdJwenqjx/ZolayWe22ePETIvlbsK3H5zkxEzjLR3zzbrfVsvYDlN8BaOlqwHw/gD9Vq7JjQLrffm31oqESUYpsItgnJ7vUPA9pkftWhdrfX7YvfW27aOcX+zwtYt1Sr6iqRMa3YSTV5psqQSUCwHfvGdssK+pkRLf8bYSc82I5W4ECFlmyFbJXyt6jJULnJ7v3NT3vRvYTIvjS8ADInKfiBSAHwKeHN5ARHYMPf0gcOIOyvd1y0YWW52YafBLnzrFs2eW8BUUlPD8mSV+8amTd8V6yRtdaNb31V9pRbSjhFaYkGSamWZElunbVpewusag4Hvsnaywa6zEX71uV5z71oMTBL73lr/fm02VXS1j4AntKOPgtupgm1aYUvDkLV2T1cdZncbal//gtqpd8yLJKHjCUjum0Uv40cf3rvv51fUlZ+Y7ZNrk62YInTijHabMNCJmGj2+8PoCxy83Vuxr32QFpRTGGIqBGrhaFFDwhO1jJXaMFVlsR5te9HizbJrFYYxJReQngE9h03F/0xjzioj8M+AFY8yTwD8QkQ8CKbAE/O3NkvfridsVwO3PEo/PNGj0UkZLPs0wpd6JGSn5g5kdIix14pva/0anlG5U8Hq13GIMtWJAMbAtsqMk49Rsm26yer755lgrmPyPjj7IU8dm2TZWXjF7hrf2/d5Kqmw5UDx/dhFB2D9RJqoZAs9DGzPInKoE6pavyerz/b7DU5ya7awZWO/Lv7VW4l37xjk91+HScheNrdo+NdtZ8fmCJ1QCxceeeYPjl5s8smcUCJhvhZyebwN20SVjYrIMUgNpaus1BGGuFXEwzZhppBQ94eC2GmcXOrTDhEQbvLzZYeAJmYGJqrWEttaKA5fX3WKdr8em1nEYYz4JfHLVaz839PhngZ+903J9PXC9Afh6AdxbSWt84umzaK05v9hFRGh2E9pRSqMXs3eygp0P2MZurTC94Qz1TqSybkTwei25T893GCl6jFXs4FoKPKJU0xxqYfFWleRaweSPPfPGNd8vSlM+fXzpllNQ+7yZVNnhc/K+w9ODz3xgjQH+Y8+8wUTt5i2a4X37Cj5/co7//GXrlvuxb9u/IiDeX5Hv4nKPB6dr7J2sMj2Scrne4917x9k7WaXRS/jMiXk+/J77ODPf5qOffZ1UGyaqAUmW8fyZZd62PeHMQteu3mcM2sCqbiR04ozlboIALxd9fvy9BwbZWvdtrQ6u/YGtFc4sdLlc71ErerYVSao5smt004sebxbX5PDrkBsNwOvNIIu5y2CtdMXhm33YhD8x06QUeJQCewMsdFJSbXh9rsNkrcCWSgFPCQVf3XCGeivWwJsddPvfPckyTs91aIYJBU9xZOfodc/n9Y61ltyVgqIVpYwl2aAjqjGGsbJ/U9foVhkeJF+bbXFk1yhbayUW2iHPn1lmpOTfcgpqnzeTKts/J3Ga8fzZJu0wJfCELwSKf/EDD19zTd5YaHOlGdEME0ZLAdtHi+zfWrvhvr9yoUHRV4yXfY5dbg6+AzA4v4d3jFIJPE5esRZfo5fyzr3jg/33r9tHP/saX3h9iV6cogSWuz5FX6EQXjxXZ/tokZLv0VitMXL62VMCnFvs8q8+dZJHdo8zVg44uK3GS+frACy0Yx6crjHfilAilAJvcL0avcSl4zo2lvUGtBsNwOvNIMtruAyWOxEf/ezrfMv9kysGmHaU8LbtozTDhJGi/RmlmSZMNAKEqWa+FbHQjhgr+jy0e/ya1gmrXV3nFzvsnajwwHTtmpTS21XsdvTINL/0qVOcXbAZPgUltMOUy41wzSDyzQzwa1kxO8bKzDRCSoE3GAz3T1YGg9XtdJkNy/jI7jG+dHaZZ19f4rEDWzh+uUWcZGQFj8+emKNW8tk+UrypFNRhbjVV9lK9h68YDOy1okeUZDxzevGa83xousofv3SRatFnpOjR6CVcrvf4roeu/l6Gr/8rlxs8snuMMwtdir6iFHgYY6yiHsrSGj6/903VmMhdQWtdryhNefq1BaIkI/AVGGhHGXGq2TZSpNWJiTJN0b8ao7hetlSgoNFLmWmE7N9aY2qkxDv3jnN6rs1sK+Lx+7fyXQ9ND5T3SMm/p9JxneK4R7negDZ8Y1i/bIdWLwFhoFzWmkF+7Jk3mFzlMjg912a2GfLFszZj5+BUdXDztcKU0VIwyFRZaMcUPEWCoeQDApk2RNrwg+/atWYxltaaU1datKOMdpTSi1vMNUO+9YGtbK2VOLfQ4exCh7/32y+xpWJnol+50KATpcw1Qx7eMzZQMjcz6B7eMcb0aJGFdmRrH0o+D+0ao+B7a37+Zgb4tSy4HWMlmr2UwztGVyjn4YDtW3WZ9QfTTx+fJfCEh3aOMj1a5rEDwiuXmzx/ZonFTkyaabr5YkIzDXh9ro0IPHbfxE0ff3jgLnqCAeJVRXzD7Bov8/mTc4OBHQARSr7iI08eZ+9EhYInCPDlCw2UwGI7ZKZhqBV9Hpiqcmq2wwe49rf+2mwrr6MwbK0VARtzGC0FK77Deud3ret1/HKLONVoA1Fi25srgTjVXKqH1Eo+Rd+jUvTwujJYa2MtDOB5donYxU48eH1qpETB93i8HAxqYA5M1d5U0eNm4xTHPcp6A9pvP3uO80tdvnK+TilQtMOEsUqBwBNE5JqCq2HWqo49t9SlUvAGM8aXztd5x54xRvMZ0vbRIq9eaQ2yVVS+bOb+rTUqBR9jDEudZDAIrJb/hTeWaORuskqg6CYZS134yvk690/V+PKFOuXAY7zs04sznj2zhDEG34MzC23OL3fZO1HhoZ2jdsW1myDODO85NMViO+L0fIevXGhQK3qMloNrrLhXLjc4vGOlG2vlqE00AAAgAElEQVR1jv98K+SZ04tsqQS8fecIRd9HKcVPfuf9NwzYXi/gvNrKEhhUXB+arg5mq9poMIqXztd5595xpkZK7BiLOLvQIc5dZFGa4CnFaMknyQxRmnF+sbPCHbRewHt44A48ePbMEgK8+74ta1pgJ2YaLLRCTs+1qRY8pkaL+ErR6CVgDEvtmANbKzx/ZgmDLajrxhmIsHu8jO8p268pzy5a/Vs/smuUZ19fohNljBTt56JUs2vcVtDHmWakFBAnGfdNXfv9+hb3cidiphGy2Im50uiRDPWQyrKrikHQvGu6ykwjZqkdk2ozWLFvPRq9hFKgmG9p/vTlGSaqAdtHS3ieWmFRvNmix83GKY57lLVmrGGS8levL/HNe8ZodGOuNEJSbQh8hacU79o3RuCtPbOGq26cL7cjojRjuWMDfePlApL7YgFOzLR474PbBrGOvt84SjPq3YTd4yUqBfvTilLNRDVYMZMdXqP5jcUO5cDD9xTaAImmF6ccn2lxer5NwVO0lLBrS5l2nFLwFJ0oIc6Tk6oFxVwzZLmTMFYJ+Jk//OoNYx67xsucnW9zaq49cKM0w5S5VsQvPnWSfZNVAs8GXd9Y7HJuscO7908MKo37A9CfvnxpEEgdKXl045QvnF7iyM4RxisF/uLkArvGyysCtsPn+noB59UB4OfPLBGlttDsy+eXqXcTvmnXKHsnKoyVC4RJhjGa/+/UPIEnzDRCyoFH0Rdaobb1JGi6sY25jJZ8Ts622VIt3jDgPTxwP3emOUjNPbPQ5fEDk4NthltuaK2pFDzqvYR6mLB3S4WRom9XvSsHnFnoUsv3M9cKKfoevoJ6L2H3lopdJa+Xrvlb31or8diBLTz3+hLLvZQtlYD7Jsu8eqU9UGhXGiGfPznP82cX8ZRCCdRKAd//zTt46tgsM/Uur811qBQ8xso+YazJG9Vew2jJZ7GT8rbtNZ7rxijJuI7BAdj26klmV/gr+sJiJ6YZpvzkd9x/TyqK1TjFcY+y1oz1xEyLLZWAWsm3A26UAoZ6V/iuh6bZWisx1+qtmWHTRxt7Rwi2TXSt6BNn1poo+jYXfbYVsdAK+dgzb7BrvMz//P5Dg0Hjp//gq6SafJariVLNvonKoJjqt589xzOnF2mHydVtEk2SaeLMINh1DUTAGPCV0E005xd7eMq2ntYIJq8L8D2hFaYYk1EtejcV8zh6ZJqf/sQVkKtLeQIUPcVSJ2bneIkvn2/kDeoKXK6HfPJrV9hSDQiUteJEbJCzFHjsmSjjKUWm4MBkgTPzXb7l/jKTNX9dWW4UcB4erI/PNPE9YbGT0Y1SSgWfZpjw3JlFKgWPg9uq/NXpRZY6MUpZyyRODVmWUvAVnspn0RqUGKZGigR5HcP/+/Jlksywd0uZ/+G9B9bMtjo+02CyWuD4TJOTV5pUCz5bKgFJPisftsCeOjaL1pqTs23GKwHaWHflUiemWvSoFgIOTlUHVt4wxjCwXDF2wF7vt170fT74jl2Dyct/evEC9W4MCJ97dQ6dWUtrvhXje0K54DE9UuB3nrvAO/eOow1MjRSJkozFdsT1moLsnaygRFjsJIyW7f0gxtBJ9ECBCNjeVMY+9hTsnahQCjyKgcf7DkzS6F1red+rOMVxj7LWjHW5m/DQjhovnqtT9BWTVTsTjTN7W6yXYdMf1J46Nsu+ySo7xkqcnu8wn8cBqkUZBHnTTFMJPALfQ5Py+ZNz/MEL55msFtkxVuK+yQqvXmmxlNl0xn0TFTxPcWi6yhNPn+X45QbtMKEZJmhtb7DUQCvKCJTgKckL2uyAZ5WXRzvKkAzCNCPNDEqg4Pv0EjtT3DlewlOKpU7E6bkOC+2Ijzx5nP/tg29fM0toz2SZRjehFdk4zZFdo7x4bpk41Zye62CMYabeY7kbk2oIPE2YJPlACIFn5QyTjAvLPfZOVDDG8ML5ZQTF8Znmin5Iq/smDQ/MP/Zt+wfb9JXxsIusHaa0wgTB0MlNLWMMYQafemWWnWMlGx9KMkSETOuBy8Vmc4HvgRLFSMknTjVL7YTEwL4ttiVGO8r445ds44bV2VZffmOJZpShjSHTVp7FTsTu3K01XGX96eOzzDZ7+EqxbbTIjjGPpU5MK0zJDByatoHiWqlDlNe0jJYDagWfei8GA8XAY99EZeBmWuu3fm6xw86xEh975g2uNLpcqYd4nhB4sNyx18wTKPg24SNKNeeXQyZrRa40o/y6+5xuhrSG2rSspuDJwOXaDJPcYvM4uM3Kdmq2ZSdBgKcUI4HCV0IjTGnlmWTt8FoFe68jxtzA5roHefTRR80LL7yw2WJsOKsHoIVWyLHLTcDWDXSilIvLPTwl3D9lq2bbYcpjB666XfozuZ9+/yF+5g+/uiITJs1sYDDNNN/7jh0UfZ/nXl/kwe01aiWfF8/VrQXSDPGUouBBlNrBzWaeaEqBT63o047SQcWtByjPBhiT1CD5TK2gQJSgNYyVPcBaE2OVgE6UEiUaIxAooeh7iAiVgnURAQSeolb0KfiKZi+hHWXs2lLmJ7/jfj7w8K4V5+zPj1+xabh5GiRY1xRAK0xYbMeEqbWEht0SAgQeaHPVIvI9RaXgkWmb2TNVK1At+ix1EsbKPltrRUbLAb/xo+9e4YLqD4IXlro0ejFJZohTTcFXNLoxD+8eZ6Ts8/mT88y3QrTpz2yFKFt53ypgpOxhMkM71ngqlzF/LzP2vD2wvUacauZbMdtGimypFgDsLB87aA9nW7022+TPj8+Szz1WuHICBdvHylQLPt//zh28eqXDmfk25xe7BL6QaavQlYh1j5UDjIGFdkSjG1MPU2oFj2/aNcrJ2c7AzVQKrKW2Om6yVlZdlKb80YsXB65LuUZGYbziI9jB/OFdo7TjjNFSwHIn5uximzRbP15R9oVK0acT2QNoowmUohB4lPMUdG0ME1V7zS8s2fstSjM8EcI0oxx4HNo+yvaRIvdN1W7YHHKzEJEXjTGP3sy2zuK4x7jeegK//ew5G5AsekzViqTaYIAkyzg522aiGvDYfVeVBqycBV2TCRN4TI1omr2Ur15o8v63T7NnsszeySpfPLtE0VcstCMbd4hTmqG1BNLMDEz2ME0HvmoRO5glQEUJngipWMXhCxixr4nSaGOX2Bwt25s+0zBZK3JwqsrJuRaNbooIJGlGueCTakOWaRbaEUoJgVKUAsVCK+Jn//gYXzi9yLcenBzMplenrRZ9O8BrY5hp9IjSDJMXeg1jzyeDJIDMGEQb2lGa++mt1bSUx4d6SUYzTGmG6eDarU5q+IvZFldaEZXC1RXgljoJnz4xS9H3KHqGVDMky7XDnAaMhsSA51ml1rfeVO732z1Rph1ldKKUSsFjvDLs+lE0w4R6L+Gx+yZYaIecnuvw8sX6QAGtJtHQ7MVUAo9/9/mzbBspcHCqxsXlnrUmBeaaIVuqRfZNVKiVfGYadt2JcsEn8BVJpkkNPH5gYpCpNVYOrptddHq+w45Rm1r7uRNLDBfjr5Yz1Yb59lWL4lK9x56JKge3VXnqWBuj11caHvY32wrtpMdT9tpnmSbwFb0kI9OGcqAYKQWkmbauXg3jJZ/ZVjQI/r860+TVGXho18g6R7u3cIrjLuJGhWY3s57A/skKV5oR55a6eCLWJ5sKcWZo9lKu1HtMj17NnBnOpDl6ZJr//OVLjJf9QfwhSTVTI4U8XgJTtSKtMB3Ub0SpRgkDyyFOzeBGvOb/0B3azYulPLHxFE9ZX/hUrchSJ6EVppR8xY7xEqm2N/Aju8c4t9Rj13iF6RHN2cUOYaIJfE2SGVKtwRiUUohvB82Sr/CUTU996VydB7fXiNNsUAXcjTOeP7PE9z6yiw99t50J/o+/+xLL3cTGWrh2YDGA5AOzyd1mqTZ4SnhwusbJ2TaZNhQ8RS+x58RX8GMffwEM7N5SHtSqLLRDZhohmTa0jWG5azN7yoG1yIxRLHUzSr4QJmZdX7xVUjaN1BchNVaeiu+jgCjLeM+hbYyUfJ4+Nc+Z+Q6n59ooZa0BXwm9xAZ9n/zKJRq9BMmzlYb9+AIDGUq+ohNnXKh36caaZi/h4nKPyWqBZpgOUloPbavhebb+Yd9klYd3jw/kHrZ412N1Vtdrsy1OXG6wb7LDxeXuup8b/u31udyI2DleYqJaxGirtNY7n75v75vRUkCWryUOmlQbeonmgW1VUg1hnHJk5yjPvL7Ini1lRIQrzdBe/1QTZ+B7htGiz39+aYb3HNp2zwfIneK4S7iZQrPfefYcZ+bbxJnNWT+4zdZUDBdzvWPvOC+eq7PQisi0phVZO3yyGrDQjvnC60tcaoS8a99Vl0A/k+bwjjG+/eAkxy43aUWpnf2LkGqYHinS6CXMNiMbg/AUYZLhiU0R1dgUxlt1fNpAoqEUBIyVbdbN1GiReifGAL7n8e79YxjglcvNFdaQJzY4bqvSi5xd7NCLDaLtzE8wxKkgorhU71HyPU7PtcmMnWFPjRRZ7sbMtiJeuWxTP48emeaDj+zik1+7zIWlnlVGa2Cw/u80E5QSKnl2Vr2bEHgKJZoozcgMjJY8otQQpnYp0fl2RDNMede+cb56oUGSW4Y6dz8ZoJtYZVMt+gSeGjTHu1gP1zyHStnW3COlgKV2jFL2O9p2LwmlwOP5s4t5rYNHmmVEKYxXAtqhdettrRV4YKrKC+frYOzMf1jZr846ilIbX/KG5A4TTStMmajaFPDA87hvqsbRI9N87Jk30KQ8d6bJXCskSjVFTygXgutmwfWttLlmj+fPLtOLM0QMl+s9unF20785hZ2AnFvsMV5p0ks1vljrcfU6S2Nl316H2CrAwLMWWT+tHeBKM+LA1ip7Jsr8/A88zC9/+tRAET75lUt0IhvjCJRQLXhEmbVM74VeVDdiM7vjOoZY3YVz9VoFJ2Ya/OXpRYwxjBR9wiTjxXN1ojQdrCcw3wp5ba5DnGZ04pR2ZDuPbq0VqPdSSoEi8OByvceffHWGZ16bJ4pTPvbMG/zyp09xYqbBDz++jwNTNR67b5JKwWaEADwwXWOsHLBnosLOsRJHdo5S76WMVwqMlXyun9W+PqmBOIN2mHCl0WO+GVkraEuZR/dP8Mmf+nb+i4OTfO1ig1cuN3ltrsV8s0eYZKTaUPI94lSz3EvwxBamYWw6ZJzZmIuvDFobFjsRJ2daLHUiMm2tjYVWTEHZmEhfWR+arrJ3okop8HKL6Cr9xxnQiTKqRZ9dY2W2VAJmGj0u1UPCJCPLLZ+Rok8rylAi1Io+E9UCzdx195Xzdd5Y7FxjmfXJNLTzwac/SFeDa2/ZfsZUmkGUZBixcRCtDVGSEaaaakENfjcnZtpMjRQpFzzrOtRQK/pUiz6vzbdJM+saW8r7Lq21cM6wvPFQvMUG2tP83Cn+z//m4YE18ZXzy/zRi5f48rllLi136UY2oSPT+rpdYS/Ve4RJyvNnlxGgVrSp29Y9uJ501yJi/5phwuV6iO/Zz6bGWr6D7QCdaeLU+rGWOrGdBGhDL9ZorRkp5rVNqebScsgvf/oUh6arg+rvKNWk+bojlaKPn1tc2piviwC5Uxx3CTdqXf3UsVm2VAJErEthoR0x2wz57Ik5Ros+5xY6vHS+TpR32awWfESEqVqRbpLhKyHJrM++4HvsHCux0I45OdvGVwwGzTPzbSqB4otnlzg528ITeNe+8RUtQKLM8PM/8DD//kfeyaP7tlhfej5ju9nbWIb+FzwhNRBlkGhDsxvz6uUWX3htnm//hc/yD//wZS43QqoFRZpq3ljscXquha8UcZrRSzRZpu1gLflA2k+LFOjEmlaUkmWGFBuDmWmEXK73MMD2sSLtKBso61OzHX7muw9xZMfIigBzybcDjwKKuaXTiRLC1CqQfuqpNgZjhDTT+MpWHxtgS8X27hopeYyWfC7WexS869+CYWIr6gVswD5v6bIWnhISbSh6+cJOSgh8qzBsfZ11f3WilNlWRLXo811vn2bfZIVd4yVmmxHNXoqXz8IxNnus7F97XYefX6PwjI1Z7ZksD+Jvv/SpU8w2QowxZMYqpn7CROCp6y7oZDO2WsRpRphmdBNrTdo4kxoM+p6sVACryQyk2rpVZ5u9QaV4/2/4+0TaWh1KbAfbei9F55ahMdBNMhbaNlvskT2jg0aJ7zs8ZS21/HpYN6AizezvwlNyT/SiuhHOVXWX0C9Ku9KKaIfpoKdQPyXxUr3H4R0jPH9mieVuQsETAk/oxRnFQPG1Sw2KgTeoSxgp+URpxmwzAqwro5e3BqkUPJv5JHaA6RdyDfel+s7D2wbrWC934nUbAl5pRgTKulKiRF83H75PoOwgFmc24KhXRaD7boM0NVzIXTMFz6aV9rN7olQT5FlcIyWPMNU2liDg+YpMG+s6ywdAT6xv3lc2UKyxg/L9U1U8pajkltWwsp4aK3NwW5XZRkg7ylZ0QxUxZMa6LeZa0SBNM9MapRQHtlY5u9CmHdssqYlqQDWf8U+Pljm8Y5R2lKGwazZkq0ffnEpB0Yk0cZaRDQVylWDbu+QnxPeEfRNVlroxSmyG2dZakTcWOyCG+XbEctemuxoxYKz75KXzdXwlzDQjDAbfEzKsFWGwsZssX4wqy2NNYGM86zWI97AtN94+VJOy0I7wfcV4wRvIIVh3WmbMuh18wboPP/Gl88SpRhQoBINNEa8VA3ZvqTBT79GKUuuizE/SejawAYw2VAs+y2kyOJ/D10BrQ5gYSoFHolO0sZZJraAoBh7Nnu0/9vj9E4NJFcCp2c7AwvrqhWW+dqlJL7GZVZVCQNH3runZdi/iFMddwupGb81ewky9x3fnP7J+EVS1aOMS7ThDa0Ot5LFzrMxSJ2a05A9iE7WSTyuMWerEg+mhwua1b6kUuFTv5ZXFinY+Is7kleZj5YCFtg3YXqrbyulyQVEp2AynuWaP/+qX/oLlTkKcZShsDyJtzCA4HijrAlkL+7q9S0ue0Elu7OZqR3b1OJO7jnSeXhpy1Z1jJ+9C2Ve0ogzPE9J8AFQiVIteXlwo7BkvcXaxw0wjItWaqVqRz56YpZdkTNaK/M6z5xgrBzywzWYJrZYwSSHDzv4zbWgaO5A0ejb9qVLw2DVe5mI95F17x5ltxsw27OBWKfg0ugkHpiqcnG2zc6w0UJDDKIFunKGUYIzYtFoMQV7Q10/NHSv7hKm28Yo4Jcts4WWrFw/cYiIMCvZ8ZS2JSsEHY1OkO/mga4wh0VettcxYJVEp5BlQGnwRfF8GNTSrKRU8lrsr+3JFeVqqzUKyCtagyAx4ItetLzq8Y4xttSJRnNGJMzIMSl1NIf7Z73mQj//VOY5fblLwFWGaESV6XWXsKcHLlW5/k/62/WSIVBvaYYLnic2Uy89FqeBz/1SNuWaPUuCvUBrDk46jR6Y5v9Tl2w5ODtqaBJ5yleOO28up2Q7v3Ds+aC09Vg54cLo2qDTtF0HNt3qDNEAldkGYU7NtDPD2nWPEqe0nlWlNlNptIDfTsR07o6RFqu0ax5nWbBu1P/7FTkzJV3zu1VnOL3XxlZCmdtaeZLb5W5ik+J5HZqAXp8SZreAuF3zKWNcRWOXgK5t5014deRxGFJ7KSG9gqug8FtIf9EAIfI8gzWwO/1AKcJqlg8eBJ4OspyjVlPKB5exCh06cocQOaLONEFHCRCVg+0iRz7w6x1StwMXlnq0fWSVPf7bdd4tlqSGMM0q+tb5aUUq54PPeQ1t5YHqU515f4EozJMmyQcHZucU2vqfYOlKk6F1bm6ENZJk9ypayrU+pdxNEcueZwIPbR0i1LU6LUs1EpcCl5R6C4VIjGQyEvhJSbQPuIsJkrcB8O6LeiYhzF+Zqd402drJRDpS9vsbQDK2C6a6h7PvHSlLbNuZff/oUb98xRsGzdTcjJVttH3gqL0q1GWRJahXwQztHB/G9pbyAc2/edaBatOnDmTHWmjS2w0DgKQ5M1Ti4rUa9mzDXikjT9ZWGldNmSK1VwubnEx6DtTBMltcZaatwRks+33JgclDzM8xwhuJwZ4DA93j8/q23fVGyzWRTFYeIHAX+Dda6/Q1jzC+ser8I/BbwLmAR+JAx5o07Leed4FK9x97J6oqmc8OBtMM7xnjf4SmeOjZDpjUF36PgKevuUHb22OglnJlvk2UZlxshvVjbVcZyJdO/mcLU4GF9/1ES0+jFvLHQGbi9bO2Arc3QQMmDSp4xlBmb2plqOy21tRmGdl4g5XtWyfhytYXEWgxXNvcHtZsh0/Z4vjIYbdtcx5keKEhtrsY3+rk+gefRSzI6UUac2kB1qq0bYvuorSROMk3N96mVbMuWKNVcaUQr5FewpiuuH0BuxzZttlr0WGrHiMCRnaMcPTLNa1eaA2WujR2gMm3IdMblem9FgLnP8CvtKGPcU1SLPr0ko1rwMMD9U1VeOl/ngakqs01rXY6WferdhDhPkS54tj7Gxl+sldcMbeV0mrGu0raWmnWDJZkeuBSz622PVdLL3ZjPnZjj6VNzBErhe9Z9iDHE2g7+JV84tH2ERi/h3fvGB/VF862QU7NtEq05sLWS9wzrWEWhhGLgDeqF0szwkSePM1b2OTBVsbUT68Q5+r8JraGbx1dWv59qOxgNTwx8sfUxWptBEkW/5qfRS9bt9XWvNjC8GTZNcYiIB/xb4P3AReBLIvKkMeb40GY/BiwbYw6KyA8B/xL40J2XduO5XrfUfn3Hn3z1kr359XCBFzS6CdMjRSqB4vRc+6rfW+U+6vz+8IaUhxm8D2lqiwRLeZwiSjW1gk+Yz7SjDOJOQsbVmy9Or6Zg2pmbzmf3Cl80pYJCSe4OMNfmyw9cBPrWUngN+bqCBjpJRpL26yiuvq9yd5bvwSD/Y+CKuBqWL/pCL9WUCx418VD5+Tw932GqVmC2GWGMjQ/1C8WGz2EfERssT7TBV3Yd6oLn5a1SUn7xqZN88ewSWf/CGQZFa0pAaUPg2QF59b69ftGkNvls3coepZqpEbvY0Xc9NM0XTi9yftm2ui/mMZ7xvK9SP8A/UvLpxjZFuADMNiPi6yhshY1BxKkm0WZFjGU9dP8c5dtmKcRoAqXxfRskLgce20aL7N5S4R8dfZDffvYcr1xu8vLFJrWSTzdKQWxfsn4Xg+mRImcXOlSLAVtKHosdu83OsSJL7Zi5Zo/FvHNtMBT7Wf3bEXv6iLKr1vjw+7AydqMNeUacylOqtS1QzGt+7sWW6LeDzbQ4/hpw2hhzBkBEPgF8HzCsOL4P+F/zx/8J+FUREfN12Cfl0HSVj37udZJMM1ktsGOshFKKd+8fH3QbnW9FAz9+qjOaoe2+KcpmfQS+x1g5oBtneeWzzSDpI8PVbEM3jRIG9RqVgkcryqgPRYJXuzCG//cfaw3i2dmtCny2j5Zs7n0rtO6D+GrB3+pK5FoxAGyO+41cVn6uB2J9VTAf64YblleA0VKBOLNpkRqoFhX7J2tcaYZ045Qw1ShlrZYkX4SqVvJphynlwGPvRIUrzR71rh4Uvq3lAvHEDuyjpQClhMlqwS6hG6Z86dwyI0WfVGsEQXOtokw1jBQV20ZLeWHiVetGm7yfVz4Qx5ltVKi04W3bRzh6ZJoz821eOl9HiXBga5VGN6ETZQPLI8yLNPuNKrUx+Xm5/m1UDARtNGFqz7snrFswN0zf4hsm1qC0rUspFzzevX+Cgu8Nuh2cW+zadNm6dYcVfNs+Zku1YBsF+opqMcJXwuVGxGgpYGqkiKeEWtljuRPRS6xVdL2fkCdQChSdeGUrmWFW3yYGYaQUMC62cPVSvcdTx2Zts8y7tH3I/8/emwfZdd33nZ/fXd/+em80lgaxEBAgmhRpSBRkiZElSmYUlxUro7EnkUdZplhJysmMUhlHKc9kKlPllOK4SqlRMpVhnKQ0sid24rFixZIpkZY1tCyIFsVNIJoAsZBYekEv6Lff9+5y5o9z3+vXjW4QENBoEDifKlR3v76477zb793fOef3+32/m81WBo4dwIW+ny8Cj250jFIqEpEKMAwsrD2ZiDwJPAkwOTm5GeO9abqyIC9dWEYQHt5V5jNHdwNaWG5b0ePU5QZT1Rpn5hs8+dh9nJpr9Cxac55DGCc0OhGOZZH3nHRvWncsl7MuGVc35rm2pWv6+56/vwNYsbLloHcQ9Gw5CK+/oapL96ZqWxa7h/NMDmU5u9CkGkSMFTO6v6DTWnXj7foZWJZeOXUTvY6oVTen/g9xF0tEq+imP693M3NtQSkdCFsdvSXm2TYXrzRpdOJeorjZ6fSCmedYDGYd2mFMLYg4um+IhyjzZ28ssFBvr5ursVgJhJ1Ip8vrQQPQSdU4Ubr4QEGwZhbczcFEsaId6W5kz7Z6PRvd15Ek2v+he51dGyYHs5y53ODv/NaLzFb1dqbv2FxpdKi1IxyBy1HMWNGn0Y50xVECHzk4wkDe44+Oz+LGCe2NSqNYUQGw04mFZYGDELxddOfqiQXov4EItDox3zu9wPv3DfP9s4vYYjFe8lmst6m0o171W5go5mttPMfCEmHXYI5OnNBeTpgcyrLcCrlcbSMiNNbZelrvb+U4FiIWvqN6uZ2Nxi6Qan7pz8ZivcP2gewtsfp9p7OVgWO9nci1f8brOUY/qNRTwFOgRQ5vbmi3nqmZCr/+9EnOLzYp+Hp/+tjZJWarbbaVfJIk4c2lFnGi8NNp9X/88wvcP17sWbSOFj0uLLV0JUyoiJMI24Kf2DnA5HAegJFiBsey0hr5pCcbIbJS5x7qvORVvN0NoVvfHycre/3d6hsr7XLeO5Jn90gex7I4NVennHNBhFLG4vR8Q/cTQE/qe1vRZzmIermOTCrTkChdGironEz3ubQUkOolMTfCSfs5sq5Nu6OVd8O02iheJwkNurluarbOuyeKdAo6N1LMODy0a4Dvnl4kTkJsW3djN9oRQXpj1b0oVpUUGMYAACAASURBVNo5Hacd07JSEqwU7T6ToC4Kndz1HVINqZh8Wr3UnbULKl11SHrN9DNeaWqlVtcW2qFKVxQr664QcETvwQ/lPXYMarOrf/apB5maqfCHr85c01OiV1HVHaeAY2v5ljhJVl379YL7eqj09ViW9qf41vFZGmFCOeswVszguTZDln799U7EoKdXaperbQbzHntHcpyZbxApxeuzdRKlSNL3eFeX7e1ejy16JdWOdOFGcI0S8rRyGxGYr7cZyvkcuW+wl8CHH8/q925gKwPHRWBX3887gekNjrkoIg5QBpZuz/BuLU8fn2Op0aGQcXqGSCLCQl27kAmKSjPEcyy8VJl2sR5S8AN2DOR6ap4iYFmCZwl5z6GYcdg3kuethQaztTbz1aCX+Hz39jIT5QzPn1uimHHSc3aI2levKtb70AlphqD3oRNdR5+WdirAd0R3Kac3txMzVU7PN1AoHtk9wC8d3c3Z+Tpf+KOTjBR8Ltfavb38gYxNPUzIOLohrpjRhk/byn5qDJWs2quWdQKX0xvP6tcQxIoBR28xBJ2YrO/QakcktlyViLZE7/9nfIePvGuMctblwHieLx87z1w1oOA7ZBwLEe0qVwsirU2VXhfHEnYMZFhodJCO7haOk6SnnltPK5B8R3pBsHt97xvO04kSgijWzY+tENsCz9L5pjDRB9qW7lnRr1kfh0hvVbne5m2C9oHvVh79Urq6BRgr+pye39gxUfr+ObZeTYaxzh+IrA4VG+V+1qILNHS5dLMTEyjI+zZRrLh4pUkrTHr5KT+tNruw1CRRCfcNZzk5q6sHP3z/MMfOXqGZSrXb8vZ5su7YhvMeUaKVBGJZ2Trd6P9kXYtt5QxLjQ4fvH94Q4HQe42tDBw/AO4XkT3AJeAXgb+65pivAZ8FjgH/DfDtd2p+49Ky7lTt7w7vKpKKaMc2z5EVGYREkaiEcwt1okSxreTzRr0NCDnPZSivjXQ8W5fPVgM9q7bSRq3ZakAn1h3TD+8aQATeXGzhuTHjvs18rZMm/q7+wPdXDyXom5IFjBY9Gp04tfnUvRq2ZRElOteSc20u19v8/MM7e5UmTz13jqxrsXckx/HpGjlP+0E4llDv6CY+bItteR/bshjIuni2Vohtp6W9vQ93OvvrnyEmyer8hqBv0IjugZip6ITx9pJPA66KkEL6+pWi0Y5pRxH/9ZXLRAkM5lzeNV7g9GVduus7WkgQVLryWTlJNYhotSOtGYXeF1eo3irOSQNe/y3XtWHvaJ6XLlT4wL5hbEv4k5PztDoWtugViufYhFG0qtclTpPlnq236Taa8bu2Pm6x0eHJx/YA8He/8gJ/dnaJJG31F7qVcSvXTqUL/aIrhDFpgj3BgtTESsg5kLaI9EqSu+8dWec9ZYv23XAsod6OsEQr5O4cyvLmQiMVU1RYttVzrXzPrjIjeZc35hscO6vVmI/sHuT+8RKvXqqy3Fpp3oO3D1xDWa0Z1om1O2EQxRs2MXb58IFRxstZTsxUybirb5cbWe3eC2yZ5IhSKgJ+GfgmMAX8J6XUayLyv4vIz6WH/TtgWEROA/8A+PzWjPbm2TGQxetzmwO95+s7Ng/v0o5kurZc0epomQnb0t2xB8YKzFbaWu3UsylmbO1HAfiuRbOjb2atVEJBBx/dGfv4oXGanYjvnl6klrrudcIEx9Y16TlPlzbaomdXfhq8ujej7mfRd4R3TZT5tZ9/gA/dP4pnW3QzDRnXYryUoRKECNLT2wrjmLPzdb7+oxlem64C2qCpmHHIuFpvKOta7BrM6uYqV3szVIKIgZzL/WNFhvMeOd/GtfWse+321NptBoW+jiPpzFLQKrnT1bbOB7G+fEa9HeM7uhGt0dZltbOVFn9yap4rzQ7lrEM10C6DtmWt2gdvR7riKUz0DLo7IXfSKi1Ib/DpykH6nvnkbJ3JgQwvnV/mD1+dYanewRZFo6MlVEYKHmKt/ph2b/Ld5HZ3ZbD2NQ3nfXKew2jR5+x8nb/7Wy/yzNRl2mGEa1t0ZZ5s0Ss3x+rqTul/iYJyzmXXYFbLZiSkDYCKVrTy/uh/bgV4juDZQsbRAn+eDRPlDJNDOfaMaA2wjGuzcyjLT+woY6XWriohlQHROZ0/PbXA2YUm7909yPZylm1p8cBCPeg9r92NfGte+9oZsSNQzvnplqie8IwVs2Qdq/f+9x19Dbpn82xhvKyrHT97dLKnQ9Utw620wruiC/zHYUv7OJRS3wC+seaxf9L3fQB8+naPazN44oFxXr24zPnFJii9tK63Y/aM5Pmlo7tZbnT43tklWqmNatbVRkXby1n2jBYYKvicmKn2VHC/f3YRO01mNzpapdaxrd5eu1JwbrHBbx17k0oQoZTeA7dEUuE/LUEymPOIlV66W6L1lGqtiE68MhdzLPjJ3YPsGsr1dJxmqwHnFxtkPe39YVta/mRyKAdot8EfvrWMl66gKi3tn9HV/2mnM/pS1uUjh8Z7/+fYmSXiRDGUc+nEilwqvufawpn5Rl/H9OoZZndLTY/XYqHe6T1Q9B2qQQQk6Yphhe4puoq5kr7exUYHx7Kw0NdY0gRxztPCdt1AW+8kvZyNAM1O0lMMdpT0tvsS9E3ZlpXnzHsO4yWPt5ZaLDe1SKNlQS2ISYCCp/tPLNGdy2klb2/cilS+Ja3W6hfyVcBsRVchDW1z+dIfn6HWJ3XeCrW+VrMTY4nq6Tg5FuwoZ1hshNQ7EaWsy3DBp5BxOTVbQylFzvdoRZ3e5MIS7StSb0dp4BGc1BPEd4T3TJbxHYeFepv5WtBbSTbbET+6WCGKdaG3ZetrGkYJjU5MJYh4ZLLMfSMFZqpt2qEO7qcvN7RcvqSFEYlaVTGo0jeEn+Y+UPrvu1Bvk/X01pht6ybIwZybih62sMTCdUWXH8da4r/fG2TvaOGeLb9dy9sGDhH5ZeC3lVJXbsN47loOTZR7Nevdqqqje4f4wP5hvnLsLabmakSJrqpRSuG7NnnfYddQlm9PzTFdadFsR5ycrVLKuDQ7EXnPZrkV9UyFFIqlZgdJq47iWJfpdm823VLO7vaM1knqpDNzGC9lOLpvmO+dXqTe0R7fTiopPTVbJ+87XEzfBdsHMjSCEM/RCqvNjg5K1SDi2NlFmu2ol+T3bIvA0jeHVpiQT6Uucp6Na1u9Jqrjl6oIenXWSUXhWh0tbKc1jfTM0E1lIJqdZNV2keeI1lBSWhpjIGOnqy+dD2qGV+/pd2eXFnoF2OxELDW6joLa2zwRAAdLYLSYWQk+omVAbEvwXZ2X8hztSOhYuhEvinV5qaVWl/M6Ao1OyJ+/eYVSxmEw72qtL6VoJTE2WqCw2YkJ4wTPEVqh6gW27sJgvOhzudHB5mqxPl3Oqzi32Oy9P5z0NQnahMqxIC36Sl+ScLneYVvJx2sJnUhRS8tjo0SR9x1816KrzajS58y4Fs2O4FpayK8TaW+QD90/TD7j8bc+eB+/dewt/vT0IjsHctSCMN0uDdI+kXQASm/ndu1Xz843ObgtYH/a6JgkCRev6JJlBGy1Mg43vS6JSu19Hb1qCCItUZKkvShRovjg3iH2jRV5fbbKSOxrOfh2SDvUDbXv3l7kX3z6oat84u/VQLGW61lxbEM3570I/Hvgm+/UPMNWc2iizD/71IO9n/srrUq+g2sJl2sdQG8R3D9WYGqmxpVmSKIS4rTjtZHuEddaEeNlH9/RN2wLPZ1N0hXNetUiidINdAr94fJsYc9IjlY6c16od9g2kNXy060o9f7WN4LvnVliOK9F5d61rUTOtTk5V6eUcagEEQ9PDjBTCai2QqaXm4wXfey0MqmUcZiu6JtEKSOMFT2CSPH3PrKPU3MNnQOKE967R1etHDuzyJVmiGtrUbui7zJf66SBQ9+1rL7E5lDOpRXGurELwVKKVhTTSGfsviP4tkUUp4UBfe/grGczOehzZmFFrVaxWi58qdHBtmBmuYVSKyZGYqUVYqnyaSEVMsx7Nkp1n1dWEtvSLUFeua6VVshoMUMpI2TcHPUg5OJyi1YnJuvZqERYDpLeuLo32Ywj5HyXHY5NKetyeq66Khci6EKKbmOpYwnKElqRQkRvWXZrD7w0GHe3wbKew3v3DPHKhSqHt5fYMZDlzHy9J2PeLRvu0i3dLWVc9o8XCVJBzfFyjnLW5dBEmZFipld88OzUHO2oRbW1uqopTHQ/T9bVBRONdsS3XpsjShSNIKKTKDxb2DWYJehEOI5FxtPVUe0oYTRns1CPek1+ri1sK3nUUtWAAddipOCx2AgZCyIOT5R54oHxdcvkTZDYmLcNHEqp/0VE/lfg48DfQDfh/Sfg3ymlzmz2AO9mupVWji0sNDppzsOinHWwxeLSckCjrc1/qoHOrDqW3k7JujbLrZBGEDExkKUVNgljhViCrGN52k+MLq+0027nZifmoZ1lxkpZylld2dQMwp7lK0CcJKnvd8wfvjrNRCnDg7vKvH/vMCdmqrx/7zDlrMuOQW05Ole1qAYxHz00zOn5Bu0wZjTWK6JGR9fc/+TkAHtHtR4XQDHj0mzH7Bkt9LZROnFC1rd5dO9Q2rgXYydqxXuDlRtk3nOwJCHrWmmDn0p1kvRNPknNiVqdkHbfHv1gzqES6IR9rFa2k/q3wvT2mNDoRNw/VmC22qYZJmTSfoMkVmQ9G88Wsq5NlCj8NEDkfIcoibDSgC5ipds4upPaSftY2lFCEMZcaXYQhOG8y0De5UcXq+v+HaNEi1B+9F3j7BktcGm5RZiEeismXQVpD4mYgm+R81yuRGGq16VLfhVQ8HWVX6JUT5pkphLwk7uH+NjhlSa3P56a4/xiExG98oi6zyUQRF1Tr4TTl+t4NhQzHt9+/TIf3D/M1EyFS8stJsq6KkmlOa0g9bmwZUU0Ua+4ddnzclP1gnaUdJP5ivl6B9vSTodhrBPvlkCtHeN7FoW0x6nWjhjK+xQyCUsNLRXi2bqasSsRsnZCZ3h7ris5nq4wZtN/ETAI/J6I/Pomju2u59Jyi2oQMl9rE4QxQRiz1GxzYanFpeUmp+ZqLDY6BOmesGN1t6RAUsG1ZpigEHYP59k1mE33ybVk9rWwLSHvO4wWfEYKPmcXmrSjqLd/24oSRgoujU7EfK2t/QjQs2VL4MKVFt87vbjKSApgpJDh/XuHeeKBcXzXxnNs9o7kWKy1WWp22F7OsHsoh2MJL1+s8Df+ww/4zsnLOBZMlHxeurDMufk6UZKwreRTzjgM5TxevlAh59lp+a0Qphv6viMMZJ1e4nsg67AcaE2trnOblVabgQ4SiBYizHoWOVcr6S6ns3LfsXHXmAPlXAvXEQ5uK7JrKI/nOvyHv/FePnpwFNvWjoTlrItvC7V2zORwLl3RaRdAx7IYyLk4tqTbSXom7Ls2Y0Uf37VphjF7hrPMVdu0OjE7BjI8uHOAC0sBnrP+x1SXylrsHsn3HkuU3u/v5oK6r6QaxKkch/RkOSZKPqWs7uTOuHbayKh1pKJYXZX8/Tsf3kvWt1PZGy0dkvUsDk0UGSn47BrMMVLwQCmutCIUig/sG6LRjvjc777CD99c4rlT88zXAgoZrbmVJFp/qpBxer1GTlr62+wkaRk0WGJhWWmlWKS3RsNUe0sn+4VGENHqxNhiUQlCau0OYZRwdqFBFCse3TOI79pcaUYMF/x7tnnvVnA9OY6/jy6JXQB+E/iflVKhiFjAG8CvbO4Q7152DGT501N6JtwOdd1/kkpLNNoxxYyjK0y6DU4Kkkj7QodxSM7VFSHv2zNEEEb84NwVhvMeCsWVRnjN5y74LpFSDKWSDgAnpmt8+OAYTzwwztdevsRCvUPOtUmUIg5VT5LctbUjQrMTc2K6xngpQy2ICON4jW9HsZc0Hch75H0HK715eo7FlUZHdwJHCRevtHh0zyCPTA4wU21jiUUzjEEk9c62aLZtmh2LnQNZxNLbR55t89iBEUaLGd6Yq/InJ+eJYoVY+tbZ7ESMlTJMlDK8udik0tJquI2ObvzIe7oIoRNFZDw79eDWwbm7bSgC5Yxu+Ko0O7y12OwJ6/3sgxMIMFdrc/FKi4PjBSaH81Rbl5itBri2he3qSrWJcpbZaoDv6Eo2R7SulWNrc67XZ+vsHtaNet1+gT85Ob+u7pIel5Y3qQWRthwuZbTlrtKz8jhRvVxDKZUsX2x0GC74vS3Cb/xomvl6B9+WXuBLFBRtePzQaO/GOjVT4dRcg7FUMypOFIN5j88enaTZUat01r5/dpFKK6SUdbEt4eRcHdAVZ7Ug4vmzS+wfy6fvd4XnOStCnAm0ooRYJWRchwTtG48IUbLSgyPolQlAO1FcXGohIhQzju55CVWvVDlOdAluOecyVtJVUiZo3BzXk+MYAT6llHqr/0GlVCIiP7s5w7o36BrUtMIYIW2qUyuNbpKWKNbb0arEp9420RIVA1mbThTzysUKtiUM5T3OX2n29RNcjZVWN8WJ3tJQSvVmiV3p53dtK3Ls7JL2w1b0bnRxsmKJ2QpjrjRD/tETe/j9F6c5t9Cg4Nt4llAPIlph0jvfP/zPr2j9LRGOnV0kjBOaaf4h62r71+fPXeHj7x5jcijH3/rgfXzud17BdumZU2U8h/fuztGMFJNDOXxbmK4EaYJe8fzZpbS3BSxLlwuHcUIYJRzZM8zkcI7vnlnERij6eobdjnRNv+9YLNY72mUwTaJaQMazyLgOpYzL9HKQysanKqlNHVgtS4vw7RzU+/nztSDtibDJpOYZjSDGL9gcHC9wcq6ucxuBFqcsZlzeta3ISxcqHJooXtVkVm+v36jn2sLRvcNU0n6GB3aWeHOxoV0CJe04t3Vw2TWU5/17hwHt9nhqrsETD4xz/FKF5WYntTvVW2merUt5f++Hl9ibGol19dI6UcLOwRwoODBe4PXZBvV2yEjB4/tnq70V9HjRox5EnL7c6Hmfz9faFDIOby02ma609HvF0dt/3a27oJOkKwu7J8eiczqrg2ecdMUI9Ts8AWx0QUUnVNi2hShdkj4xkMWxhFcuVPnY4fF7uhrqVnE9OY5/co3fTd3a4dx9dJVt1zqbdR/PuhZLvX6DNFikH/hWJ6aQcai2wlUfGx049PbJUCHDaDFDOePwxuW69q2IErKORWhpCev+jZeMo70zChmXd28vcXq+QT3QS/0P7h/ufaDyGVcL6C00OTtfp9WJ8GxL6wjZVs/x7oP7h/lLD+7gD1+ZodGOWGx0yPs2D+0oM17O9iQZ+tV/F2oBM8utVC9Ksdzs9Kpfnj4+x/3jRUBbj1aaIbV2RCnj8sCOEkN5n5lKwG98+qGrrm893ct3U7l5Eb3/vpzW3P+9j97PQM7jtemqzgcJNFJ/7G5Py48uVdOSU4eRgsvlqjbIWm52qLcjwkiR9yxt7mNZzFbbHJoo8efnlvjooTFAq+uWc3qFNVfTku1Zz0ormHTQu3SlSTvSMuvvGdVy+m8uNpmaqTFeWmkqu38sz3QluGoSYIvOn3S1zvQ1iPjQ/SOcmK7i2DZDeZeZSkDOs9k/trKd1e14PjRR5lOPbOfsfJ3lZtgrrVUI87U2jSDiK8fe0u+vVC+t24MRhDGztTaHJ0qcnqtx5nKDYsah6Dss1TucXwq4byRHNQgp+g5Xmh2qQdo/Yuki5eF8hm0THi9eqKAUREpRzrlsH8im0iaKtxZ1fuxaYjjd97eu7tPhxLMtMqm74EM7ywwX9PvmXhUlvNUYI6dNZGqmwlPPndPbCH3CaI8fGuXZqXnKWZf37x3mG6/O0AxjPFtvkYSp6J0tWjJCe3qrXsULopffGddithrwX1+5xGJDm/bkPItWRyuzdpO7nqP9GETAdWwGc462n3VsHt0z1PMS6Jek6N7oj+4dZv9onu+dXmCpGZJxLYbzHr6z0oMyNVPhtZka20p6v74dJZxbbFLOuVy80uSLz5zitekKF6+02Fb0WGx0aEWql2BrhQmtMPUejxMmSj5PPXeO0YLPjoHcKqn5SivcsFtXawsJnmNTAIIwIUKvqrpbE790dHfvb3JiukKcKK0vZQm1dsTOwSwPbC/xa2my9OuvXuLffOcsF660dHOkpwPvTCVgW9mnGujyUYXqbRnVg6iXcN6XzthbnYjpSsCOAZeMa3O51qaYdRkr6p6J+4HD24v82emlVR4PQ4UM94/leXOh0VMOdiwh61kcnij2An3/DLo/mCLCRMlf5VTX3/F8aq7BXzg4xldfukiQetMrpZPdYZzw9R/NMJDzGC/6LNTbjBR8QK8O5qsBJ5Ti7EID19EBQRd36NxYkPYJVYOIpUZIwbeZqbS04KRjoVA0I8Vo0U8VhWMKvt46VEpRa0cM5z06caCVntdpDe8GDdta8QLx09Jh37XYUc4yWsxc831juHFM4NhEnj4+Rznr9m583a9fPna+18hXzrqUMg71trbEVKwox8YKqql8afcDokiNZpR2tNMli7rU0LOtnodGf1WVdgLU8iQ5Swu8/YUDw70KqvWambqOg5Amc1P9rI7orbOje4d7JYtffOYUgzn92kSklzP54VtXsC2LnYM5Dk3o8t0/O7NIJ4zTGbNubGulpapRrMh60hN+FKzeNkwQRr3S5A+lVTrAqsA8UvCYrmgLVr39pAPxh/aPrLrBdp3Z5mptHKG3VZN1tYfGn55eZGqmwqGJMn/pwR2cmmtgpboWC/U2UZzgiLBQ63DfSJ5aoGVdumMt+HbacAgP7ND+7MfOLBHFusorSLfoir7bk50B8B2HD+3Xf5fXpitUg4hy1mEw77OjnKEVaae/UsZlW8lfZfrVT7ffoKvG/N3Ti7y52OTw9iK+46wyHOpWOnVSFz5LJHV8TIhS//Txok81iKi09IphIOf1VhC+q7fjhnIui40OUZIwUsywZyTH+SstyjmXakuvaGtBmJpuCZ6tt67COGGk4LPUCBkqeLTTMl793tbyM5VAV7KdXWgQJwmVVtR7f3dXSdk+p8k9IzkuLgW0o4Sdg5lel3e/yZLh5jCB4xbTP9s7MV3loV0lYGXGXMw4XFjS3eP1dowlsNBo4zortf7r5SX6H+vq5PV3mCul7TA7sfZ9yLi6tl7oSnELgzmXVpjQDkM+kG4xbUT3BvuVY2/xvTNLDOZcPvnw9t6Np7/O/dJyi8Pbi7x0Xt/MfUc3Mc5V2/z0wdFewNwzWuDETA1QDLt6m6dfv8myhJ2DWYIw5vXZGpPDOf7Bxw6sGsNP7R/CdeyeBlZ/YH7swCjfOj5LM4xphpBxbN61rcjff/z+q17boYkyC7WAr/9oFtuyeuZHl2sdxkv+KtXT/tdnCT37VJXAoW2FXrIV9GShlHWpBhEHxwskSnFiuqabAm2YrbaZGMjyvvsGObfY7AlSdm9u3fOcX2qyczBHMeNwfrHBi+eXeXjXAO/rWyFeS+6if7X7gX1DTM3U+LPTS3xo/3DvOb74zClOTFd5Y66GY0nPnZC+LvSsZ7N/rMCL55cpZRzma1rGfKkRUkor6SbKemsp57v4rs3+0TyvTVcRhMMTZf7ye7bza994nSjRPTiebZH17J5p2EQ5Q7UVsa3oc2qu3nPZy/s2dpoDqwYRgi4acazVQpWuLbriD/BcC7CYHM7RjmLOLwUc3FY2eY1bjAkct5C1W1NvzNV4/uwVju4TRgoZ5msBP3zrCkuNDp1Il5tOLwe9jm4n7dReJSO+5jk8S3p5C9cSPrB/hNPzDdxGm9lqG9e2KGUsOpHCc5K0+kQL1QVphc17dpV7XubX4tBEmdG+pq1+nj4+1/v62nQFz9ZChouNsFdVNVLwenLvXYbyLnPVFnGiV2DtKOnZzhYyDnlfP0870jPLa43h+XOLPH5o5eY5Usjw+LvHVzWtXcvnuVsy6zndlki9FxhGMV97+RLfOjGLoCU69o7k2TuS4/lzVyD1vQbhR9M1nnxs+Koto6mZSq9TejDn8uGDozTbMS9dWGai5DM5nMe1V6Tn+6UtvvjMqVUBsbuymKm28Vz7uuQu+le75bS66filKi9eWEah3f92DeV4aFeJ589qOQAtJqlfm4hO7u8Y0Fs9j0wOcPpynSBKKGXdtGnR5/5xPbauvMx8LWCx1kYBj+4dpNIKeXZqnl0D2XSloFiot1Mfd5X24Fj8/CMT/PHrCyw02jTaUVoFpvWppiva3S9IbROzro3naFFKnRPUvRy+Y/EXDoxyYJte5SVKmbzGJmECxy1k7dbUu7eXeP7sEscvVTk0ofjBuSssNjqMFz2qQcz5pSaNdtxzevPs1d7b3aanVYiWVtAzdO1xMZxzWay1mShlqLZD6m1t5PRT+4Z5fbamk9+Oxb7RAvvH8gzl/euWg+5v2upSzDicmKlwfqlJOevy0M4yPzh3hddn6zy6d7C3Ksm5Vm/fv8u2UoaTc3VdhmlbWqW0E/f6DlQqILdQ77Dc7OhZ8UyFd6U3g/4xCHLV+X3H4WOHx3uy6P/l5UuMlzJ89ujkVSusTqwT1Qv1Dp1YpS5zwsxygO857M44KHTZ7/NvtillHIZyDrFS+Ap2DmaIEvjqSzM8dmBsVdHDpeUW55eaPLyrfNWW0ky1jevY7Bkt8Hd+et9VAWC9az45nMd17F5RwI383Vbphil4bbpKLYjYVvYZK2Y5uk/4wbkrXFpukvcd8p7ds799cKce22gxg+fYHN03wuc+doAvPnOKNxfqvfJr1xIt25HmTx7YUVqVV7m0LBwYLzBbaxMlehJjIWwfzPZyfocnSqBUb5tv/1iB0WKGoYJPJ9JBNwhj3ShrW2wrCcvNDq3UzvX+0XwvaMC9rV672ZjAcQtZ+4EfLWZ4755BXrlY4ZWLlbSaRGmpCoF2FK+otUpXSn3lfMLVPgdd0Tbfsch7Nu1OxKlam8G8y8HxIpfrbS4stjgwXmD3SJ5Ly3qvdyjnUQ1CTl9usK0Ubbg/vpaNvNAr6FILcwAAIABJREFUraiXuC5nXR7dK7w2XV1V8gj08iTdZK9tW7xrvEC9EzFXbYOCgbxL1rVphUna0RtR9G3GSnp/+sJii5xrrxqzziuUe3mF7vkrrZCJsscX/ugked9hrOBRbYV84Y9OAqwKHjsGsnTCmFjVeyWjZ+cbJAhjBZ+Mpz8ewyIEnZjFZkd3NdvCWClDznNQSsuWd1dg/SvOl88vU2lqD5buTXT3SB7PvXYAuJb//PXSf45uSSxAMev0kvenLzcYKWQYKWT4mQe28fpslcMTZS4tt64qde5e2+7f9cB4nt9/8SJ536Hoa7vhVhizfzTPe/cMY/Up1hYzjlZDsC0OT5Qo9m23dfNN3fdRLV1ttKNEe78XMxQzDjOViI8f3nbVden+3M3J9RcWmLzG5rFlsup3IzsGstSC1TX3zXaMI8KFpSbnFhramlXpMtC0cnDFiGfNvlSYrA4aJV8nIm1L5yvGShk6iaKQcTg4XuTXPvUg//a/fy9f/MWH2DNaYKYSsL2cIU4UYaIoeDaVVsiL55c5ML56C2kjnnhgfF05aW28tDLvGC1meOzAKIe3l/jcxw708ghPPraHclaXhZazLk8+tof37xvhvfcN89ce3c1fe3Q3Hz+8Dc/W2yIjBZ/Ror5p7x8rpKZKuvdh7Rg+c3T3Ved//NAoXz52nkorpBZov+1y1iXvO3z52PmrXpttWxwYK+A7FksNnbwdKXgM5PpXMRZiwVgxw2DOZ/dwnlwaVNqR9ojv+lB3b4CWCEMFDxGt5trlegLARtf8RiS8+89RaXUgLabYP5rXXdrQS8p3x3V4osznPnaA3/j0Q/zapx7kV544eNXfrrs6OjXX4OFdA5SyLo2O3r56eNcASboK7Kd77vXeC4fSQNV9L5Uybk96p56ep3vNrnVdNnqvmbzG5mBWHLeQ/kqkYsbhrYUGf/7mErYlhJH2VwgjCKL1u7r744Znr3a1K2VsihmX+ZpOpA8XfT7yLn0j6e7ldulX8fziM6dwbYvZWpt6oGWyD4wVrivH0T1Xd1bYX4H19PG565oVb6Qo2n+dPMdmcjjH9nKG755ZZLzo97YpAHK+TRBG/H+nLlMLIooZh6NpM1v/+bs5plorouBrHajp5YDtAxmKvs1cNVg1hv7X5rk2R/dlWagFHJ+u6obDtDqs65vywPYSP3xruZfQ1n0YCfcN59gxkNWVUC29VVjIOAznXKrNkIV6e91Z+41e8xu5CfafwxLdef3IZLl3TZ8/u0QhVSbYaFzXUoO9tNxi90iePaMrq8BEKZphvO4qsDv+9c7XvzraP5bnh28t9yoG+yui+gs2np1a7AkSXs94DbcWEzhuIWs/8LO1djqDiunEb29v2f/7MF5Rf3Ut/XWkoOWfh/Mu/SoU15rFbvQBvxHLy+u5+d/I1sB6N8ZfeeJgLzHcH5DmawE/OHcFz9YGSQXPptqK+NM3FnjxrWX+3kf39bafujP+YtahEyb4rt4uWWp0KGZcxkuZdceytgeiq1i81jflM0d384H9w3zp22e0dEfe477hHJZlcWA8zzdfmwXQQSWMObcYMl7yaEU6sN9IALgVN8HuOboTmu62U3+gvtFxddloO+3d2/Xz3UjQ659wDeX9Xnd9Kbu6aKBLK0x4dM/wKpfJG11dbNSYa7g+ZCsU0kVkCPhd4D7gTeC/Xc/vQ0Ri4Efpj+eVUj+39pj1OHLkiHrhhRduzWBvgn/4n1/hO6/PsdgIrxk0+q1au9ii9ZzaYaK9D0R4ZHKQ4bzL67N1ihmHxw6MrtorXu+Nv/ZGDCv7wjdbbbIZH77+yrRixuG5U/PUAu09EkQxS2mjo+sI5awHwBd/4aFVsian52r82ZnF1IIWGh29XfX5v3jwmiXI/WO4lsz2eq/76eNzvLlQ5+TcSr6kmhpdfPEXH9rym9Kt/lut/Tu93fvwVo3vVryfb/XY7xZE5IdKqSPXc+xWrTg+D/yxUuoLIvL59Od/tM5xLaXUe27v0G4NUzMVXj5/hYW3ERuEq4NGVybcsYTEFu1GVvR6NfztSDFe8q9rtrh2++xWJg03Y2tg7WokjBWP7h3k5QuVtIZfS210Yr2V0Z+YPr/U5OXzywwVPH5ie4lzS01qrYhi1rnuoNEdw7Vkttd73f/uu28yOazzB91Ko1LGoZxz74ib0a3+Wx2aKPfySXPVoFe59uM+x/WOb6MqvxtZQW/UmNvfu2O4NlsVOD4JfDj9/svAd1g/cLwj6c5olurtDYUGQa80LNGBI+27WuUxEcWJrrRCdyN3q15+4WcOXPcb/Fbsl99u1uZodDLeZbYSkE09Jrr+7cN5j9emdWnwRMmn0uxQbYWg4MjkILZt3ZaZZHfrplulBFw1M76bmJqp9Epou7I1z07Ns3e0sKnX+lZUnN2K4HOvs1WBY1wpNQOglJoRkbENjsuIyAtoD5AvKKX+y0YnFJEngScBJicnb/V4r4vucvuZE3O4tm7Uc6wVx7Z+bMCxtXGQb2tzpq4lqZsKvHXSKDKQdfjQgbFVFSQ3wjs5adhdMW0r+Zy+XNOe3pb0mgfvG84xU233lGm7M/6FepvZWpt/+nOHb8tr38yV3Z3I9czaN2Mr81Zc51sRfO51Nq0cV0SeFZHj6/z75A2cZjLdc/urwL8UkX0bHaiUekopdUQpdWR0dPSmx3+jdFcZlVZIsxMyU2nR7MREiV5VrL3QYyVd8rmtlGW04Gm/jfR3cSpl7ooOIg9sX9G16m7L3Ct0V0z3jRS4f6yAnQr8DeRcDo4XsCztmLjWSOoTPzHB5FDutgXMe60ctL+Etkv/rL3/89Av8NnVGPtxuRXX+VaUO9/rbNqKQyn1+Ea/E5E5EZlIVxsTwOUNzjGdfj0rIt8BHgbuSLva7gwsjGNqgW7QyLsWtY42aLLo2l6SWoL6DOdcTszUWG6FWKLd7JJYG9AotAjg9nKGVLvtnl1Or6yYDmyYmL4TZpDv5JXdjfJ2s/bNzCPc7HV+J27f3mls1VbV19Cugl9Iv/7B2gNEZBBoKqXaIjIC/BRwx1rVdvdN//xclZGCx0K9Q9ZzaEWdXjf4QMZhuOADWub8xGyVINKBBdGy0Xa6OhHRxxR856pGqHuZW10abPjxeLstozs9j3AvBfnNYKs6x78AfExE3gA+lv6MiBwRkd9MjzkEvCAirwB/gs5xnNiS0V4H3a7xaqBnYdsHMohA3ncYzHmMFlzGS1kc28K2oBaE1FoxnTBJpUh0r0Y3HxInutP39OU67TA0y+lrcK9tE90JvN01X09FwUx87h62pI9js9mKPo7unu7Z+TpKKUQk1SnSSqJXmh0Gci47BnLsH8vzzdfmCKOERifCtrRF5lqfGsfSHuSeI/yVR3au6iUwGO5kTK/EO48b6eMwWlW3iO4M7IHtJZZbEc1OBKLVbnOeTcF3aLZjLtcCTl9uEMUJviNkXQdLBNuSVRavjqUNkcpZB9e2UWA+cIZ3DGYVeHdjJEdukrXJ2s8c3c3kcJYvffsMQepmtn8sz2w1ptXWW1kzywFxkhBYuvt5pJBlvt5muRESJorBnEM2FdGL4gQPeOnC8ta+UIPhBjF5hLsXEzhugq6uUdeY6Y25Gs+emGWh3sa2hNGCRyeG45eqCPT0pcRSqbmSwnOERGkXtcmhHFMztZ6jX5woogSGctp7wmAwGO4ETOC4Cb5y7C3OLzYppBLj2kOgjmuJttdsaie8JNHlta4NKEEsC99SJIliYiDH5FCuV1r6pWff4Htnl2iF2ha2lLWJE3jfLjNzMxgMdwYmx3ETvHRhmYJvk0l9vzOuTRwrmp2EIIzpRIpa6vCngCgGsYDUH1wB5azD3/rgfYDWOxrMexwYLzA5lKOUdch5Tk+Z1WAwGO4EzIrjJhDkKh2qRGntKVFXe4aLQBQrQhSOCFnf5o25Gr/xzVPsGsoxUc6kfhMu+0s+nVgZyWeDwXDHYQLHTfDwrjLHzi4hIvip6J70RQtZEznibge5AAIDOQ/XFhbqbR7YoQNDOevCUO6WyJ4bDAbDZmC2qm6CzxzdzZ4RbcHateHMuBa+NlzrdYyncQLQPhKeY5PzbD6wbxjHsmhH8arz3kkdtgaDwbAWs+K4CQ5NlPmHP3OAp4/PadvQIGIgp+VGLJXguRZhnBAn2pTJtrSEyLZyBt+1GS1m8JyrY7fpsDUYDHcyJnDcJN3cw/mlJjsHc+wbzfHs1GUWah0SlWALKAHX1uW5l2sd6kHEoYkilVbIUN7DEqHSCo3OksFgeEdgAsctoF8Z9+x8k6LvUG2FqASKWZdyRjf0BWHCfSN5HtxZphMrylmXX3niYO8cRqnTYDC8EzCB4xZwabmFa8NL5yv4jsVIwcexhIvLAR/YO0zOt5maqdGJEx7dM7iu5pQJFAaD4Z2CCRw3QVdu5LXpCov1NoM5j4xrA5D1HHYNZjmz0KAWRAzmXH5qv7bY/NzvvMKu4SyHJ8qm1NZgMLzjMIHjBukPFhevtDg4XuChnWW+/uoMnSjAtQXb0qW5D08O8MrFCh951xjlrMt8LeDU5ToIVJphzxXNiL8ZDIZ3EiZw3ABff/USX/rjM0SJotmJ8B2Lk3N1fnL3AJPDeWYrAbPVNntHC7x7ewnPsRGkZ7F5er6B71j4jkWtHd1SVzSDwWC4XWxJH4eIfFpEXhORREQ21H8XkSdE5KSInBaRz9/OMa5laqbCl759BgSG8i6tMNaexYnilQsVgk5ELQiptkKaQUijraujHt5V7hna1IOo1yhYyuigYXo2DAbDO42tagA8DnwKeG6jA0TEBv418BeBw8B/JyKHb8/wrubp43OEcUIp4yAiOBa0wpi3lpqcnK1xudom61o4tvDmUotXL1Z4/NAonzm6u+feV/BtqkFEO0rYP6YbB03PhsFgeKexJYFDKTWllDr5Noe9DzitlDqrlOoAvwN8cvNHtz6vTVdodWLeuFzn7HyNVpiQKEUnSr/GCQph70iBfaN5SlmXU3ONVYY2pXRr6uB4gaG8b+xgDQbDO5I7OcexA7jQ9/NF4NGtGMjUjE6E+65FEMbU2zGJUriWEIrCsgTb1npVed9BKUU1CHtbUP2GNmuNn0zPhsFgeKexaYFDRJ4Ftq3zq19VSv3B9Zxincc2NEgXkSeBJwEmJyeva4zXy9PH5zg4XuDkXJ3RosX5pRZJkhCLsHckz1KjQ6MdstToUAuWAYVrW5SyHlMzlVWBwbiiGQyGdzqbFjiUUo/f5CkuArv6ft4JTF/j+Z4CngI4cuTIhgHmx+HScovJ4TxhnPDKxQphrEUJu5GtFoREiRY2bEfa5k9EkXXElNsaDIa7jjtZHfcHwP0iskdEPOAXga9txUB8W3j6+AzPvbHAQj1AKW0DG8YJS40OSumlUFcN17WFrGvTirSsyNPH57Zi2AaDwbApbFU57s+LyEXgKPB1Eflm+vh2EfkGgFIqAn4Z+CYwBfwnpdRrt3usUzMVpisBs5U2nSgmikEpsC0QEapBSKLAtcCxwLIgUYp2GFEPIlNuazAY7jq2JDmulPoq8NV1Hp8GPtH38zeAb9zGoV3F08fn2D2c56XzyziWRZwkWAKebVHIOCw1QmxLESuwU+MmBQSRQsSU2xoMhruPO3mr6o7g0nKLYsbBsYVCxiHn2XiuhSI1aBJtIYsC27K0v7jSj4dRYsptDQbDXYcJHG/DjoEstSBiWylDJ9YluHGsUCjaYULRd3Aci1LGRtKViOdY5Dwbx7ZMYtxgMNx13Ml9HHcETzwwzlPPnWPfaJ5aEFJvxzhWgu/aeI7FI5MDnJipknEdShmHdpTQjhIOjhe4b6RggobBYLjrMIHjbeh2fj99fI5mGDOzHFBvh+Q8l4d3lfnM0d2cna/zpW+fYbHRYTjvcd9wDsuyzBaVwWC4KzGBYwPWdngfGM8zXws4OVtnMOdxeHsR17F7fRpf/IWHVh1vfDYMBsPdigkc6zA1U+Gp585RzrpMlDO8uVDn91+8SNa1yTjCbCXg3EKDyaEc+0bzPH18js997IAJFAaD4Z7ABI4+uquMZ07M4drCu7eXsMTlzHyDZidmrhqAAstSgMW5hToL9YBXL1WuWmmsXbGYFYjBYLhbMFVVKd1VRqUVkqgElOLF88ucmq1y/koT24IwVoSJoh1paZFmJ2G+1qHSDJkoZ3qOfl9/9VLvXP2PT81UtvplGgwGw01jVhwpTx+fo5x1KWddHMtittqmEydcvNLEsy3iRPdtJOnxXZ2qRGnpEUuk5+j35WPnOTxR6v1snP4MBsPdhFlxpHQb/eZrAY22NltyLOhECZ0oIQgjEH3BukHDEi01EiUrmorFjMNcNejZxfY/bqRHDAbD3YAJHCndRr/T8w1KWZddQ1kssbAtIePaTJSzuLYOGYLWqgIIE3CsFQX4WhAxXsr07GL7HzfSIwaD4W7ABI6UJx4Yp9IKWap38GzBEmEw7/HTB0fxHOHClWbaMZ6agihwbC1qiILLtVbP0e+zRyd73ydKGac/g8FwV2ECR0q30W+o4LHUDIkTRRQn/OkbC0xX2oSxwnOsXuCwLCHvuYwVPCYGsrxyoUo56/LkY3v4Sw/u6NnFzlSC3uMmv2EwGO4GTHK8j0MTZf7pzx3m158+yRtzdSrNDrWONm1SCoIwwdF6hsSJQkR4z65B9o8XmakEfO5jB1adywQKg8FwN2JWHGs4NFFmezlDsxPRiuLe44qViirfsXBsYedglsVmaPIXBoPhnsIEjjVMzVQ4dnaRK01tB7uWROkqKsey8Gxhqd4x+QuDwXBPsVUOgJ8WkddEJBGRI9c47k0R+ZGIvCwiL2z2uKZmKvz60ydZaoRsZFqeoBsBLYGZasBQwTP5C4PBcE+xVTmO48CngP/rOo79aaXUwiaPB4DfOvYWJ6arxMk6S40+Cr6FIISx4rNHJ03QMBgM9xRbsuJQSk0ppU5uxXNvxNRMhWdfv0ylFdKJNw4cDuDYNlnPppRxODXXuH2DNBgMhjuAOz3HoYBvicgPReTJax0oIk+KyAsi8sL8/PwNP9FXjr1FO0qIk4R2tP5GlWdDKetgiTBa9HAsy3SDGwyGe45N26oSkWeBbev86leVUn9wnaf5KaXUtIiMAc+IyOtKqefWO1Ap9RTwFMCRI0c2SlFsyEsXlin5NpVmqMcPvTxH11u84DlkfYcoVszX2mwrZ001lcFguOfYtMChlHr8FpxjOv16WUS+CrwPWDdw3CyCEIQJri0kiSLuCz2ljE2UQIIQRlrtsNmJGcp7pprKYDDcc9yxW1UikheRYvd74OPopPqm8PCuMs1OjGMJvmuR82x8R8i6Fohw30ienYNZXVUVKXYP5/mVJw6axLjBYLjn2JKqKhH5eeBLwCjwdRF5WSn1MyKyHfhNpdQngHHgqyLSHef/o5R6erPG9Jmju/nemUWWGh0a7ag7TlwHojhh30ie3SN5akFEpRWaElyDwXDPIkrdcDrgjufIkSPqhRduvO3j669e4je+eYq5WhuU7tVwHJsd5QwP7RqgEyvj5mcwGO5KROSHSqkN++r6MVpVKVMzFU7NNVAoMo5FGCuiJMFWiihWCPAbn35oq4dpMBgMW84dm+O4nfTbxooInSgmiGI8WyfDz8zX+YNXpvn6q5e2eqgGg8Gw5ZjAwYptbCeKma+1aYUJcaKoBXGa5xDiRPGlPz5jfMMNBsM9jwkcaNvYIIx48fyy7t9QWswwAeIkJk4UcaK4tNzif/vaCRM8DAbDPY0JHGjb2KmZGr5jUcy6eM6KFWwUAUqwLCj4Nov1Nk89d84ED4PBcM9iAgfaNvZKM0QpxWDOxbUtLMCzQAkg4Ds2xYzLSMGnnHV5+vjcVg/bYDAYtgQTONDmTR/aP4yIkCjYPZxnx2CmJzlSyjiMFzNYlrB/LE8x4xiNKoPBcM9iynFTPnN0N089d45y1qWYcagFEReWmsxVA1phzGDeY/9YnpFChkorNBpVBoPhnsUEjpRDE2UePzTKl4+dZ64aMF7K8Nmjk+wdLawKKJVWSKUV8gvv3bnVQzYYDIYtwWxVpUzNVHh2ap7DEyV+7qHtHJ4o8eyUlmd/8rE9lLMuM5WActY1ciMGg+Gexqw4Uvp7OZ4/V6UeRLi28JVjb/HPPvWgCRQGg8GQYlYcKf29HO0wpuDboBTfPb1oSm8NBoOhDxM4Uvp7OTKujYiACIM5U3prMBgM/ZjAkdLfy6GUIghj2lHC4e1FU3prMBgMfZjAkdLfy1FrR2Rcm5/cPYDvOKb01mAwGPowgaOPzxzdzd7RAo/uGeZ9e4ZwbZtKKzT2sAaDwdDHlgQOEfkXIvK6iLwqIl8VkYENjntCRE6KyGkR+fxmj+vQRNmU3hoMBsPbsFXluM8A/1gpFYnIPwf+MfCP+g8QERv418DHgIvAD0Tka0qpE5s5sEMTZRMoDAaD4RpsyYpDKfUtpVSU/vh9YL027PcBp5VSZ5VSHeB3gE/erjEaDAaDYX3uhAbAvwn87jqP7wAu9P18EXh0o5OIyJPAkwCTk5M3PIipmQpPH5/j0nLL+IobDAbDNdi0FYeIPCsix9f598m+Y34ViIDfXu8U6zym1nlM/0Kpp5RSR5RSR0ZHR29orP3WsRNlLWJoPDcMBoNhfTZtxaGUevxavxeRzwI/C3xUKbVeQLgI7Or7eScwfetGuEJXbqScdQF6X58+PmdWHQaDwbCGraqqegKdDP85pVRzg8N+ANwvIntExAN+EfjaZozn0nKLYmZ1DDWeGwaDwbA+W5Xj+FeADzwjIgDfV0r9bRHZDvymUuoTacXVLwPfBGzg3yulXtuMwewYyPLmQp3ZaptqEFLKuGwr+dw3UtiMpzMYDIZ3NFsSOJRS+zd4fBr4RN/P3wC+sdnjOTCe5/dfvEjedyj6uulvernFx99tGv8MBoNhLaZzHDg11+DhXQOUsi6NTkIp6/LwrgFOzTW2emgGg8Fwx3EnlONuOZeWW+weybNndGVrKlHK5DgMBoNhHcyKA53jqAXRqsdqQWTEDQ0Gg2EdTOBAS6p3vcQTpXrfG3FDg8FguBoTODDihgaDwXAjmBxHihE3NBgMhuvDrDgMBoPBcEOYwGEwGAyGG8IEDoPBYDDcECZwGAwGg+GGMIHDYDAYDDeErK9o/s5GROaBt37M/z4CLNzC4dwqzLhuDDOuG8OM68a5U8f2445rt1LqusyM7srAcTOIyAtKqSNbPY61mHHdGGZcN4YZ141zp47tdozLbFUZDAaD4YYwgcNgMBgMN4QJHFfz1FYPYAPMuG4MM64bw4zrxrlTx7bp4zI5DoPBYDDcEGbFYTAYDIYbwgQOg8FgMNwQ93zgEJF/ISKvi8irIvJVERnY4LgnROSkiJwWkc/fhnF9WkReE5FERDYsrRORN0XkRyLysoi8cAeN63ZfryEReUZE3ki/Dm5wXJxeq5dF5GubOJ5rvn4R8UXkd9PfPy8i923WWG5wXH9dROb7rtH/cJvG9e9F5LKIHN/g9yIi/0c67ldF5JE7ZFwfFpFK3/X6J7dpXLtE5E9EZCr9PP6P6xyzeddMKXVP/wM+Djjp9/8c+OfrHGMDZ4C9gAe8Ahze5HEdAg4C3wGOXOO4N4GR23i93nZcW3S9fh34fPr959f7O6a/q9+Ga/S2rx/4u8C/Sb//ReB375Bx/XXgX92u91Pf8z4GPAIc3+D3nwD+CBDg/cDzd8i4Pgz84RZcrwngkfT7InBqnb/lpl2ze37FoZT6llKq6xv7fWDnOoe9DzitlDqrlOoAvwN8cpPHNaWUOrmZz/HjcJ3juu3XKz3/l9Pvvwz85U1+vmtxPa+/f7y/B3xUROQOGNeWoJR6Dli6xiGfBP5vpfk+MCAiE3fAuLYEpdSMUurF9PsaMAXsWHPYpl2zez5wrOFvoiP0WnYAF/p+vsjVf6StQgHfEpEfisiTWz2YlK24XuNKqRnQHypgbIPjMiLygoh8X0Q2K7hcz+vvHZNOXCrA8CaN50bGBfBX0q2N3xORXZs8puvlTv4MHhWRV0Tkj0Tk3bf7ydNtzoeB59f8atOu2T3hACgizwLb1vnVryql/iA95leBCPjt9U6xzmM3Xcd8PeO6Dn5KKTUtImPAMyLyejpL2spx3fbrdQOnmUyv117g2yLyI6XUmZsd2xqu5/VvyjV6G67nOf8r8B+VUm0R+dvoVdFHNnlc18NWXK/r4UW0xlNdRD4B/Bfg/tv15CJSAP5f4H9SSlXX/nqd/3JLrtk9ETiUUo9f6/ci8lngZ4GPqnRzcA0Xgf6Z105gerPHdZ3nmE6/XhaRr6K3I24qcNyCcd326yUicyIyoZSaSZfjlzc4R/d6nRWR76Bnarc6cFzP6+8ec1FEHKDM5m+JvO24lFKLfT/+W3Te705gU95TN0v/zVop9Q0R+T9FZEQptenihyLiooPGbyulfn+dQzbtmt3zW1Ui8gTw/7d3x6BRBFEAhv/XqCAWiRZiIwpCCgsFEYmClYIpAoK1KUIghY2F2NsYUOw0IGprYZdCFFFbUQRjUAsVMRZiYWEjSIqxmDkIiQmZ092N5P9gud29Pe7tsNzbnX03exEYTSn9XGGzF8C+iNgTEZvINzMbq8hZq4jYGhHbevPkG/1/rP5oWRftNQOMlfkxYNmVUUQMRMTmMr8DOAq8bSCWtez/4njPAE9WOGlpNa4lfeCj5L7z9WAGOFsqhY4AP3pdk12KiJ29e1MRcZj8m/p99U/9k+8N4DbwLqV0bYXNmmuztqsB1tsEfCD3A74qU6/SZRdwf9F2I+TKhY/kLpum4zpNPmP4BXwDHi6Ni1wdM1umN+slro7aazvwGHhfXgfL+kPArTI/DMyV9poDxhuMZ9n+A5fIJygAW4B75fh7Duxtuo3WGNflcizNAk+BoZbiugt8BRbK8TUOTAKT5f0Arpe451il0rC1n/ekAAABAklEQVTluM4taq9nwHBLcR0jdzu9XvTbNdJWmznkiCSpyobvqpIk1TFxSJKqmDgkSVVMHJKkKiYOSVIVE4fUsDKS6aeIGCzLA2V5d9exSf0wcUgNSyl9AaaBqbJqCriZUvrcXVRS//wfh9SCMjzES+AOMAEcTHmEWum/syHGqpK6llJaiIgLwAPgpElD/zO7qqT2nCIPX7G/60Ckv2HikFoQEQeAE+QnsZ1v4yFEUlNMHFLDykim0+RnJswDV4Cr3UYl9c/EITVvAphPKT0qyzeAoYg43mFMUt+sqpIkVfGKQ5JUxcQhSapi4pAkVTFxSJKqmDgkSVVMHJKkKiYOSVKV3/5rk3+qJxrbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(df_train.X,df_train.Y,'o',alpha = 0.5)\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('y')\n",
    "plt.title('Train Data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-1.448637</td>\n",
       "      <td>-0.238349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.956413</td>\n",
       "      <td>-0.027134</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.282732</td>\n",
       "      <td>-0.350744</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-1.700539</td>\n",
       "      <td>-0.731897</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.938974</td>\n",
       "      <td>-1.320550</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          X         Y\n",
       "0 -1.448637 -0.238349\n",
       "1  0.956413 -0.027134\n",
       "2 -0.282732 -0.350744\n",
       "3 -1.700539 -0.731897\n",
       "4 -1.938974 -1.320550"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_test.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Test Data')"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X10XHd95/H3V5rRgy1Z8oMiy7LjOE9gxyeQoCQN0GxKEmoCJYWWJu3SQlvqzbZptzmlbbps027P2ZYCLW2BwrrAKXBYEqBNcYtxSEJpli1Jo4Q4KFZiHCfxkyLLdjSSbI00M/ruH3Mlj+UZaUbzcK+sz+scHc3DT/d+5+rO/d7f073m7oiIiBSrLuwARERkcVHiEBGRkihxiIhISZQ4RESkJEocIiJSEiUOEREpiRKHiIiURIlDliwzG8v5mTKz8Zzn/7mM5T5mZu+d4/3XmpnnrOsVM9tpZj9RwjruNLOHFxqjSDmUOGTJcveW6R/gIPBTOa99ucqrz+Ss+yrgUeBfzOyOKq9XpGxKHCIFmFm9mf2hmR0ws+Nm9mUzaw/eW25m95nZSTMbNrPHzWylmf0FcA3w2aA28RfzrcfdB9z9Y8CfAR/NWf+9ZvaimY2aWZ+ZvT14/Srgr4Abp2sswevvMrM9ZjZiZi+b2X+v/FYRUeIQmcvvAm8F3gysB1LAx4P3PgDEgG5gDXAXMOnuvwM8AXwgqFH8Tgnr+0dgvZltCp4/D7wRaAP+HLjPzNa4+w+A3wa+G6xjbVB+BPgFoB14F/BBM9u2gM8tMiclDpHC/gtwj7sfdfck8D+B283MyCaRDuASd0+7+xPufqrM9R0Nfq8CcPf7g9rIlLt/CTgCvKHQH7v7I+7+bFD+KeCrwH8qMyaRcyhxiOQRJIcNwK6gKWoY+AHZ78xq4HPAvwFfN7PDZvanZlZf5mq7g98ngxh+1cyeyVn/pWRrN4VifpOZ/ZuZDZlZAnj/XOVFFkqJQyQPz142+gjwFndvz/lpcvfj7j7h7ve6+2uBG4D3ANMd2wu95PS7gMPu/qKZXQ58AtgOrHL3dmA/YHOs46vA/cAGd28D/j6nvEjFKHGIFPYZ4MNmtgHAzC4ws58KHt9sZlvMrI5s30IayAR/NwhcXOxKzGytmd0N/AFwT/ByCzAFDAF1ZnYn2RrHtEFgg5nFg2VY8Dcn3D1pZm8km8xEKk6JQ6SwjwAPA98xs1Hg34Grg/e6gW8Ao0AfsIvsGT9kO9B/ycxeNbOPFFh2fTAi6hSwB7gJuG16GHDQR/EZoBcYADYFj6ftBl4CjpnZ4aCGdCfwsSDW3wO+VubnF8nLdCMnEREphWocIiJSEiUOEREpiRKHiIiURIlDRERKEgs7gGpYs2aNX3TRRWGHISKyaDz55JPH3b2jmLKhJY5gbPwXgbVkx6vvcPe/nlXGgL8GbgVOA+8PhinO6aKLLqK3t3e+YiIiEjCzl4stG2aNIw38jrs/ZWatwJNm9pC7780p8zbgsuDnOuDTwW8REQlJaH0cwcXbngoejwL9nLlWz7TbgC961mNAu5l11ThUERHJEYnOcTO7iOzNbB6f9VY3cCjn+WHOTS7Ty9huZr1m1js0NFSNMEVEhAgkDjNrAf4B+G13H5n9dp4/yTvV3d13uHuPu/d0dBTVvyMiIgsQauIILtD2D8CX3f0f8xQ5TPbS1tPWc+aeBSIiEoIwR1UZ2Xsa9Lv7XxYothO4y8zuI9spnnD3gVrFKCKyGPQPJNjdN8iR4XG625vZtrWTzV1tVVtfmDWONwG/CLzFzJ4Ofm41szuDS0hD9oqjB8jeh+DvgF8PKVYRkUjqH0iw49EXSYyn6GprIjGeYsejL9I/kKjaOkOrcbj795jnJjPBpaJ/ozYRiYgsPrv7BmlrjtPWHAeY+b27b7BqtY7QO8dFRGThjgyP09p0dh2gtSnGkeHxqq1TiUNEZBHrbm9mNJk+67XRZJru9uaqrVOJQ0RkEdu2tZPEeIrEeIop95nH27Z2Vm2dShwiIovY5q42tt+wibbmOAOJJG3NcbbfsKmqo6rOy6vjiogsJZu72qqaKGZTjUNEREqiGoeIyCK3lCYAiohImcKYAKjEISKyiOVOAKwzm3m8u2+wautU4hARWcQ0AVBEREqiCYAiIlISTQAUEZGSaAKgiIiUrNYTAJU4REQWoVrP3cgV9q1jP29mx8ysr8D7N5pZIudGT/fWOkYRkagJY+5GrrBrHH8PfBL44hxl/q+7v6M24YiIRN+Xvv8yB4bGSGWclqYYl3Ysn5m7UYtaR6g1Dnd/FDgZZgwiIotJ/0CC7+0/Ae60NNYzkcrw1MFhkql0Vedu5FoMo6quN7M9ZvYtM7si7GBERMK0u2+QlcviYIaZ0RSvpzFWR//AaFXnbuSKeuJ4Ctjo7q8DPgH8U6GCZrbdzHrNrHdoaKhmAYqI1NKR4XG2rGtlIj1FMpXB3XF3Xj1d3bkbuSKdONx9xN3Hgse7gLiZrSlQdoe797h7T0dHR03jFBGple72ZhpjMd6wsZ2meD2jE2nMjB+/dPXSGFU1HzNba2YWPL6WbLwnwo1KRCQ80zPF4/X1XLtpFddtWs3FHS289/qNNYsh1FFVZvYV4EZgjZkdBv4IiAO4+2eAnwX+q5mlgXHgDnf3kMIVEQnd9Ezx3Dkct1+zvqYTAO18PA739PR4b29v2GGIiCwaZvaku/cUUzbSTVUiIhI9ShwiIlISJQ4RESmJEoeIiJREiUNEREqixCEiIiVR4hARkZIocYiISEnCvh+HiIgEwryrXylU4xARiYCw7+pXCiUOEZEI2N03SFtznLbmOHVmM4939w2GHdo5lDhERCLgyPA4rU1n9x60NsVqdle/UihxiIhEQHd7M6PJ9FmvjSbTNburXymUOEREImD6PhuJ8RRT7jOPa3VXv1IocYiIRMD0fTbamuMMJJK0Nce5eXMHu/sG+eDX9vDxh/ZFpqNcw3FFRCJic1fbzPDb6VFWbc3xs0ZZbb9hU+hDdEOtcZjZ583smJn1FXjfzOxvzGy/mT1jZlfXOkYRkTBEeZRV2E1Vfw9sm+P9twGXBT/bgU/XICYRkdBFeZRVqInD3R8FTs5R5Dbgi571GNBuZl21iU5EJDxRHmUVdo1jPt3AoZznh4PXzmFm282s18x6h4aGahKciEi1RHmUVdQTh+V5zfMVdPcd7t7j7j0dHR1VDktEpLryjbKKQsc4RH9U1WFgQ87z9cDRkGIREamp3FFWURL1GsdO4JeC0VU/BiTcfSDsoERElrJQaxxm9hXgRmCNmR0G/giIA7j7Z4BdwK3AfuA08MvhRCoiUr7Fctn0+YSaONz95+d534HfqFE4IiJVE+UJfaWKeh+HiMh5YXffIFNTU/QPjDCSTLGiKc7aFY3s7htU4hARkXPtHUhw8MRpmuL1tDbGSKYyPPfKKKdTmbBDK1nUO8dFRM4LifE0ZkZTvP6s34nx9Px/HDFKHCIiNbCiKQYOyVQGdyeZyoAHry8yiy9iEZFF6Ip1bSyL1/PK6ARjyTQtTTE2rlrGpo6WsEMrmWocIiI1sG1rJ/X1dWzpWsFNmy9gS9cK6uvrInEJkVIpcYiI1ECULyFSKjVViYjUSFQvIVIq1ThERKQkShwiIlISJQ4RESmJEoeIiJREneMiIhGxWK6eqxqHiEgETF89NzGeOuvquf0DibBDO4dqHCIiRah2bWB33yCZzBR7B0ZmZpavbY3m1XNDrXGY2TYze97M9pvZPXnef7+ZDZnZ08HPB8KIU0SWtlrUBp49mmDf4BgTqQwtjfVMpDLsGxzj2aOqccwws3rgU8AtZO8t/oSZ7XT3vbOK3u/ud9U8QBGRwO6+Qdqa47Q1xwFmfleyNjCSTINBU7weyP6eSE9lX4+YMGsc1wL73f2Au08C9wG3hRiPiEheR4bHaZ11FdvWphhHhscrto625tjMVXNzf7c1R69HIcyIuoFDOc8PA9flKfczZnYDsA+4290P5SkjIlI13e3NvHR8jFdGJs66e99Fa1oq1vexpSu4em7OOi5avYyL1kTv6rlh1jgsz2s+6/k/Axe5+5XAw8AXCi7MbLuZ9ZpZ79DQUAXDFJGl7vLO5Tx1cJjEeIqWhnoS4ymeOjjMsgarWN/Htq2d1NXVsblrBTdv7mRz1wrq6qJ59dwwE8dhYEPO8/XA0dwC7n7C3SeCp38HvKHQwtx9h7v3uHtPR0dHxYMVkaVr3+AprtrQzormOKcmp1jRHOeqDe088tzxmb6POrOZx7v7Bktex2K6em6YTVVPAJeZ2SbgCHAH8Au5Bcysy90HgqfvBPprG6Is1GKZyCThWiz7yZHhcTauWX7WTZem3NlzJMF1m1bNvHZ8LMmPBscYHM2e75b6eRbL1XNDq3G4exq4C3iQbEL4qrs/a2Z/YmbvDIr9lpk9a2Z7gN8C3h9OtFKKxTSRScKzmPaT7vZmRmeNbhpNpulc0TTz+vGxJE++PMxIMk1na2OkP0+5Qp3H4e673P1yd7/E3f9X8Nq97r4zePwH7n6Fu7/O3X/C3Z8LM14pTu7QxXKr73L+Wkz7ybatnSTGUyTGU0y5zzx+3/UXzjz+0eAYE+kMQ6MTHBtN0j8wwtTUVCQ/T7miN85LFr0jw+N0tTWd9Vqlhy4uRFSbRQrFFdV4KyWq+0k+0/0Puf+P269Zz+auNi7uaGF33yCHXx1nIp1h1fIGVi5rIJnK8Nwro5xOZcIOv+KUOKTiutubSYynZiZJQbZa393eHFpM080ibc3xs5pFwu58zBfXR3Y/T3O8jr6jo6xcFmdzV2tk4q2kKO4ncynU/zD9+vcPnGAk5/NMT+BLjEdvAl+5lDik4rZt7WTHoy8C2TPI0WSaxHiK269ZD4Rz5l+Lmb+ViGsyneHgidOMTqRYuyJ7Nv70oQRXX9g+04xzviSO+faTMJSzb65oijFyOkUylaExVsdEego8+/r5RlfHlYqba1hhJTpE+wcSfPyhfXzwa3v4+EP7ivrbWsz8XYjZce0fOpW9TlF6iqZ4PU3xehpjdewfOhWJeCspasNPy903r1jXxuWdLTTG6xmbyNAYr+fyzhauWHd+JPpc518qlEgoVK0v98x/oU1OUW0WmR3XWDJNrA6WN8ZmkkdjrI6xZHpB8Uapn6RQLFGpQc23b863LbM1qNNs6VpxVg0qihP4yqUahwALO4tfiHLP/Bc6EqfQqJhSv9SV3k6z44rXG2MTGV6/vo2J9BTJVIaJVIZ4vZUcb5SGu0YplkLm2jeLiT9qNahqUo1DqtJxXOjsrNwz/4WOxJlrVEwpn6nS2yk3rr0DCdJTzunJNEeGk2xa3czRxASvjqd586Wr+cXrN5a0nij160QplkLm2jeLjT9KNahqUuKQor8U/QMJvvT9l/nBoWEM46oNbbw3z8FsrgNsuR2i5SSecr/U1Tr4Tf/twZOnueaiZSRTafoHRnl2YIwfv3R13m1cjCgNd41SLIXMtW9+7nsvRT7+WlJTlRTVfNQ/kOAju5/n8QMnaagzYnXw/QMn+diD5zbXzNWcVG51vlJNTgtRzQ723G3WuaKZG19zAW957QWsaW1acFIqNNs5jH6dKMVSyFz75mKIv5ZU45CizuJ39w1y8tQkLU2xmRvNmBnHxybOOePOPbscGk2yf+gUo+MpMMruEK1Ek9NCVbODvRpn5NUa7rqQDvdyY6lVJ3+hfTOKQ4fDpBqHFHUWf2R4nMn0FI2xM7tMdqx65pyD2/TZ2dBokqcODs907jbU11WkQ3RzVxt333I5H3vP67j7lstr1qZczdpONc5oq9FZu9BO7ulYJtMZHu4f5D9ePMmyeHGHnyh0rC+lju9iqMYhRZ3Fd7c386PB0ZkhogAT6SkaY/XnHNymz84ODI3RWJ+97cpkxnnDxjbi9fVV6xCt9lnp7I7sxHiaFU2xmRFd5ayrWme0le6sLbefZzw1xXWbVs98xmIGFxSzzlrUSJZKx3cxVOMQYP6z+G1bO1m1vIGxZJrkZJrxyey8gjUtjeeccU8fYFMZZyKTTTRv2NjOmpamqnUo1uqsdHNXG9u2dtLSGKdrRSMj4yn+5Zmj3H3/Hr75zJGylrsYzmjL6edZ6FDq+dYZhRrJUqMahxRlc1cbv7ftNWeNqrr+4lUFR/xs7mrjli2dNZt0V8vhnrv7BpmamuL5wTEaY3WsXt7ASDLNJx55gYs7Wsrqv4laopitnH6ehfbjzLfOxTDU93yjxCFF29zVxp+++8qiy+drfnn5xCnWtTXxwa/tqWiTQi2Hex4ZHmcgkaQxVjfTbLeiKcbJU6nz/mBVTpPaQpPOfOtcDEN9zzdqqpKqmd38MpnOUGdGPFZf8SaFWg6X7G5v5sSpybMGCkykp1i1PF7ywapWM/Yrtf5ymtQWOrhgvo51DZWtPXP38FZutg34a6Ae+Ky7f3jW+43AF8nea/wEcLu7vzTfcnt6ery3t7fyAUdMrYYoVmo9H39o31lnnMfHkvQdGWEyM8Vbt6wtK/7cSYe5Z6W5B7VKfY7+gQR3378HyNY0JtJTJE5P0hCrJ1Zfxy1bOotadjExV1M561/otizl73LLNtQbgyMTbFi17JxYgVC34/nCzJ50956iyoaVOMysHtgH3AIcJnsP8p939705ZX4duNLd7zSzO4B3ufvt8y17KSSOWh10KrmeD35tD11tTdSZzdxms6HeSGWc6y5eXXb8cx2U5vscpR4Iv/nMET7xyAukp5zGmJEYT9MQq+O6i1fSGIsV9VlmJ1Jg5vndt1w+72cqVzHrz6ea+9705332aILDr47zms4WLly9nEf3DXF8dIJVLQ1kpqClKcba1kY2dbRw9y2XR+pijotVKYkjzD6Oa4H97n4AwMzuA24D9uaUuQ344+Dx14FPmpl5mNWkiKhVh+BC15Pvi5zbxr3/2KmZpp7W5lhF4p+rc3muzwEUvERKoeW9/crumTu/fXvvK6xuaWBr9wrWtJxpa5/vs8zXNl/tm08ttG+gWvte7ucdGU8B8PzgGC1NMUbGU5yazJAeSbJx9XImUhn2DY7N3F1vMQwsOJ+EmTi6gUM5zw8D1xUq4+5pM0sAq4HjsxdmZtuB7QAXXnhhNeKNlFp1CM5ez9Bokv3HxhgcnQDIe2ZX6IB38+YOHu4fAiAxPkljfR0TGeeKdStKir/Ys8vccnuPjvC6DSuAM2fX0+tbyIEwd9mGsWVd61lJYzGMFlpoZ3W19r3pzzuZznDg+CkA6s3YcyjBRGaKWB1kPHvFgum7640k899dTzWQ6gqzc9zyvDa7JlFMmeyL7jvcvcfdezo6OsoOLgrm6risVYdg7nqmZ4KPJNN0tjYW7NwuNF5/3+CpmY7VOqsDMzatXsb+oVM8tHeQR/cN0VCf719+RrFj9meXi9cbjx94leNjyZky09ur1LkJpSw7399O/0+HRpMcOnm6YGdxOXMmiun0XmhndbX2vSPD4yRTaZ46OEydZQ9OjnPw1dMYkM449Wa4O8lUBnenrfncc1/N66i+MBPHYWBDzvP1wNFCZcwsBrQBJ2sSXcjm2/lrdbG/3PXsPzY28/plnS0FJ3DNdcCbnmj4Fz93JauWN7BvcIyJyezNi0aTaQZHJub8ghc7iWx2uSvWrcCAviMj52yvUg+EpSw71+z/aUOsnil3UulM3hFKCz1AF3vgXOgIqWrte93tzfQPjNIYq+OC1iYyDpkpaIrVMT6ZITPlTAS1kakp57VrW9kyz83CSploWIywR8FFRZiJ4wngMjPbZGYNwB3AzllldgLvCx7/LPCdpdK/Md/OP/tLn0pnaI7X8bnvvVTRHTp3PYOjE6xois3MAof8Z8DFHPA2d7Wxrq2JlqYYk1NOc0OM6y9ZxYZVy+b8ghd7Fr53IMHeowm+vfcVHjtwAjO4ZtNKJjNT5xwkSz0Qzo6ho7Wp4LJz5fufbly9nDWtTXln7C/0AF3KgXMh1/2q1iz3bVs7efV0CndnWUM9q5c3ABCrMybSGdqWNbBx1TLWtTVxajLD2EQm77ao1lWMVZM5I7Q+jqDP4i7gQbLDcT/v7s+a2Z8Ave6+E/gc8CUz20+2pnFHWPHWWjHtyNMdgrl9CqtbYhXvRJ1exvcPnODE2AR2LNv+vKaliYMnTjEwMnHWhL5iJ4lNZJwbLu+gzs40T025z/kFb6w3Ht03xGRmihVNcS69YDnx+rOvl9U/kODQiXGw7HDZZCrDky8P85rOFt66Ze05I4ZKveJuvr6Bpngs77Jzldo3sNArAdei/6sandGbu9r48UtX03d0hNGJNCuXN3DNppX0HRmho7WJK9atYP/QKcaSaVqbYnSuaMwbQ7WuYqwZ6meEOnPc3XcBu2a9dm/O4yTwnlrHFQWl7PzV3qGnE9Pa1kZGTmfPentfepWutiZeGDrFVRvazxn1U8wBr9Qv+DefOcL/23+c46cmWdZQTyYzxfdfmGDTmuXc/pNnDti7+wa5vLOFfcfGggsx1jGRzl4i5M4bL8m77FIOhAudPb2QA9pCDtBRvb96Md57/cZzhvq+ejrFGy9ZRUdrEx2t2YQ45c5AIpl3GdW6YKRmqJ8xb1OVmd1lZitrEYycUUozRTVvMARnEtOmjhau3tjOiuY4mSnYd2yMqza0s6mjJe8Nm+ZrAinlM/YPJPjEIy/Q1BBj46plGMbg6ATxejvnzPPI8Dgb1yzn6gvbaYzXMzaRYUVTjPUrm8+a17HQtuqo9Q2EtZ5qyLdt33zpapriZ+/fcyXCajWlaYb6GcXUONYCT5jZU8DngQeXSj9DmEpppqj2GWbumdb0Wd+UOzv3HGXjmuVnlS0lYZVymfLdfYOkp5xVy+OYGS1NcZKpDI2xOiYzZ++O09sj9ww1d/tUYn7EQmoCC216KlWt1lMts7ft9P8Liq9BVKMpTTdzOmPexOHu/8PM/hB4K/DLZCfhfRX4nLu/UO0Al7Jid/5q79CFElPniiZGk+myEtb05zt48jTd7dnLSeQ7kB8ZHmfV8vhZ9wNpjNVx4tQk11+y5qxlzrc9SmnaW8ickbnKVeqANt/6ajkhrpb3QQkzEUYljigoqo/D3d3MXgFeAdLASuDrZvaQu/9eNQOUM2Zfu8fIdjB3tzdz8+YO9g2eqsoOXehA/L7rL5yZ0FdOwirmQN7d3sxkKsO+YEhwY6yOkWSaeH1dwfuBFPqCF9tWXWzNpNozvGer9fqiEEtUZoZHJY6wzZs4zOy3yA6JPQ58Fvhdd0+ZWR3wI0CJowZyv6CxOnj8wEkcuO7ilSTGUzzcP1S1A8dcB+Lpy26Uk7CKOZBnk9dpLr+ghVdGkpw8lSJWZ/zmTZeUfGZfbNNesTWTWo+2idLonkrHohnfi0MxNY41wLvd/eXcF919yszeUZ2wZLYvff9lDgyNkco4w+OTLG+oZ1lDjANDp/mxi1cD1T1wFDoQV+IMrJgDeW7yaojXc/0lCz+oFNu0V2zNpNajbaI0uqeSsUyfHGUyU7wykuTpg8M82PcKv3nTJbz9yu5KhSwVUEwfx71zvNdf2XAkn/6BBN/bf4KVzTFaGus5dDLFq6cmaYjVEa+v45KO5axuaeTI8PiiPGMr9kBeqWaCYtuqi62Z1Hr4a5SG21Yylt19g2QyU+w7lr2z4qrlcYZGJ7j3G3t54AdHuWJd26LYn5cC3chpEdjdN8jKZXFOpzK8ePwUpyczTKSnSE6mqTd46uAwLx8/RUO9LcqZrdUaPjnfOis1XLjWw1+jNNy2krEcGR7nlZEzd1Y8PZlhJJkmlZliJFjuYtiflwLdOnYRyDYHNPLvL5wilZkiVpe9h8VkBloasyOM9g2O8YaN7ZFp+y5VlDodc2ttzfG64FpS6YI1k1qPtonS6J5KxtLd3szTB4dZtTy73548PQlk9/Gxicyi2p/Pd0oci0B3ezPffX6UrrYmDp48TV2d0WAGOIlkhovXNNG2LM5ExlnVUr2JgEvB7FFCxd6kqNaJL0qJtlKxbNvayYN9rzCSzM7nGZ/MUGfQ2tRASzDBVftzNKipahGYvvhbY6yOlcsaWNYQY3ljjEs6WrigtYkt69rY0tWmma0VUM0rq8rcNne18Zs3ZS8Jc+LUJE3xOlY0Zyd8XtqRnWiq/TkalDgWgemLv2VvYJP9l61e3kCszojX20ybcpTavheral++Reb29iu7+fjtr+MdV65ja3cbjfX1XH5BC6tbGrU/R4iaqhaJ3Iu/TaTT7D06yqunU7z50tX84vUbZ5oKotL2vVhFacTSUpXb9DV7lKD252iw8/GyUz09Pd7b2xt2GBW3GIfaLja5fRy5Q4PDmJUtUktm9qS79xRTVjWORSRKHaLns2XxOv7jxZM4zlUb2pU0RGYJJXGY2SrgfuAi4CXg59z91TzlMsAPg6cH3f2dtYpRlp7c2sZNmy+YqW2IyNnC6hy/B3jE3S8DHgme5zPu7q8PfpQ0pKo0okqkOGE1Vd0G3Bg8/gLwXeD3Q4qlJtQ/EX1RugaUSJSFVePodPcBgOD3BQXKNZlZr5k9ZmY/PdcCzWx7ULZ3aGio0vGWRTe5Xxw0D0akOFVLHGb2sJn15fm5rYTFXBj08v8C8Fdmlv+G0YC773D3Hnfv6ejoKDv+SlITyOKgeTAixalaU5W731zoPTMbNLMudx8wsy7gWIFlHA1+HzCz7wJXAYvuroNqAlkconQNKJEoC6uPYyfZm0N9OPj9jdkFzGwlcNrdJ8xsDfAm4CM1jbJCNKls8dCQZ5H5hdXH8WHgFjP7EXBL8Bwz6zGzzwZlNgO9ZrYH+Ffgw+6+N5Roy6QmEBE5n2jmeI1oVJWIRJlmjkeQmkBE5Hyhq+OKiEhJlDhERKQkaqqqMvVtiMj5RjWOKtKMcRE5HylxVJFmjIvI+UiJo4p0G1IROR+pj6NMc/VhaMa4iJyPVOMow3x9GJoxLiLnIyWOMszXhzF90by25jgDiSRtzXHdhlREFj01VZWhmKveasa4iJxvVOMog278IyJLkRJHGdSHISJLkRJHGdSHISJLkfo4yqQ+DBFZalTjEBGRkoSSOMzsPWb2rJlNmVnBG4eY2TYze97M9pvZPbWMUUR/XizJAAALMElEQVRE8gurqaoPeDfwvwsVMLN64FNkby17GHjCzHZG6faxuvKtiCxFodQ43L3f3Z+fp9i1wH53P+Duk8B9wG3Vj644uvKtiCxVUe7j6AYO5Tw/HLyWl5ltN7NeM+sdGhqqenC68q2ILFVVSxxm9rCZ9eX5KbbWYHle80KF3X2Hu/e4e09HR8fCgi6BrnwrIktV1fo43P3mMhdxGNiQ83w9cLTMZVaMrnwrIktVlJuqngAuM7NNZtYA3AHsDDmmGZo1LiJLVVjDcd9lZoeB64FvmtmDwevrzGwXgLungbuAB4F+4Kvu/mwY8eajWeMislSFMhzX3R8AHsjz+lHg1pznu4BdNQxNRETmEeWmqkjTcFwRWaqUOBZIw3FFZKlS4lggDccVkaVKiWOBdBMnEVmqlDgWSMNxRWSpUuJYIA3HFZGlSjdyKoNu4iQiS5ESRwG6ZLqISH5qqspDczRERApT4shDczRERApT4shDczRERApT4shDczRERApT4shDczRERApT4shDczRERArTcNwCNEdDRCS/UBKHmb0H+GNgM3Ctu/cWKPcSMApkgLS799QqRs3jEBHJL6ymqj7g3cCjRZT9CXd/fa2ThuZxiIjkF0ricPd+d38+jHUXQ/M4REQKi3rnuAPfNrMnzWz7XAXNbLuZ9ZpZ79DQUFkr1TwOEZHCqtbHYWYPA2vzvPUhd/9GkYt5k7sfNbMLgIfM7Dl3z9u85e47gB0APT09vqCgA93tzSTGU7Q1x2de0zwOEZGsqiUOd7+5Ass4Gvw+ZmYPANdSXL9IWbZt7WTHoy8C2ZrGaDJNYjzF7desr/aqRUQiL7JNVWa23Mxapx8DbyXbqV51mschIlJYWMNx3wV8AugAvmlmT7v7T5rZOuCz7n4r0Ak8YGbTcf4fd99dqxg1j0NEJL9QEoe7PwA8kOf1o8CtweMDwOtqHJqIiMwjsk1VIiISTUocIiJSEiUOEREpiRKHiIiURIlDRERKosQhIiIlUeIQEZGSKHGIiEhJlDhERKQkShwiIlISJQ4RESmJEoeIiJREiUNEREqixCEiIiVR4hARkZIocYiISElCSRxm9lEze87MnjGzB8ysvUC5bWb2vJntN7N7ah2niIicK6wax0PAVne/EtgH/MHsAmZWD3wKeBuwBfh5M9tS0yhFROQcoSQOd/+2u6eDp48B6/MUuxbY7+4H3H0SuA+4rVYxiohIflHo4/gV4Ft5Xu8GDuU8Pxy8lpeZbTezXjPrHRoaqnCIIiIyLVatBZvZw8DaPG99yN2/EZT5EJAGvpxvEXle80Lrc/cdwA6Anp6eguVERKQ8VUsc7n7zXO+b2fuAdwA3uXu+A/1hYEPO8/XA0cpFKCIiCxHWqKptwO8D73T30wWKPQFcZmabzKwBuAPYWasYRUQkv7D6OD4JtAIPmdnTZvYZADNbZ2a7AILO87uAB4F+4Kvu/mxI8YqISKBqTVVzcfdLC7x+FLg15/kuYFctYuofSLC7b5Ajw+N0tzezbWsnm7vaarFqEZFFJQqjqkLXP5Bgx6MvkhhP0dXWRGI8xY5HX6R/IBF2aCIikaPEAezuG6StOU5bc5w6s5nHu/sGww5NRCRylDiAI8PjtDad3WrX2hTjyPB4SBGJiESXEgfQ3d7MaDJ91mujyTTd7c0hRSQiEl1KHMC2rZ0kxlMkxlNMuc883ra1M+zQREQiR4kD2NzVxvYbNtHWHGcgkaStOc72GzZpVJWISB6hDMeNos1dbUoUIiJFUI1DRERKosQhIiIlUeIQEZGSKHGIiEhJlDhERKQklv9WGIubmQ0BLy/gT9cAxyscTiVENS6IbmxRjQuiG5viKl1UY1tIXBvdvaOYgudl4lgoM+t1956w45gtqnFBdGOLalwQ3dgUV+miGlu141JTlYiIlESJQ0RESqLEcbYdYQdQQFTjgujGFtW4ILqxKa7SRTW2qsalPg4RESmJahwiIlISJQ4RESnJkk4cZvZRM3vOzJ4xswfMrL1AuW1m9ryZ7Teze2oQ13vM7FkzmzKzgkPqzOwlM/uhmT1tZr3VjqvE2Gq9zVaZ2UNm9qPg98oC5TLB9nrazHZWMZ45P7+ZNZrZ/cH7j5vZRdWKZQGxvd/MhnK20wdqFNfnzeyYmfUVeN/M7G+CuJ8xs6sjEteNZpbI2V731iiuDWb2r2bWH3wn/1ueMtXZZu6+ZH+AtwKx4PGfA3+ep0w98AJwMdAA7AG2VDmuzcBrgO8CPXOUewlYU+NtNm9sIW2zjwD3BI/vyfe/DN4bq8E2mvfzA78OfCZ4fAdwf43+f8XE9n7gk7Xcr4L13gBcDfQVeP9W4FuAAT8GPB6RuG4E/iWE7dUFXB08bgX25flfVmWbLekah7t/292n7xn7GLA+T7Frgf3ufsDdJ4H7gNuqHFe/uz9fzXUsVJGx1XybBcv/QvD4C8BPV3l9cynm8+fG+3XgJjOziMQWCnd/FDg5R5HbgC961mNAu5l1RSCuULj7gLs/FTweBfqB7lnFqrLNlnTimOVXyGbm2bqBQznPD3PuPycsDnzbzJ40s+1hB5MjjG3W6e4DkP1CARcUKNdkZr1m9piZVSu5FPP5Z8oEJy8JYHWV4ik1NoCfCZo2vm5mG2oQVzGi/F283sz2mNm3zOyKWq88aOq8Cnh81ltV2Wbn/R0AzexhYG2etz7k7t8IynwISANfzreIPK+VPYa5mLiK8CZ3P2pmFwAPmdlzwdlR2LHVfJuVsJgLg212MfAdM/uhu79QbmyzFPP5q7KNilDMev8Z+Iq7T5jZnWRrRm+pemTzC2ubzecpstd5GjOzW4F/Ai6r1crNrAX4B+C33X1k9tt5/qTsbXbeJw53v3mu983sfcA7gJs8aBSc5TCQe8a1Hjha7biKXMbR4PcxM3uAbDNE2YmjArHVfJuZ2aCZdbn7QFAVP1ZgGdPb7ICZfZfsWVqlE0cxn3+6zGEziwFt1KY5ZN7Y3P1EztO/I9v/FwVV2a/KlXuwdvddZva3ZrbG3at+8UMzi5NNGl9293/MU6Qq22xJN1WZ2Tbg94F3uvvpAsWeAC4zs01m1kC2I7Nqo3GKZWbLzax1+jHZjv68oz5CEMY22wm8L3j8PuCcmpGZrTSzxuDxGuBNwN4qxFLM58+N92eB7xQ4cal5bLPawN9Jtu08CnYCvxSMFPoxIDHdPBkmM1s73T9lZteSPa6emPuvKrJeAz4H9Lv7XxYoVp1tVuuRAFH6AfaTbf97OviZHuWyDtiVU+5WsiMWXiDbXFPtuN5F9kxhAhgEHpwdF9lRMXuCn2drEVexsYW0zVYDjwA/Cn6vCl7vAT4bPH4j8MNgm/0Q+NUqxnPO5wf+hOxJCkAT8LVgH/wP4OJa/P+KjO3Pgn1qD/CvwGtrFNdXgAEgFexjvwrcCdwZvG/Ap4K4f8gcIw5rHNddOdvrMeCNNYrrzWSbnZ7JOYbdWottpkuOiIhISZZ0U5WIiJROiUNEREqixCEiIiVR4hARkZIocYiISEmUOESqLLiK6Ytmtip4vjJ4vjHs2EQWQolDpMrc/RDwaeDDwUsfBna4+8vhRSWycJrHIVIDwaUhngQ+D/wacJVnr04rsuic99eqEokCd0+Z2e8Cu4G3KmnIYqamKpHaeRvZS1dsDTsQkXIocYjUgJm9HriF7F3Y7q7FDYhEqkWJQ6TKgquYfprs/RIOAh8FPhZuVCILp8QhUn2/Bhx094eC538LvNbM/lOIMYksmEZViYhISVTjEBGRkihxiIhISZQ4RESkJEocIiJSEiUOEREpiRKHiIiURIlDRERK8v8BW4Jh59b88h8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(df_test.X,df_test.Y,'o',alpha = 0.5)\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('y')\n",
    "plt.title('Test Data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use Polynomial Regression with $p \\in \\{1,...,10\\}$ such as:\n",
    "\n",
    "$$\n",
    "Y = \\beta_0 + \\beta_1X + \\beta_2X^2 + ... + \\beta_pX^p + \\epsilon \n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathcal{D} := \\{(x_i,y_i)\\}_{i=1}^n \\in \\mathbb{R}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\begin{bmatrix}\n",
    "    y_1 \\\\\n",
    "    y_1 \\\\\n",
    "    \\vdots \\\\\n",
    "    y_n \\\\\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "    x_{1}       & x_1^2 & \\dots & x_1^p \\\\\n",
    "    x_{2}       & x_2^2 & \\dots & x_2^p \\\\\n",
    "    \\vdots      &\\vdots & \\ddots & \\vdots \\\\\n",
    "    x_{n}       & x_n^2 & \\dots & x_n^p\n",
    "\\end{bmatrix} \n",
    "\\begin{bmatrix}\n",
    "    \\beta_0 \\\\\n",
    "    \\beta_1 \\\\\n",
    "    \\vdots \\\\\n",
    "    \\beta_p \\\\\n",
    "\\end{bmatrix} + \n",
    "\\begin{bmatrix}\n",
    "    \\epsilon_1 \\\\\n",
    "    \\epsilon_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    \\epsilon_n \\\\ \n",
    "\\end{bmatrix}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Least Square solution of this equation is:\n",
    "$$\n",
    "\\hat{\\beta} = (X^TX)^{-1}X^Ty\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will write a function for solving least-square estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [],
   "source": [
    "def least_square_estimate(X,y):\n",
    "    '''\n",
    "    X:Design matrix\n",
    "    y:True response\n",
    "    '''\n",
    "    X = X if isinstance(X,np.ndarray) else np.asarray(X)\n",
    "    y = y if isinstance(y,np.ndarray) else np.asarray(y)\n",
    "    b_hat = inv(X.T@X)@X.T@y\n",
    "    return b_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vandermond matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [],
   "source": [
    "def X_design(X, p=1):\n",
    "    return np.vander(X, p+1, increasing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(X,Beta,y_hat = None):\n",
    "    if y_hat: # mean of Dtrain\n",
    "        return X@Beta + y_hat\n",
    "    return X@Beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MSE(y_hat,y):\n",
    "    return (y_hat - y).T@(y_hat - y)/(y.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will calculate the $\\hat{y}$ for $p \\in{1,...,10}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_MSE(X_tr,y_tr,X_te,y_te,p):\n",
    "    training_errors = []\n",
    "    test_errors = []\n",
    "    \n",
    "    for i in range(p):\n",
    "        X_des = X_design(X_tr,i) #vandermond\n",
    "        beta_hat = least_square_estimate(X_des,y_tr) #beta hat \n",
    "        training_errors.append(MSE(predict(X_des,beta_hat),y_tr)) #train error\n",
    "        test_errors.append(MSE(predict(X_design(X_te,i),beta_hat),y_te))  #test error\n",
    "        \n",
    "    plt.plot(training_errors, marker='.', color='blue', linewidth=2, label='train error')\n",
    "    plt.plot(test_errors, marker='.', color='red', linewidth=2, label='test error')\n",
    "    # Add legend\n",
    "    plt.legend()\n",
    "\n",
    "    # Add titles\n",
    "    plt.title(\"Train&Test Errors\")\n",
    "    plt.ylabel(\"MSE\")\n",
    "    plt.xlabel(r'$p$')\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEYCAYAAACtEtpmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VfW19/HPIsyTjCqCCCoqgwwS0YgiODCIgtZLFcWpWqyPQ22vXvW2VqpPW1vuvQ+1tbVeBayiaGkVrFSmiiiiEhQHQAQRJSIlTAIyJlnPH3uHHEKSk+HsbHLyfb9e55U973UOetbZv/37rW3ujoiISFnqxB2AiIgc/pQsREQkKSULERFJSslCRESSUrIQEZGklCxERCQpJQtJe2aWYWY7zaxj3LGI1FRKFnLYCb/YC18FZrY7Yf7qih7P3fPdvam7f1nO85uZTTCzzeFrarH1sxPi2W9m+xLmf1/R+BKO+19m9liSbTaZ2a5in9GvK3tOkfKqG3cAIsW5e9PCaTNbC9zk7nNL297M6rp7XgpDuAi4AugJbAbOLhbf4IRzPwOsdvdxKTx/Mue5+9vJNirpc6noZxXBZys1lK4spMYxs/9rZs+b2XNmtgMYY2ZZZva2mW0zs6/N7BEzqxduX9fM3Mw6hfPPhOv/YWY7zGyRmXVOOEUesAvY4O57ykpUZcR4uZl9FMbzupmdkrBuXBjjdjNbYWZnmdnlwB3AjeHVwluVOOdtZjbHzB4zs63AXaUsywg/w3VmtsHMnjCzpuExepjZHjP7gZnlADPMrKmZvWBmW8xsa/g5H1HR+KRmU7KQmuoy4FngCOB5gi/4HwJtgP7AUODmMva/CrgfaAV8CTyUsG4ZcBTwJzOzigZmZmcDvwOuA1qH8b0Yfkn3Ba4muGo5ArgY+Mrd/wo8AjwZNpmdVdHzhgYB7xJ8Dr8tZdmtBJ9ff+BkoD0wPuEY9YE+QBfgcoLP0YFjgLYESW1fJeOTGkrJQmqqN939ZXcvcPfd7r7Y3d9x9zx3XwM8Dpxbxv7T3D3b3fcDU4DeAGZWH5gFfB9oBzxWmDDM7B0zG1aO2G4Gfuvu74X3S/4ANA/PkQc0BroDGe7+mbt/UcH3Pie8Yil8XZWw7lN3nxied3cpy64Gfu3uX7r7N8BPgTEJxzDg/vBz3Q3sJ0gSx4ef77sJx5ZaQslCaqp1iTNmdoqZvRI2q2wHHiT4JV2aDQnTu4DC+yQXAo3c/Tng34CuBAmjBXA8sLAcsR0H/CzxCx1oCbR39w8Ivpx/BWw0s6fNrG05jpnoQndvkfB6NmHduhK2L77sGCAxQX0BNE1oWtrn7hsT1j8OvEVwdbQubMLSd0cto39wqamKl0v+E/AxcKK7Nwd+RvALuaLqEvz6J/z1fAlwOvAO8IS7by/HMdYB/1nsC72xu88IjzvJ3bOAE4BmBImtpPdUGSUdo/iy9QQJrVBHYGd4lXHI9uF9m5+6+8kETVpXAaNSEKvUIEoWki6aAd8A35pZV8q+X1GWBUBzM3vAzBqFy+YDJwEF5TzGn4A7zey0sBtuMzO71MwahjeQB5hZA4Irmj1Afrjfv4DOlblPUkHPAXebWQcza05wv2ZKaRub2YVm1jW8mthOkEzzS9te0pOShaSLfye4obyD4Mv6+cocxN23AoOBcwh+gX8A1AP6Ajeb2Q3lOMYC4EfAE8A2YCVBV1wHGgH/j6BL7tdAA2BcuOuzQAtgi5m9UcYp/llsnMXTFXybvwdeJrhaWkWQpO4uY/tjw+13EHwe04G/VfCcUsOZHn4kIiLJ6MpCRESSUrIQEZGklCxERCQpJQsREUkqbQoJtmnTxjt16hR3GCIiNcqSJUs2uXvSgaFpkyw6depEdnZ23GGIiNQoZlaucjNqhhIRkaSULEREJCklCxERSSpt7lmISHrYv38/OTk57NmzJ+5Q0krDhg3p0KED9erVq9T+ShYicljJycmhWbNmdOrUiehrKtYO7s7mzZvJycmhc+fOyXcogZqhROSwsmfPHlq3bq1EkUJmRuvWrat0taZkAbBoEfzqV8FfEYmdEkXqVfUzVTPUhAnw7/8O7tCwIcybB1lZcUclInJY0ZXFe+9BQUGQLPbtg/nz445IRGK0bds2/vCHP1Rq34suuoht27alOKLDg5LFddcVTderBwMHxhaKiMSvrGSRn1/2AwJnzpxJixYtUhpPXl5emfOlSRZrRSlZnH8+ZGYG0/fcoyYokRoolbcd7733Xj777DN69+7N3Xffzfz58xk0aBBXXXUVp556KgCXXnopffv2pXv37jz++OMH9u3UqRObNm1i7dq1dO3ale9///t0796dwYMHs3v37kPOlZuby+WXX87pp5/O6aefzsKFCwEYN24cY8eOZfDgwVx77bVMnjyZUaNGcckllzB48GDcnbvvvpsePXpw6qmn8vzzwYMhS4o1Zdw9LV59+/b1Svvv/3YH96uvrvwxRCQlli9ffmA6aB9O/assn3/+uXfv3v3A/GuvveaNGzf2NWvWHFi2efNmd3fftWuXd+/e3Tdt2uTu7scdd5zn5ub6559/7hkZGf7++++7u/uoUaP86aefPuRco0eP9jfeeMPd3b/44gs/5ZRT3N39gQce8NNOO8137drl7u6TJk3y9u3bHzjvtGnT/IILLvC8vDzfsGGDH3vssb5+/foSYy3tsy36jMn2cnzH6gY3wEUXBTe5X30V8vMhIyPuiETkMNKvX7+Dxic88sgjvPjiiwCsW7eOVatW0bp164P26dy5M7179wagb9++rF279pDjzp07l+XLlx+Y3759Ozt27ABgxIgRNGrU6MC6Cy+8kFatWgHw5ptvMnr0aDIyMjjqqKM499xzWbx4Mc2bNz8k1lRRsgA4+WTo1AnWroXsbDjjjLgjEhGC64BkFi0KWpP37YP69aPp0NikSZMD0/Pnz2fu3LksWrSIxo0bM3DgwBLHLzRo0ODAdEZGRonNUAUFBSxatOigpFDSOYvPexkfTPH9UkX3LADMgqsLgJkz441FRCokKytIEA89lJpE0axZswO/7kvyzTff0LJlSxo3bswnn3zC22+/XelzDR48mN///vcH5pcuXVqu/QYMGMDzzz9Pfn4+ubm5LFiwgH79+lU6jvJQsig0bFjw9x//iDcOEamwrCy4777UXFG0bt2a/v3706NHD+6+++5D1g8dOpS8vDx69uzJ/fffz5lnnlnpcz3yyCNkZ2fTs2dPunXrxmOPPVau/S677DJ69uxJr169OO+88/jNb37D0UcfXek4ysPKupypSTIzM71KDz/69lto1Qr274cNG+DII1MXnIiU24oVK+jatWvcYaSlkj5bM1vi7pnJ9tWVRaEmTYIxFu4wa1bc0YiIHFYiTRZmNtTMVprZajO7t4T1PzCzj8xsqZm9aWbdEtbdF+630syGRBnnAWqKEhEpUWTJwswygEeBYUA3YHRiMgg96+6nuntv4DfA/4T7dgOuBLoDQ4E/hMeLVuFN7lmzgi60IiICRHtl0Q9Y7e5r3H0fMBUYmbiBu29PmG0CFN5AGQlMdfe97v45sDo8XrS6dIHjj4ctW+DddyM/nYhITRFlsmgPrEuYzwmXHcTMbjWzzwiuLO6o4L5jzSzbzLJzc3OrHrGZmqJEREoQZbIoqXj6IV2v3P1Rdz8BuAf4aQX3fdzdM909s23btlUK9gCNtxAROUSUySIHODZhvgOwvoztpwKXVnLfSlu1Cm66CcaMCYuQDRwIDRrAkiXwr39FcUoROYxVpUQ5wIQJE9i1a1cKIzo8RJksFgNdzKyzmdUnuGE9I3EDM+uSMDscWBVOzwCuNLMGZtYZ6AJEchNh0iR48kmYMiUoGbDog8YwaFCwUl1oRWqduJNFZUuSl3e7yoqsNpS755nZbcAsIAOY6O7LzOxBgiqHM4DbzOwCYD+wFbgu3HeZmb0ALAfygFvdPZLuSYlPGix89lHWsGFBUcGZM+Haa6M4rYik0qJFwf+8AwdWeRh3YonyCy+8kPHjxzN+/HheeOEF9u7dy2WXXcbPf/5zvv32W7773e+Sk5NDfn4+999/P//6179Yv349gwYNok2bNrz22msHHXvJkiX8+Mc/ZufOnbRp04bJkyfTrl07Bg4cyFlnncXChQsZMWIEH330Ea1ateL999/ntNNO4yc/+Qnf+973WLNmDY0bN+bxxx+nZ8+ejBs3jvXr17N27VratGnDs88+W6X3XqbylKatCa/Klih/6y33OnWCssX16wfz/umnwYKWLd3376/UcUWkcg4qox1DjfLiJcpnzZrl3//+972goMDz8/N9+PDh/vrrr/u0adP8pptuOrDdtm3b3L2oTHlx+/bt86ysLN+4caO7u0+dOtVvuOEGd3c/99xz/ZZbbjmw7XXXXefDhw/3vLw8d3e/7bbbfNy4ce7uPm/ePO/Vq5e7H1rKvEKfbQiVKC+frCz4wQ/gD38IWp+CHyVd4MQTYfVqeOcd6N8/7jBFJCazZ89m9uzZ9OnTB4CdO3eyatUqzjnnHO666y7uueceLr74Ys4555wyj7Ny5Uo+/vhjLrzwQiB4kl27du0OrL/iiisO2n7UqFFkhI9LePPNN/nrX/8KwHnnncfmzZv55ptvgENLmUdF5T6AO8IOu+++G5SGAtSFVuRwUJ7rhLfegkaNgufQNGoUzCfbp0IhOPfddx9Lly5l6dKlrF69mhtvvJGTTjqJJUuWcOqpp3Lffffx4IMPJj1O9+7dDxzno48+Yvbs2QfWV7QkuYVt6FGVJC9OyYLgcRYnnwxbt8Kbb4YLC7vQKlmIHN5SXKO8eInyIUOGMHHiRHbu3AnAV199xcaNG1m/fj2NGzdmzJgx3HXXXbz33nsl7l/o5JNPJjc3l0Xhs1/379/PsmXLyhXTgAEDmDJlChA8T6NNmzY0b968Su+zomp9M1ShESNg/HiYMSPsDHXuudCwIbz3XlCFNuLyvyJSBVlZKXviUWKJ8mHDhjF+/HhWrFhBVnj8pk2b8swzz7B69Wruvvtu6tSpQ7169fjjH/8IwNixYxk2bBjt2rU76AZ3/fr1mTZtGnfccQfffPMNeXl53HnnnXTv3j1pTOPGjeOGG26gZ8+eNG7cmKeeeiol77UiVKI8tHAhnH02dO4Mn30W9pK66KLgymLSJLj++pTFKiKlU4ny6KhEeQqceSa0bQuffw4ffxwu1GhuERFAyeKAjAy4+OJgekbh0MHCm9xz5kDEA15ERA5nShYJRoY1cadPDxeccEJQiXbbNqjCc3ZFpGLSpXn8cFLVz1TJIsGFFwb3tBcvhvWFlajUFCVSrRo2bMjmzZuVMFLI3dm8eTMNGzas9DHUGypB48ZBwnj55eB1880ETVG//W1wo/uXv4w7RJG016FDB3JyckjJYwfkgIYNG9KhQ4dK769kUczIkUGimD49TBbnnhsM9Fm6NLjcOOaYuEMUSWv16tWjc+fOcYchxagZqpiLLw66zc6bBzt2ELRLnXdesPLVV2ONTUQkLkoWxRx1VDC2Z9++hArlKv0hIrWckkUJRowI/h7ShXb27ITiUSIitYeSRQkKu9C+8ko4vOL444PiUdu3h4/TExGpXZQsSnDKKXDSSbBliwoLioiAkkWpSm2K0ngLEamFlCxKkTia2x0YMCAYiPHhh/DVV7HGJiJS3ZQsSpGVBW3awJo1sGwZ0KABnH9+sFJdaEWkllGyKEWZhQXVFCUitYySRRkOKSyYWIVWXWhFpBZRsihDYWHBd98NCwt26gRduwZDuxcujDs8EZFqo2RRhiZN4IILgum//z1cqNHcIlILKVkkcUhTlMZbiEgtpGSRxCWXFBUW3LmT4EHdTZrARx/BunVxhyciUi0iTRZmNtTMVprZajO7t4T1Pzaz5Wb2oZnNM7PjEtblm9nS8DWj+L7V5aijgudz790blIaiQYOitil1oRWRWiKyZGFmGcCjwDCgGzDazLoV2+x9INPdewLTgN8krNvt7r3D14io4iyPwtHch/SKUhdaEakloryy6Aesdvc17r4PmAqMTNzA3V9z913h7NtA5R/jFKHC+xZ//3tYWLAwWcydG9QyFxFJc1Emi/ZAYqN+TrisNDcCiXeNG5pZtpm9bWaXRhFgeZ1yCnTpEhQWXLgQ6NgRuncPbmKoC62I1AJRJgsrYVmJT2A3szFAJjA+YXFHd88ErgImmNkJJew3Nkwo2VE+r9dMhQVFpHaLMlnkAMcmzHcA1hffyMwuAH4CjHD3vYXL3X19+HcNMB/oU3xfd3/c3TPdPbNt27apjb6YQwoLqgutiNQiUSaLxUAXM+tsZvWBK4GDejWZWR/gTwSJYmPC8pZm1iCcbgP0B5ZHGGtSZ50FrVvDZ5/B8uVA//7QtGlQZfDLL+MMTUQkcpElC3fPA24DZgErgBfcfZmZPWhmhb2bxgNNgb8U6yLbFcg2sw+A14CH3T3WZHFIYcH69Yu60OrqQkTSXKTjLNx9pruf5O4nuPsvwmU/c/cZ4fQF7n5U8S6y7v6Wu5/q7r3Cv09GGWd5aTS3iNRWGsFdAYMHB2Py3nkHvv6ag7vQ7t1b5r4iIjWZkkUFHFJYsEMH6NEDvv024WHdIiLpR8migtQUJSK1kZJFBRUWFpw7NywsqPEWIlILKFlU0NFHwxlnBLco5swh6ELbrBmsWAFr18YdnohIJJQsKuGgwoL16gWP1AM1RYlI2lKyqIRSCwsqWYhImlKyqISuXeHEE2HzZli0iKJkMW+eutCKSFpSsqiExMKC06cD7dtDz56waxcsWBBrbCIiUVCyqCQVFhSR2kTJopIKCwuuXh10hNJ9CxFJZ0oWlVS3LgwfHkzPmAFkZUHz5vDJJ/D557HGJiKSakoWVXDQaO569YLiUaCrCxFJO0oWVZBYWHDDBjSaW0TSlpJFFTRtGhQWdA8LCw4dGqz45z9hz55YYxMRSSUliyo6qAvtMcdA796we7e60IpIWlGyqKJLLgn+zp0bVCpXU5SIpCMliypq1y4oLLhnT1hYUOMtRCQNKVmkwEFNUWeeCS1awKefwmefxRqXiEiqKFmkQGJhwXyrqy60IpJ2lCxSoFs3OOEE2LSpWGFBJQsRSRNKFilwSGHBxC60u3fHFpeISKooWaTIQYUFjzoaTjstuOv9+uvxBiYikgJKFinSvz+0agWrVsHKlagpSkTSipJFiiQWFpw+HY23EJG0omSRQgcVFjzjDGjZMqhhvmpVrHGJiFRVpMnCzIaa2UozW21m95aw/sdmttzMPjSzeWZ2XMK668xsVfi6Lso4U2XIkKCw4Ntvw782qwutiKSPyJKFmWUAjwLDgG7AaDPrVmyz94FMd+8JTAN+E+7bCngAOAPoBzxgZi2jijVVmjaF889PKCyo0dwikiaivLLoB6x29zXuvg+YCoxM3MDdX3P3XeHs20CHcHoIMMfdt7j7VmAOMDTCWFPmoC60Q4YEM/PnB8/nFhGpoaJMFu2BdQnzOeGy0twIFP4EL9e+ZjbWzLLNLDs3N7eK4aZGYWHBOXPg26ZHQd++QRfa+fNjjUtEpCqiTBZWwjIvcUOzMUAmML4i+7r74+6e6e6Zbdu2rXSgqXTMMdCvX5Af5s5FTVEikhaiTBY5wLEJ8x2A9cU3MrMLgJ8AI9x9b0X2PVwd1BSV2IXWS8yVIiKHvSiTxWKgi5l1NrP6wJXAjMQNzKwP8CeCRLExYdUsYLCZtQxvbA8Ol9UIBxUW7NsvGK23Zo260IpIjRVZsnD3POA2gi/5FcAL7r7MzB40s/C3N+OBpsBfzGypmc0I990CPESQcBYDD4bLaoTu3eH44yE3F95enFF0o1tNUSJSQ0U6zsLdZ7r7Se5+grv/Ilz2M3cvTAoXuPtR7t47fI1I2Heiu58YviZFGWeqHVJYUKO5RaSG0wjuiBw0mnvIkCCDvP66utCKSI2kZBGRs88OblV8+ims3HokZGbC3r3w2mtxhyYiUmFKFhE5pLBgYRdaNUWJSA1UZrIIxz8UTvcvtu62qIJKF+pCKyLpItmVxY8Tpn9XbN33UhxL2hkyBOrXDx61urFjJrRuDWvXhg+8EBGpOZIlCytluqR5KaZZs4TCgv/IKHrcqrrQikgNkyxZeCnTJc1LCUpsilKyEJEaxryM9nMz2wWsJriKOCGcJpw/3t2bRB5hOWVmZnp2dnbcYRxi/Xpo3x4aNYJNn2yicacjoV492Lw5qGkuIhIjM1vi7pnJtqubZH3XFMVTax1zDJx+OixeDHOXtmFEv37wzjtBF9rCErUiIoe5Mpuh3P2LxBewEzgNaBPOSzmoKUpEarpkXWf/bmY9wul2wMcEvaCeNrM7qyG+tFA4mvvllyF/SMJ4C3WhFZEaItkN7s7u/nE4fQPB0+suIXjcqbrOllOPHtC5c1BY8J28vtC2LXzxBXzySdyhiYiUS7JksT9h+nxgJoC77wAKogoq3RxUWPDlOkVdaDWaW0RqiGTJYp2Z3W5mlxHcq3gVwMwaAfWiDi6dHFRYUPctRKSGSZYsbgS6A9cDV7j7tnD5mUCNKhset3POgZYtg8HbqzoPhjp1YMEC2Lkz7tBERJJK1htqo7v/wN1HuvvshOWvuft/RR9e+kgsLPjSG63hjDNg/36YNy/ewEREyqHMcRaFT64rTeLDiiS5ESPgmWeCpqi7hw0Likb98pdw5JGQlVW9wSxaBPPnw8CB1X9uEalxkg3KywLWAc8B76B6UFUydGhQWPCtt2D7qPY0B3j3XRgwILipcfTRJe9oZXzsZa0rbf3XX8OLLwZdd+vXD65ulDBEpAzJksXRwIXAaOAq4BXgOXdfFnVg6ahZMzjvPHj1Vfh0wQYOjK/Py4O//jWeoPbtC64wlCxEpAxlJgt3zyfoAfWqmTUgSBrzzexBdy9eslzKYcSIIFn8JXcQmQ0bBl/WGRlw++1w/PFFG5Y2YK+k5RXZFoLxHY89FkzXrx80RYmIlKHMQoIAYZIYTpAoOgEzgInu/lXk0VXA4VpIsLivvoIOHYLCgltnLqLBovnx3Dc47TR4/3245x54+OHqPbeIHDbKW0gwWbmPp4C3CMZY/NzdT3f3hw63RFGTtG8fPI57926YvSML7rsvniagO8NqLQsWVP+5RaTGSTbO4hrgJOCHwFtmtj187TCz7dGHl54OKiwYl8svD26iLFoEK1bEGIiI1ATJxlnUcfdm4at5wquZuzevriDTzUGFBfNjCqJJE/jud4PpyZNjCkJEaopkVxYSgVNPhU6dYOPGoOdsbL4X1oL885+DHlkiIqVQsojBQYUF42yKysqCk0+GDRuCLloiIqVINs6iSsxsKPBbIAN4wt0fLrZ+ADAB6Alc6e7TEtblAx+Fs1+m22jxkSPhkUeCHqxbtkDHjkXrEsfRlTZd3u1K2+frr6FJE+PGQddz3Mr7YNIkuPjiyr0ZEUl7kSULM8sAHiUY1JcDLDazGe6+PGGzLwmKFN5VwiF2u3vvqOKLW72wZu8338D//m88MZjBlAbXsrrOT7CXX4ZNm6BNm3iCEZHDWpRXFv2A1e6+BsDMpgIjgQPJwt3Xhutq3bMx3nwzKDxbUBB8aZ97Lpx1VtH64sNfEufLM13WurfeCjpBucMX+4/hsy5DOXHlTJgyBX74w8q/KRFJW1Emi/YEdaUK5RA8Ya+8GppZNpAHPOzuLxXfwMzGAmMBOia249QAAwdCgwbBAO769YN6gtU13GLRIjj77KJElX/NDfDTmTBxItxxR/J6UyJS60R5g7ukb5yKPHS6Yziq8CpggpmdcMjB3B9390x3z2zbtm1l44xFVlZQv++hh6q/jl9WFkyYEEw3bw6dbr8EWrWCDz8MRnWLiBQTZbLIAY5NmO8ArC/vzu6+Pvy7BpgP9EllcIeDrBgHcN92W9CFd8sWeP6lBjBmTLBikp5pJSKHijJZLAa6mFlnM6sPXElQVyopM2sZ1qTCzNoA/Um41yFVZ1ZU8WPCBPDrbwhmpkyBPXviC0xEDkuRJQt3zwNuA2YBK4AX3H2ZmT1oZiMAzOx0M8sBRgF/MrPC0uddgWwz+wB4jeCehZJFil11FbRtG7Q8LdjeG3r3hq1bYUa5crqI1CKRDspz95nufpK7n+DuvwiX/czdZ4TTi929g7s3cffW7t49XP6Wu5/q7r3Cv09GGWdt1bAh3HJLMD1hAkUjutUUJSLFaAR3LXfLLUFvrOnTYe1ZVwUzs2dDTk7coYnIYUTJopY7+mgYPToYczHh6dZBHZKCgqBelIhISMlCDtzonjgRdl0R3uieNKn0J+2JSK2jZCH07h0MEtyxAx5fOxiOOQZWr4aFC+MOTUQOE0oWAhRdXfz20boUjLk2mJk4Mb6AROSwomQhQFBw9oQTYO1amNcxbIp64QXYuTPWuETk8KBkIQBkZARloQAeev4k6N8fvv0Wpk0re0cRqRWULOSAG24IakW98QasHZRwo1tEaj0lCzmgWTO46aZg+hervguNG8OCBcHNbhGp1ZQs5CC33x48Z+OpvzVj10X/FiycPDnWmEQkfkoWcpBOneCyy2D/fniuUVj+46mnID8/1rhEJF5KFnKIwm60984cQEHn44PSH3PnxhuUiMRKyUIO0b8/9O0LmzYb7/W8PlioG90itZqShRzCDH70o2D63uXX4Wbw0ktB+XIRqZWULKREo0ZBu3Ywb1VHtvS5APbuheeeizssEYmJkoWUqH794NGrAE8WhGMuVP5DpNZSspBSjR0bPCDpgaWXkt/sCFiyBD76KO6wRCQGShZSqjZt4NprYQ+NWNjxqmChbnSL1EpKFlKmH/4w+PuT1WFT1NNPw7598QUkIrFQspAydesGgwfDm3szyT2qO2zaBK+8EndYIlLNlCwkqaAbrfHHPeGIbjVFidQ6ShaS1ODBcMop8PtvxlCQURdmzoQNG+IOS0SqkZKFJFWnTlACJJcjeaP58KBO1NNPxx2WiFQjJQspl2uugZYt4X+2Jjznwj3eoESk2ihZSLk0bgw33wwzuYhtDY6EFSvg3XfjDktEqomShZTbrbcCdesxcd81wQKN6BapNSJNFmY21MxWmtlqM7uSuU8cAAAO90lEQVS3hPUDzOw9M8szs38rtu46M1sVvq6LMk4pnw4dgppRT3rYFDV1KuzaFW9QIlItIksWZpYBPAoMA7oBo82sW7HNvgSuB54ttm8r4AHgDKAf8ICZtYwqVim/O++E5XQnO6MfbN8OL74Yd0giUg2ivLLoB6x29zXuvg+YCoxM3MDd17r7h0BBsX2HAHPcfYu7bwXmAEMjjFXKqV8/OOsseCI/4Ua3iKS9KJNFe2BdwnxOuCxl+5rZWDPLNrPs3NzcSgcqFfOjH8FUrmSPNYR582Dt2rhDEpGIRZksrIRl5e1rWa593f1xd89098y2bdtWKDipvEsvhSM6tuBvflmw4Kmn4g1IRCIXZbLIAY5NmO8ArK+GfSVidevC7bfDRMLyH5MnQ0HxlkQRSSdRJovFQBcz62xm9YErgRnl3HcWMNjMWoY3tgeHy+QwcdNN8E7j8/iCjkEz1Pz5cYckIhGKLFm4ex5wG8GX/ArgBXdfZmYPmtkIADM73cxygFHAn8xsWbjvFuAhgoSzGHgwXCaHiRYt4Prv1eEpwl7NutEtktbM06RkQ2ZmpmdnZ8cdRq2yahUMPWkNn3EC3qgR9vXXcMQRcYclIhVgZkvcPTPZdhrBLZXWpQt0v+R4XmMgtns3vPBC3CGJSESULKRK7rwTJhGMuSh4UuU/RNKVkoVUyaBB8GmPy9lOM+q883ZQYFBE0o6ShVSJGdz84yY8zxUA+KTJ8QYkIpFQspAqGz0aXmoZNEXte/LPkJcXc0QikmpKFlJlDRvC6bdn8Qkn02DLBnj11bhDEpEUU7KQlLjl/xh/zgiuLnb+XmMuRNKNkoWkxFFHwe7LryGfOjSc8zKosKNIWlGykJS5/j+P4VWGUrdgP7ufmBJ3OCKSQkoWkjK9esHi7mFT1O8mQZpUBxARJQtJscyfX8ImWtP26w/Jz34/7nBEJEWULCSlLrqsAX8/4moA1v5MI7pF0oWShaRUnTrQ4OagKart3Gdhz56YIxKRVFCykJS7+Ke9+aBOH5rnbWXNhPI+wkREDmdKFpJyzZrBmnODq4sdv9OYC5F0oGQhkej731exl/qcun4WG7Jz4g5HRKpIyUIi0bFPa95rP4I6OB/d/ee4wxGRKlKykMg0v/N7AJywYBK7d2nMhUhNpmQhkel252A21juG4wtWM3fcm3GHIyJVoGQhkbG6GeQOvRaA/f87SQO6RWowJQuJVJdfBr2iBm97gdde3hlzNCJSWUoWEqn6PU5i3XH9acq3LP3ptLjDEZFKUrKQyLX8UXB1kfnRRD75JOZgRKRSlCwkck2/91321m3MAN7g2QdXxx2OiFSCkoVEr1kzdg8fBUCTv0xmy5aY4xGRCqsbdwBSO7T40Q0w/SmuzpvMqO/8nFtuy+C004J1Zgdvmzhf2nRFt8vOhsWL4aKLICur8u9DpLYyj7A/o5kNBX4LZABPuPvDxdY3AP4M9AU2A1e4+1oz6wSsAFaGm77t7j8o61yZmZmenZ2d2jcgqePOtrYn0mLzGobwKrMZEksYZnD99XDbbdCnz6GJRaS2MbMl7p6ZbLvIrizMLAN4FLgQyAEWm9kMd1+esNmNwFZ3P9HMrgR+DVwRrvvM3XtHFZ9UMzM+6HMD5869n8cZy0wuIrdhRxo0NDz8xnY3nHA64W/wKlxmuBctL75tAYXf/kaBB9N79hrt9n5OZz5nsl/HpEnDmTQJjj0WRowIXgMHQv361fVhiNQ8UTZD9QNWu/saADObCowEEpPFSGBcOD0N+L2Zfuulq6PO74HPheP4klt4DPYQvKrRKP7CypZn8ju/ncnrRvLoo0149FFo3hyGDYORI4O/LVpUb1wih7sok0V7YF3CfA5wRmnbuHuemX0DtA7XdTaz94HtwE/d/Y3iJzCzscBYgI4dO6Y2ekm5U3wFbnXAC3Az7Jxz4Mwzi57V7X7wdEnLKrP+vffwxdnh9QmcsvVtHuVtHmnUhA9P/A5/3D6GiV+cx/PP1+X556Fu3eBKY+TI4KpD/2mJRHjPwsxGAUPc/aZw/hqgn7vfnrDNsnCbnHD+M4Irkp1AU3ffbGZ9gZeA7u6+vbTz6Z5FDbBoEZx/PuzbF7T5zJtXPXebE89brx7ceiu89VawPJTX9mg+6Daax3aMYeL7fQ40YUFwb2PkyODVq5fuc0h6Ke89iyiTRRYwzt2HhPP3Abj7rxK2mRVus8jM6gIbgLZeLCgzmw/c5e6lZgMlixpi0SKYPz/46V6d3ZJKOu9nn8GUKfDMM7Bq1YFN87qcwkc9x/CnnVfzzJud+PbbosMcd1xwtTFyJAwYEOQekZrscEgWdYFPgfOBr4DFwFXuvixhm1uBU939B+EN7u+4+3fNrC2wxd3zzex44I1wu1J76CtZSKW5B31rn3kGnnsOcnMPrCo462yWnTaGidtHMXV2KzZsKNrtiCOCrriF9zmaN48hdpEqij1ZhEFcBEwg6Do70d1/YWYPAtnuPsPMGgJPA32ALcCV7r7GzC4HHgTygHzgAXd/uaxzKVlISuzfD3PnBonjxRdh9+5geb16+EXDWXXGGJ7ZOpy/vtKQ5QldNerVg0GDiu5zdOgQT/giFXVYJIvqpGQhKbdjB7z0UpA45s6FgoJg+RFHwKhR5Awcw/Prz2H6y3VYuLBoNUDfvkX3OXbuhNdfPzxa3tL93HrPFd9fyUIklb7+GqZODRLHe+8VLT/2WLjqKrYOH8OMNT2YPh1mzYJdu4o2MQtauurUgdNPh5Ytow9369ZgxHpBQXDefv2gVauDY0r8m2xZRfbZsgUWLID8fMjIgHPPhdatKVVFOgyUte2mTcGXZuF5Bw2CNm3K3jdVy3JzYc6conMPHgxt25Yea6rk5sLs2cF5GzWqXJ8RJQuRqKxYUXRj/Isvipb36gVjxrD70tHM+6Q906cH+WWnHuMh1SAjAx56CO67r2L7KVmIRK2gIOiC+8wz8MILwc95CH5ynncejBnD28d8h4EjmrN/fzB+41e/gq5dow9txYrgSyMvLzjvL35RdN7iQ1GSLavoPp9+CrPGLeLs/Pm8mTGQwQ9kcdJJJcdZka+fZNt++inMfajovOf/tOi8Je2bymWrV8Ovf130ef/Hf8CJJ5YdbyqsXg2/+U1w3gYNdGVRLkoWEqu9e+Ef/wgSx8svB2M6ABo2ZFuPs9iyrS6NehxPuzOOq7aQvn7nC3Yt+5zG3TsH502WCVI1vW4dBX9+GsvPxzMyqHPtNcEd/6p8O5dn26++ouDZ54rOe/110Llz8JO7bt3gVZXpJNutmLqULX9fRKuLz6LrNUm/e0tWiUE8y/+8mK2vvEXz0cM5dexZlTilkoVIPLZtg2nTgsTx+utxRyO1RSVvWsReSFCk1mrRAm66KXjdcw+MHx/8AjaDc86pvlHrb7xx6HnLU9e9qtPr1sFTTxXd7b3++qAjQPFtS9q/rGXJtv3yS3jyyaBNJiMDrrkG2rUL5vPzg7+VmS7Pdlu3Br3nCjVpErwqojI/3Hft4sCo0X37gjv8Ef33pWQhEqVLL4Xf/a6oxMnDD1d/iZPqPG+hG2+Mpx/ptdfG13818fOeMyeef+eBAyM7lZqhRKJWUzvgS8XU0H9n3bMQEZGkypss9AxuERFJSslCRESSUrIQEZGklCxERCQpJQsREUlKyUJERJJKm66zZpYLfJF0w9K1ATalKJyaora959r2fkHvubaoyns+zt2TFlRPm2RRVWaWXZ6+xumktr3n2vZ+Qe+5tqiO96xmKBERSUrJQkREklKyKPJ43AHEoLa959r2fkHvubaI/D3rnoWIiCSlKwsREUlKyUJERJKq9cnCzIaa2UozW21m98YdT9TM7Fgze83MVpjZMjP7YdwxVRczyzCz983s73HHUh3MrIWZTTOzT8J/77R/qIWZ/Sj87/pjM3vOzBrGHVOqmdlEM9toZh8nLGtlZnPMbFX4t2Wqz1urk4WZZQCPAsOAbsBoM+sWb1SRywP+3d27AmcCt9aC91zoh8CKuIOoRr8FXnX3U4BepPl7N7P2wB1Aprv3ADKAK+ONKhKTgaHFlt0LzHP3LsC8cD6lanWyAPoBq919jbvvA6YCI2OOKVLu/rW7vxdO7yD4Amkfb1TRM7MOwHDgibhjqQ5m1hwYADwJ4O773H1bvFFVi7pAIzOrCzQG1sccT8q5+wJgS7HFI4GnwumngEtTfd7anizaA+sS5nOoBV+chcysE9AHeCfeSKrFBOA/gIK4A6kmxwO5wKSw6e0JM2sSd1BRcvevgP8CvgS+Br5x99nxRlVtjnL3ryH4QQgcmeoT1PZkYSUsqxV9ic2sKfBX4E533x53PFEys4uBje6+JO5YqlFd4DTgj+7eB/iWCJomDidhO/1IoDNwDNDEzMbEG1X6qO3JIgc4NmG+A2l42VqcmdUjSBRT3P1vccdTDfoDI8xsLUFT43lm9ky8IUUuB8hx98KrxmkEySOdXQB87u657r4f+BtwVswxVZd/mVk7gPDvxlSfoLYni8VAFzPrbGb1CW6GzYg5pkiZmRG0Y69w9/+JO57q4O73uXsHd+9E8G/8T3dP61+c7r4BWGdmJ4eLzgeWxxhSdfgSONPMGof/nZ9Pmt/UTzADuC6cvg6YnuoT1E31AWsSd88zs9uAWQQ9Jya6+7KYw4paf+Aa4CMzWxou+093nxljTBKN24Ep4Q+hNcANMccTKXd/x8ymAe8R9Pp7nzQs/WFmzwEDgTZmlgM8ADwMvGBmNxIkzVEpP6/KfYiISDK1vRlKRETKQclCRESSUrIQEZGklCxERCQpJQsREUlKyUJERJJSshARkaRq9aA8kaiZ2VSCGmSdgKOB/+Pur8QalEgl6MpCJFq9gDXufgZwNcFoW5EaRyO4RSJiZo0ISi8c6+57zKwV8E74gBqRGkVXFiLR6QGscvc94fxpwAcxxiNSabpnIRKdXkDH8DnQGcDPCR7AJFLjKFmIRKcXMAWYDzQHfunuC2ONSKSSdM9CJCJmtgD4vruvjDsWkapSshCJiJl9RXBzu7Y891vSmJKFiIgkpd5QIiKSlJKFiIgkpWQhIiJJKVmIiEhSShYiIpKUkoWIiCSlZCEiIkn9f8g45TWEKd5rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_MSE(df_train.X,df_train.Y,df_test.X,df_test.Y,p=11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$p=5$ is the optimal value because it has the lowest test error, performance does not  \n",
    "improve beyond it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standardize(X,y,mean=None,std=None):\n",
    "    '''\n",
    "    Standardize X,\n",
    "    Centralize y\n",
    "    '''\n",
    "    if(mean and std): # standardize test data with mean and std of Dtrain\n",
    "        X_std = (X-mean)/std\n",
    "        mean_y = np.mean(y)\n",
    "    else:   \n",
    "        mean_x = np.mean(X, axis=0)\n",
    "        std_x = np.std(X, axis=0)\n",
    "        mean_y = np.mean(y)\n",
    "\n",
    "        X_std = (X-mean_x)/std_x\n",
    "    \n",
    "    return (X_std, mean_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "def least_square_est_ridge(X,y,lmd):\n",
    "    '''\n",
    "    X:Design matrix\n",
    "    y:True response\n",
    "    lmd:Lambda,regularization parameter\n",
    "    '''\n",
    "        \n",
    "    X = X if isinstance(X,np.ndarray) else np.asarray(X)\n",
    "    y = y if isinstance(y,np.ndarray) else np.asarray(y)\n",
    "    \n",
    "    #dont standardize the first column of vandermond(all 1s)\n",
    "    b = np.ones((X.shape[0],X.shape[1]))\n",
    "    b[:,1:] = standardize(X[:,1:],y)[0]\n",
    "    X = b\n",
    "       \n",
    "    b_hat = inv(X.T@X+np.identity(X.shape[1])*lmd)@X.T@y\n",
    "    return b_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standardize_test(X,y,mean,std):\n",
    "        \n",
    "    X = X if isinstance(X,np.ndarray) else np.asarray(X)\n",
    "    y = y if isinstance(y,np.ndarray) else np.asarray(y)\n",
    "    \n",
    "    #dont standardize the first column of vandermond(all 1s)\n",
    "    b = np.ones((X.shape[0],X.shape[1]))\n",
    "    b[:,1:] = standardize(X[:,1:],y,mean,std)[0]\n",
    "    X = b\n",
    "    y = standardize(X[:,1:],y,mean,std)[1]\n",
    "    \n",
    "    return (X,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_MSE_ridge(X_tr,y_tr,X_te,y_te,p):\n",
    "    \n",
    "    y_hat = y = standardize(X_tr,y_tr)[1] # mean of response on Dtrain\n",
    "    lmd = range(101) #lambda\n",
    "    training_errors = []\n",
    "    test_errors = []\n",
    "    \n",
    "    for i in lmd:\n",
    "        X_des = X_design(X_tr,p) #vandermond train data\n",
    "        X_des_te = X_design(X_te,p) #vandermond test data\n",
    "        beta_hat = least_square_est_ridge(X_des,y_tr,i) #beta hat \n",
    "        \n",
    "        training_errors.append(MSE(predict(X_des,beta_hat,y_hat),y_tr)) #train error\n",
    "        \n",
    "    \n",
    "#         X_des_te = standardize_test(X_des_te, y_tr, mean_tr, std_tr)\n",
    "        \n",
    "        test_errors.append(MSE(predict(X_des_te,beta_hat,y_hat),y_te))  #test error\n",
    "        \n",
    "    plt.plot(training_errors, marker='.', color='blue', linewidth=2, label='train error')\n",
    "    plt.plot(test_errors, marker='.', color='red', linewidth=2, label='test error')\n",
    "    # Add legend\n",
    "    plt.legend()\n",
    "    \n",
    "    # Add titles\n",
    "    plt.title(\"Train&Test Errors\")\n",
    "    plt.ylabel(\"MSE\")\n",
    "    plt.xlabel(r'$\\lambda$')\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEYCAYAAABRB/GsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9//HXJzuRTTZFQMHWIrJUBZGAaFgFiohV9GqtS23p/bXXrbdutVqt97a29rpdWy3XulSrxWLrSq2CCWtQoVpEXEBAQdSEKJsYSDKf3x/nJAwhO5nMZOb9fDzmkTmznPM5Mzpvvt/vOd9j7o6IiKSutHgXICIi8aUgEBFJcQoCEZEUpyAQEUlxCgIRkRSnIBARSXEKAmmzzCzdzHaa2eHxrkWkLVMQSKsJf7SrbhEz+zJq+VtNXZ+7V7p7e3f/sJHbNzO708xKw9ufazz/YlQ95Wa2J2r5nqbWF7Xe35jZfQ28ZouZ7arxGf2qudsUaYqMeBcgqcPd21fdN7MNwHfdfV5drzezDHevaMESpgDnAEOAUuCkGvVNjNr2o8Bad7+pBbffkLHuvqyhF9X2uTT1s4rBZyttmFoEkjDM7L/MbLaZPW5mO4DzzSzPzJaZ2VYz+9jM7jazzPD1GWbmZtY3XH40fP7vZrbDzIrMrF/UJiqAXcAn7l5WXwjVU+OZZvZmWM8CMzs66rmbwhq3m9nbZjbSzM4ELgMuCf+Vv7QZ2/wPM3vJzO4zs8+BH9fxWHr4GW40s0/M7H4zax+uY5CZlZnZv5vZJuAZM2tvZk+Y2Wdm9nn4OXdqan3S9ikIJNGcATwGdAJmE/x4Xw50A0YBk4Dv1/P+84AbgC7Ah8AtUc+9BRwC/N7MrKmFmdlJwP8CFwJdw/r+Fv4ADwW+RdDa6ARMBT5y9yeBu4E/hN1YI5u63dAY4FWCz+GuOh77IcHnNwroD/QCbotaRxZwHHAUcCbB5+jAYUB3gsDa08z6pA1TEEiiWezuz7p7xN2/dPfX3P0Vd69w93XALOCUet4/x92Xu3s58CfgWAAzywL+AXwP6AncVxUGZvaKmU1uRG3fB+5y93+G4xO/AzqG26gAcoGBQLq7v+/uHzRx318KWxpVt/OinnvP3R8It/tlHY99C/iVu3/o7tuAnwLnR63DgBvCz/VLoJwgAI4MP99Xo9YtKURBIIlmY/SCmR1tZs+HXR3bgZ8T/Au4Lp9E3d8FVI1LTADaufvjwFnAAIIw6AwcCSxpRG1HADdG/1gDBwO93P1fBD+8vwSKzewRM+veiHVGm+DunaNuj0U9t7GW19d87DAgOnw+ANpHdffscffiqOdnAUsJWjUbw24l/SakIH3pkmhqTof7e2AV8FV37wjcSPAv26bKIPhXO+G/ek8DTgBeAe539+2NWMdG4Cc1fqxz3f2ZcL0Punse8BWgA0Fo1bZPzVHbOmo+tpkgrKocDuwMWwf7vT4cJ/mpu/cn6GY6D5jRArVKG6MgkETXAdgGfGFmA6h/fKA+C4GOZvYzM2sXPlYIfA2INHIdvweuMLPjw0NRO5jZdDPLCQdjTzazbIKWSBlQGb7vU6Bfc8Ylmuhx4Coz621mHQnGR/5U14vNbIKZDQhbAdsJgrKyrtdL8lIQSKL7T4LB2R0EP8Szm7MSd/8cmAiMJviX87+ATGAo8H0zu7gR61gIXAncD2wF3iU4HNWBdsAdBIelfgxkAzeFb30M6Ax8ZmaL6tnEyzXOI3ikibt5D/AsQStnDUEAXVXP6/uEr99B8Hk8Dfy1iduUJGC6MI2ISGpTi0BEJMUpCEREUpyCQEQkxSkIRERSXJuYdK5bt27et2/feJchItKmrFixYou7N3hiY5sIgr59+7J8+fJ4lyEi0qaYWaOmOVHXkIhIilMQiIikOAWBiEiKaxNjBCKSHMrLy9m0aRNlZWXxLiWp5OTk0Lt3bzIzM5v1fgWBiLSaTZs20aFDB/r27Uvs5+BLDe5OaWkpmzZtol+/fg2/oRbqGhKRVlNWVkbXrl0VAi3IzOjatesBtbKSOgiKiuCXvwz+ikhiUAi0vAP9TJO2a2jxYjjlFIhEoF07mD8f8vLiXZWISOJJ2hbBokVBCADs2QOFhXEtR0QSwNatW/nd737XrPdOmTKFrVu3tnBFiSFpgyA/H6paS5mZwbKIpLb6gqCysv6Ls82dO5fOnTu3aD0VFRX1LteloVqbKmmDIG94JQ91uZJ5jOOxS4vULSTSRrXkWN+1117L+++/z7HHHstVV11FYWEhY8aM4bzzzmPw4MEATJ8+naFDhzJw4EBmzZpV/d6+ffuyZcsWNmzYwIABA/je977HwIEDmThxIl9++eV+2yopKeHMM8/khBNO4IQTTmDJkiUA3HTTTcycOZOJEydywQUX8NBDDzFjxgxOO+00Jk6ciLtz1VVXMWjQIAYPHszs2cFF+WqrtcW4e8Lfhg4d6k22ZIlHwB28PKud+9KlTV+HiLSo1atXV98P//ds8Vt91q9f7wMHDqxeLigo8NzcXF+3bl31Y6Wlpe7uvmvXLh84cKBv2bLF3d2POOIILykp8fXr13t6erq//vrr7u4+Y8YMf+SRR/bb1rnnnuuLFi1yd/cPPvjAjz76aHd3/9nPfubHH3+879q1y93dH3zwQe/Vq1f1dufMmePjx4/3iooK/+STT7xPnz6+efPmWmut67Pd+xmz3BvxG5u0g8UsWFB9N608HCRQs0BEahg+fPg+x9/ffffd/O1vfwNg48aNrFmzhq5du+7znn79+nHssccCMHToUDZs2LDfeufNm8fq1aurl7dv386OHTsAmDZtGu3atat+bsKECXTp0gWAxYsXc+6555Kens4hhxzCKaecwmuvvUbHjh33q7WlJG8Q5OfjGIYTScsgTYMEIgmlMZdLLyqCceOCAz6ysmJz9N9BBx1Ufb+wsJB58+ZRVFREbm4u+fn5tR6fn52dXX0/PT291q6hSCRCUVHRPj/4tW2z5rLX88HUfF9LSdoxAvLy+PSIEwF49sT/UmtApA3Kywt+/G+5pWVCoEOHDtX/Kq/Ntm3bOPjgg8nNzeWdd95h2bJlzd7WxIkTueeee6qX33jjjUa97+STT2b27NlUVlZSUlLCwoULGT58eLPraIzkDQKg7Ij+AHyyp0ucKxGR5srLg+uua5l/y3Xt2pVRo0YxaNAgrrrqqv2enzRpEhUVFQwZMoQbbriBESNGNHtbd999N8uXL2fIkCEcc8wx3HfffY163xlnnMGQIUP4+te/ztixY/n1r3/NoYce2uw6GsPqa4YkimHDhnlzLkzzwTlXc8QTtzGr3y+Zue7aGFQmIk3x9ttvM2DAgHiXkZRq+2zNbIW7D2vovUndIsjq3QOA7G3Fca5ERCRxJXUQ5PYNgiD3CwWBiEhdkjoI2h8ZBEGn3cU08oQ9EZGUk9RBkN4zCIIeFFNaGudiREQSVFIHAT32BkGxeodERGqV3EHQvXvwhxJKihP/6CgRkXhI7iDIzuaLzE5kUsHn65Nz+lgRabwDmYYa4M4772TXrl0tWFFiiGkQmNkGM3vTzN4ws+XhY13M7CUzWxP+PTiWNXyRG7QKvlivviGRVBfvIGjutNONfV1ztcZcQ2PcfUvU8rXAfHe/1cyuDZevidXGd3fsAdvWsntjMdA/VpsRkVgpKgomjczPP+DTi6OnoZ4wYQK33XYbt912G0888QS7d+/mjDPO4Oabb+aLL77g7LPPZtOmTVRWVnLDDTfw6aefsnnzZsaMGUO3bt0oKCjYZ90rVqzgRz/6ETt37qRbt2489NBD9OzZk/z8fEaOHMmSJUuYNm0ab775Jl26dOH111/n+OOP5/rrr+c73/kO69atIzc3l1mzZjFkyBBuuukmNm/ezIYNG+jWrRuPPfbYAe17feIx6dzpQH54/2GgkBgGQUWXHrARKj9Wi0AkocTq2sX1zJZw6623smrVqup5f1588UXWrFnDq6++irszbdo0Fi5cSElJCYcddhjPP/88EMxB1KlTJ26//XYKCgro1q3bPustLy/n0ksv5emnn6Z79+7Mnj2b66+/ngceeAAIWiILwhmRL7roIt577z3mzZtHeno6l156KccddxxPPfUUL7/8MhdccEF1fStWrGDx4sW1TlzXkmIdBA68aGYO/N7dZwGHuPvHAO7+sZn1qO2NZjYTmAlw+OGHN7+C8MghHTYkIjW9+OKLvPjiixx33HEA7Ny5kzVr1jB69Gh+/OMfc8011zB16lRGjx5d73reffddVq1axYQJE4DgCmI9e/asfv6cc87Z5/UzZswgPT0dCKadfvLJJwEYO3YspaWlbNu2Ddh/uupYiXUQjHL3zeGP/Utm9k5j3xiGxiwI5hpqbgHphwVBkPG5gkAkoSTAPNTuznXXXcf3v//9/Z5bsWIFc+fO5brrrmPixInceOON9a5n4MCBFNVxGbWmTjttYWspVtNO1xTTwWJ33xz+LQb+BgwHPjWzngDh35j+QmeH8w3lbFcQiLQ5LTwPdc1pqE899VQeeOABdu7cCcBHH31EcXExmzdvJjc3l/PPP58f//jH/POf/6z1/VX69+9PSUlJdRCUl5fz1ltvNaqmk08+mT/96U9AcD2Ebt260bFjxwPaz6aKWYvAzA4C0tx9R3h/IvBz4BngQuDW8O/TsaoB4KB+QRC013xDIm1TXl6LtQKip6GePHkyt912G2+//TZ54frbt2/Po48+ytq1a7nqqqtIS0sjMzOTe++9F4CZM2cyefJkevbsuc9gcVZWFnPmzOGyyy5j27ZtVFRUcMUVVzBw4MAGa7rpppu4+OKLGTJkCLm5uTz88MMtsq9NEbNpqM3sSIJWAASB85i7/7eZdQWeAA4HPgRmuPtn9a2rudNQA/jLBdi4sRRyCiO+LCQnp1mrEZEWoGmoY+dApqGOWYvA3dcBX6/l8VJgXKy2W5MdsneaiZIS6NOntbYsItI2JPeZxVA9zYTmGxIRqV3yB0HXrkQwulHKlk80F7VIvLWFqyK2NQf6mSZ/EKSnsyM7OPlj2/tbGnixiMRSTk4OpaWlCoMW5O6UlpaScwADoPE4s7jV7Wrfg067S9i1oRiI7UWgRaRuvXv3ZtOmTZSUlMS7lKSSk5ND7969m/3+lAiCPZ16QOlb7NmkQQKReMrMzKRfv37xLkNqSP6uIaCya3DkUOQTBYGISE0pEQRVh5BaiYJARKSmlAiCzF5BEGRtVRCIiNSUEkGQc3g439AODVCJiNSUEkHQ/sggCA7aVczSpXEuRkQkwaREEKzdFl7E3osZPz6Y2VZERAIpEQSL39s739CePcFV70REJJASQXD8pCAIerOJkVZEfn586xERSSQpEQTDD3oLB7LZwzzGkYf6hkREqqREEBBeNBogvVJ9QyIi0VIjCPLzwYJdjVgG6hsSEdkrNYIgL4+teZMAeLTb5S168WsRkbYuNYIAyBx1IgBbt6fMLouINErK/CoeNCiY8fDQsg1s2xbnYkREEkjKBIH16wtAP9azfn08KxERSSwpEwSEc6D3Yz0bNsS3FBGRRJI6QXDYYZSnZXEIxWx854t4VyMikjBSJwjS0tjR5QgAdry5Ib61iIgkkNQJAqC8V9A9VL5mQ3wLERFJICkVBGlf6QtA5iaNFouIVEmpIKg6hLR9yXrc41yMiEiCSKkgaDcgCILeFevZsiXOxYiIJIiUCgI7UoeQiojUlFJBEH0ugU4qExEJpFYQdOvG7oxcOrONzau3xrsaEZGEkFpBYMbObkGr4ItVahKIiECqBQFQ3icIgsj7CgIREUjBIMg8KggCnUsgIhJIuSBoPzgIgg6l61myJM7FiIgkgJQLgvUEQXCEb2D8eCjSdexFJMWlXBAs+agvEBxCukfXsRcRSb0gGDwtaBEcxRpG2lJdx15EUl7KBcHw9qtxIIty5vk48lDfkIiktpgHgZmlm9nrZvZcuNzPzF4xszVmNtvMsmJdwz6i+oIyIuobEhFpjRbB5cDbUcu/Au5w96OAz4FLWqGGvfLzIT0dgErSqTgpv1U3LyKSaGIaBGbWG/gGcH+4bMBYYE74koeB6bGsYT95edjllwPwD07l3S55rbp5EZFEE+sWwZ3A1UAkXO4KbHX3inB5E9Crtjea2UwzW25my0tKSlq2qulB9hzKJ6xc2bKrFhFpa2IWBGY2FSh29xXRD9fy0lovEePus9x9mLsP6969e8sWN3gwAINYxap/VbbsukVE2phYtghGAdPMbAPwZ4IuoTuBzmaWEb6mN7A5hjXUrnNnvuh2OO0oY8uyta2+eRGRRBKzIHD369y9t7v3Bf4NeNndvwUUAGeFL7sQeDpWNdRb36AhAKStUt+QiKS2eJxHcA3wIzNbSzBm8Ic41EDuiCAIDitdyVZdmkBEUlhGwy85cO5eCBSG99cBw1tju/VJO+7rAAxhJW++CaNHx7kgEZE4Sbkzi6sNCVoEQ1ipI4dEJKWlbhB89atUZObQjw28t3xbvKsREYmb1A2CjAx29RsIQNlrq+JcjIhI/KRuEACZQ4Puoax3VuoiNSKSslI6CD7pEQTBMZUrdZEaEUlZKR0ERTv3DhjrIjUikqpSOgiO+mYw1cRQljPSl+giNSKSklI6CE44eC0O5LCHl3y8LlIjIikppYMgui8okz3sfL6wzpeKiCSr1A6C/HwsMxMAx1jRPj++9YiIxEFqB0FeHvzudwBsohdPF+siNSKSelI7CAAuvJDKnFz68iGrXi6OdzUiIq1OQZCZiY8IWgKd3lzEzp1xrkdEpJUpCICM/GDq0VGRRTqpTERSjoIAquegHs0iFi2Kcy0iIq1MQQAwYgSR9AyO5Q1WFGyPdzUiIq1KQQCQm0vlccNIJwJLl7JwYbwLEhFpPQqCUHH/oHtoZGQREydqAjoRSR0KgtDS9JOBYJxAE9CJSCqpNwjM7Pyo+6NqPPcfsSoqHo44bxQRYCRLOdkLNQGdiKSMhloEP4q6/781nvtOC9cSV8M7voOZkUElLzCJXh+qb0hEUkNDQWB13K9tuW0rLKzeoSz28MHDhfGsRkSk1TQUBF7H/dqW27b8fMjKql58fPMp8atFRKQVNRQER5vZSjN7M+p+1XL/Vqiv9eTlwbx5RA46iDScJSs7UFoa76JERGIvo4HnB7RKFYnipJNIO+cceOABpvozzJ07mG9/O95FiYjEVr0tAnf/IPoG7ASOB7qFy8nn9NODPzzNM8/EuRYRkVbQ0OGjz5nZoPB+T2AVwdFCj5jZFa1QX+sbP55ITjuG8xqvPb2ZBQviXZCISGw1NEbQz91XhfcvBl5y99OAE0myw0er5eayddgEAE4tf1ZnGYtI0msoCMqj7o8D5gK4+w4gEqui4q2oR9A9dCW3c/yeIp1lLCJJraEg2Ghml5rZGQRjAy8AmFk7IDPWxcXLkXmH4EB/3mM+4zglS00CEUleDQXBJcBA4CLgHHffGj4+AngwhnXF1YDylUBwxlw2u/GCwrjWIyISS/UePuruxcC/1/J4AVAQq6LiLj8fy8yE8qBn7L538hnpYMl1LrWICNBAEJhZvQdQuvu0li0nQeTlwXPP4d/4BlZRwaL3e/L663D88fEuTESk5TV0QlkesBF4HHiFZJtfqD4TJ2Jnn4099hgX8jAPPvgzBYGIJKWGxggOBX4CDALuAiYAW9x9gbsn/xH2F18c/OFBHrg/onMKRCQpNXRmcaW7v+DuFxIMEK8FCs3s0lapLt7GjmX3IYfTlw84saxQ5xSISFJq8AplZpZtZt8EHgV+CNwN/DXWhSWEtDReHRi0Cv6HH3H8niIKkneIXERSVENTTDwMLCU4h+Bmdz/B3W9x948aWrGZ5ZjZq2b2LzN7y8xuDh/vZ2avmNkaM5ttZlkNrSueeowfAsCx/Iv5jOMrxWoSiEhyaahF8G3ga8DlwFIz2x7edpjZ9gbeuxsY6+5fB44FJpnZCOBXwB3ufhTwOcG5CgmrP+/iVJ1TUMbO5wrjXJGISMtqaIwgzd07hLeOUbcO7t6xgfe6u+8MFzPDmwNjgTnh4w8D0w9wH2IrPx8LL1hjwEPvj2LFiviWJCLSkhocIzgQZpZuZm8AxcBLwPvAVnevCF+yCehVx3tnmtlyM1teUlISyzLrl5cHL78MPXqQhnMk67n4Yg0ai0jyiGkQhEcdHQv0BoZT+4Vuar3kpbvPcvdh7j6se/fusSyzYaNGwa9/DcDV/JpVb0YYO1ZhICLJIaZBUCWco6iQ4BDUzmZWdSJbb2Bza9RwwM49l+2dejOQ1XyD59m9G81KKiJJIWZBYGbdzaxzeL8dMB54m2COorPCl10IPB2rGlpUVhafXXAlAHdzKSd6ET17xrkmEZEWEMsWQU+gwMxWAq8RXNTmOeAa4EdmthboCvwhhjW0qL7ThuBAPz6ggDGsuKcIr7VjS0Sk7WhorqFmc/eVwHG1PL6OYLyg7XntNcwM3MlmN+1XFPDCC3lMnhzvwkREmq9VxgiSRn4+5OQAwaGkhnPxxbB4cVyrEhE5IAqCpsjLg/nz4fTgUpbf4UE++3SPjiASkTZNQdBUeXkwZw5buh/NV3mfJzmToeVFPFPvlRtERBKXgqA5MjLYefZ3AJjKc8xnHNteUJNARNomBUEz9T2svHoOohzK6PRGIX9NjTlZRSTJKAiaa8wYLBw4TsP5nM5ceCG8+GKc6xIRaSIFQXNVzUE0ZgwAV3In5TvLmDIFli6Nc20iIk2gIDgQeXnw97+zpfsA+vMef2cyJ1QW8ZvfxLswEZHGUxAcqOxstl18BQ7kU0gBY9jybBFvvRXvwkREGkdB0AK+0rkUsPDiNbsZW/EPxo5Fl7UUkTZBQdAS8vOxdnvPOM5jKcXFzoQJGi8QkcSnIGgJVWccX3kl5Rk5nMpLzOEsTqgs4sYb412ciEj9zNvA9JnDhg3z5cuXx7uMRtn4nRvp8+AtOFBGNmMpYMrP88jICKYqysuLd4UikirMbIW7D2vodTGbfTRV9TmqXdSJZruZynP89MY80tIgOztoOCgMRCSRqGuopeXnY+3aAUEY/FvOU+TwJZEI7Nmjq5qJSOJRELS0qvGCn/wEDj2Ur5St5jmmch2/4ITKIsrL4Ze/1GylIpI4NEYQS6tWwYknwq5dRDDKyGE883klLU/dRCISc40dI1CLIJYGDYKzzwaC+YhyKOMUCtVNJCIJRUEQazNnQlYWEITBUJZzLb/khMoivvxS3UQiEn/qGmoNRUVwxx3wl78AVHcTjWM+r6qbSERiRF1DiSQvD554AqZPB/Z2E+VTQCQCu3erm0hE4kdB0JquvnqfbqIJzOM6fsHwSBHLlsEvfqFuIhFpfeoaam1FRfDb38Kf/wyVlTjwJe0Yx3yWkUe7duomEpGWoTOLE1VeXnDLzoYHHsCAdnzJ97mXfAop/DKfu+7Ko7BQU1KISOtQiyBeiopg7FgoKwPAgUrS2EO2BpFFpEVosDjRVV3q8ic/gR49MCCDCNnsJj8816CsDP7v/3SIqYjElrqG4qmqm2jCBBg/HiorSSfCN3ieDMqZ5xN48EFNWCcisaUWQSLIz4dFi+Cb34SMDE5iCT/nZxQwhhEUVbcO/vAHtQ5EpOVpjCDRXHop3HMPEIwbFDGC5ziNAsawDLUORKTxNEbQVp13HrRrBxZcA3kky/hvrqeAsfu0Du69V60DEWkZGiNINFXTWBcWwsKF8MILGJBNGXdwJc8ylZd9HI88kocZ5OSodSAiB0ZdQ4ms6hDT3bsh/J4c2EM2+RSwjODXf8KEYJhhzBgFgojs1diuIQVBoisqCloHS5bA889XP/weRzGHs3iW06oDITsbCgoUBiISUBAkm6IiGDcuGCCIah1UWgaX+R10YgeF5LOxVx6nnw7nn69AEEl1CoJkVNU6WLUKHnus+mEnmNp6dzi1ddXRRVdfDR07aqoKkVSlIEhmVa2D3bshEql+2IF5jKOAsdWHm0LQZfTyyzByZJzqFZG4UBAku6rWQZcucNllwbUvQ8GAchaXchdd+ZxC8lnbLY/x4+EHP4DRo+NWtYi0IgVBKikqCo4hnTcPX7AACx+urcvILDiB+ZhjYPJkdRmJJLO4B4GZ9QH+CBwKRIBZ7n6XmXUBZgN9gQ3A2e7+eX3rUhA0Uj1dRss4kRc5lReYVN1llJ4O114LBx2kcQSRZJQIQdAT6Onu/zSzDsAKYDpwEfCZu99qZtcCB7v7NfWtS0HQBFVdRp07wxVX7NdlVEEGN3EjETIoJL86FDIz4X/+B3buVCiIJIu4B8F+GzJ7GrgnvOW7+8dhWBS6e//63qsgaKaiouDEgiVL8Llzq7uMIGiildcYR6jqOsrKgrvvhtJShYJIW5ZQQWBmfYGFwCDgQ3fvHPXc5+5+cC3vmQnMBDj88MOHfvDBBzGvM2lVdRnt2QORCO6+3zjCHrIZy8vVLQQAs71HHCkMRNqehJl0zszaA08CV7j79sa+z91nufswdx/WvXv32BWYCqrmL7rlFrjvPiwnJ/iVBwxIx8mhjNmcw6N8i9EsAILz1srKgsHl00+HuXPjuA8iEjMxbRGYWSbwHPAPd789fOxd1DUUX9GHnl5+edBSiPrvwIGFnMJqBvBHvs0y9p6AMGBAcPve92DSpNYvXUQaL+5dQ2ZmwMMEA8NXRD1+G1AaNVjcxd2vrm9dCoIYqgqFd9+FP/5xn0AAcIy5TOI9+vMEZ+/TdXT00cHtkktg6tRWrltEGpQIQXASsAh4k2BsEuAnwCvAE8DhwIfADHf/rL51KQhaQfQ4AkBl5X4viWA8y1Te56v8hRn7hMKRR8LXvgYXXABHHAELFmigWSTe4h4ELUlB0EqqWgddu+5z6KlXVu5zxBEEofA8U1hPPx7n3H26j8LhB7Ky4PbbYds2hYJIPCgI5MA0MRQKGMM6+vFHLmQx+89hkZkJN98c3FcoiLQOBYG0nCaEggMrGMp6+vEgF/F3vrHf6jIy4Be/gOOOg9deUzCIxIqCQGKjCaEAsIEjWMeRzOYc3uIYRrN4nzOaIQiGyy+HTp1g/HiFgkhLURBI7NUWCmbBQHMt/105VdNcZPJFR9pSAAALrUlEQVSDDo/y1o7e5LNgn2BISwuOQDr0UDjvvGCcobBQrQaR5lAQSOtqYihAMLYATjlZ/JSfk0Fkv9ZC9MDzbbdpLiSRplAQSPzUFgrp6cFzFRXBr3tGRjBLahQHKknnIS6glO48xXQA8incJyAyMuDii4Pz4U4/PXivWg0i+1MQSGKoCoX8/GC56n4kss95C9HzH1WJAE4a4Owhm8u5k658tl+rIS0taHRkZQUzqG7frlAQAQWBtAW1tRzc8Uik1oHn6DGGu7iMnXTgRSbuEwpVMjLg298OBqDPOisIC7UaJNUoCKRtqS0Uqv6pX+/gs7GI0RTTgzmcyUb6cAoL6xxryMyEG28MlseMUShIclMQSNtVW3dS1657J8iDBlsNlWQwi++ync48wzRg/7GGtLSgd+rgg2H6dOjTBxYtUqtBkoeCQJJPHecwEInUeWQS7B1rMJxyMvkpt5BBBYWMqbVbKT0dzjkHOnaEGTOgXTt1K0nbpCCQ5NbQkUnU3WqAqgvypPEPTuVzDuYppvMRh9XZreQejDt897vQvn3QgtC4gyQ6BYGkjrq6kqIDIhKpnlG1tiOUIAgHCK7r/DjnsJOO/IWzKKPdft1KsHfcISMDfvADyM0NToYzU0BIYlAQiDQUEOydGsOhnjGH4BWVZHA/l7CLg/grZxAhvcGA+OEPIScHTjtNASGtT0EgUpcmdCs1JiAqyOT3zKSMdvyN6ThpDQbEzJnB9aCnTw+WFRASCwoCkcZoTLcSNGrcAfYPiPu5JAyIM6ggo86AqBqDuPrq4MI+GzfClCl7S1JISHMoCEQORGMCouocB5oWEJVk8Ajfpowc/sZ0dtKh1oCAvSGRmQkXXbS3FZGbq4CQhikIRGKhhQMC9g5SV5LO4x3/nTe39+FQPuEpplNOVp0hAcHmpk4NDnGdOBF69YIVK/af0UNhkZoUBCKtqZldTHWNQcDek+MgOEHuSb5JGTn8ncls5jBOYgmFBNurLSyqTszOyIDLLgtKmDZNh72mEgWBSCJoSkCkp0P//viqVY0KCQim8nYMC8ckbuUaPDw/AmoPiOgxidNPD7qbxo+Hww4LWhNjxuwtVWHRtikIRBJZXbOywt5ZWetpRUDdIQH7j0n8iXMpJ5u5TKaYHpzMoka3Ji65JCjltNOgQwdYsEAB0VYoCETaqsa0Isyqr+0QcbDKyiAZ6jhZLlp0a4KoI5zu4AqcNJ5lap2HwFZJS4OTTgrGJvLygivKrVsXjFfocNjEoSAQSTb1tSIaOWhdX3dTlejWRIR0nmcKe8imgHw+4VAG8DYvMw6ovTVRJT096HLKyoLRo6FHD1i7FiZP1sl1rUVBIJKKmtKaSEsjEvEmtSYgehB7b9fTE5xFJZnMZyyfcgjH8TqFBIMN9Y1TpKcHRztlZQUtjO7dYc2a4BwKhcWBUxCIyF5NbU1As8Ymqjh7WxYWhsXjnEOEDOYxnmK6M5QV9YZFlarZYLt3h+LioOxu3WD16mA4JXp3FBj7UhCISNM0pTWRnk4k4lARtCbMWyYsgpZFBgWMoYTuDGElBYxlGXmMoKjB1sWkScHJdyNHBqVXtS6ixy2qdi0VgkNBICIto6mtiRYPC2Nz1uEcumcTaUSoIINH+RaOMZ9xbKFbo1sXZvsePnvJJcHjEyZAly7BriZTWCgIRKT1xDEsoPZxi2eZSiUZLGYUpXSjP+8wlymAVYcFNHyuRdVJeFWtjHfeCVoZ6emJ38pQEIhIYmluWGRmwmWXEbnjTiivOODAqFJ1Ml4aTkU4xbgBLzGez+hKHkX1hkWVmq2M888P7p9yStDKWLly/7GM6PuxDA4FgYi0PXWFRV5e3dOH12xdhEdBWTMGuaPV1sp4nslESGcZI/iMLnyN95jLFHaT0+hWRlpacD8S2ds9FYkEZ3QffHDLnt2tIBCR5HUAXVHA3sNnHaisDK4V4d7s4Kh5kl4w+J3OXzgTw1jMKD6jCwNZzd+ZtM9FjaD+4MjKCnapOWGgIBCR1NZQWDRlDKOeVkZjTtKLVrN7qio4KkjnCc7GgCWM4nM6M5DV/INJTPnFSVx3XdM/AgWBiEhTHEgrY5+T9NgbHM0cy4i2h0ze/f0CBs9sepOgsUGQ0Yy6RESST17evv0vtd0fPLjBVkZa1f2o4LDauqcqK6vP7qYyUmdwZFolg0sLoZaB6paiFoGISKy0RPdUVhbMn9+sQQJ1DYmItDX1HTXVDOoaEhFpa+rrnoqhtFit2MweMLNiM1sV9VgXM3vJzNaEfw+O1fZFRKRxYhYEwEPApBqPXQvMd/ejgPnhsoiIxFHMgsDdFwKf1Xj4dODh8P7DwPRYbV9ERBonli2C2hzi7h8DhH971PVCM5tpZsvNbHlJSUmrFSgikmpaOwgazd1nufswdx/WvXv3eJcjIpK0WjsIPjWzngDh3+JW3r6IiNTQ2oePPgNcCNwa/n26MW9asWLFFjP7oJnb7AZsaeZ72yrtc2rQPie/A93fIxrzopidUGZmjwP5BDvyKfAz4CngCeBw4ENghrvXHFBu6TqWN+aEimSifU4N2ufk11r7G7MWgbufW8dT42K1TRERabqEHSwWEZHWkQpBMCveBcSB9jk1aJ+TX6vsb5uYdE5ERGInFVoEIiJSDwWBiEiKS+ogMLNJZvauma01s6Sb4M7M+phZgZm9bWZvmdnl4eNJP8urmaWb2etm9ly43M/MXgn3ebaZZcW7xpZkZp3NbI6ZvRN+33nJ/j2b2ZXhf9erzOxxM8tJtu+5KbM0W+Du8PdspZkd31J1JG0QmFk68FtgMnAMcK6ZHRPfqlpcBfCf7j4AGAH8MNzHVJjl9XLg7ajlXwF3hPv8OXBJXKqKnbuAF9z9aODrBPuetN+zmfUCLgOGufsgIB34N5Lve36Ixs/SPBk4KrzNBO5tqSKSNgiA4cBad1/n7nuAPxPMfpo03P1jd/9neH8HwY9DL5J8llcz6w18A7g/XDZgLDAnfElS7bOZdQROBv4A4O573H0rSf49E5zn1M7MMoBc4GOS7Htu4izNpwN/9MAyoHPVlD0HKpmDoBewMWp5U/hYUjKzvsBxwCs0YZbXNupO4GogEi53Bba6e0W4nGzf9ZFACfBg2B12v5kdRBJ/z+7+EfAbghkIPga2AStI7u+5Sl3fa8x+05I5CKyWx5LyWFkzaw88CVzh7tvjXU8smdlUoNjdV0Q/XMtLk+m7zgCOB+519+OAL0iibqDahP3ipwP9gMOAgwi6RmpKpu+5ITH77zyZg2AT0CdquTewOU61xIyZZRKEwJ/c/a/hw8k8y+soYJqZbSDo7htL0ELoHHYhQPJ915uATe7+Srg8hyAYkvl7Hg+sd/cSdy8H/gqMJLm/5yp1fa8x+01L5iB4DTgqPMogi2Cg6Zk419Siwr7xPwBvu/vtUU9VzfIKTZjltS1w9+vcvbe79yX4Tl92928BBcBZ4cuSbZ8/ATaaWf/woXHAapL4eyboEhphZrnhf+dV+5y033OUur7XZ4ALwqOHRgDbqrqQDpi7J+0NmAK8B7wPXB/vemKwfycRNA1XAm+EtykEfebzgTXh3y7xrjVG+58PPBfePxJ4FVgL/AXIjnd9LbyvxwLLw+/6KeDgZP+egZuBd4BVwCNAdrJ9z8DjBGMg5QT/4r+kru+VoGvot+Hv2ZsER1S1SB2aYkJEJMUlc9eQiIg0goJARCTFKQhERFKcgkBEJMUpCEREUpyCQEQkxSkIRERSnIJApJnMbLCZfWBm/y/etYgcCAWBSDO5+5sE01xcEO9aRA6EgkDkwBQDA+NdhMiBUBCIHJhbgWwzOyLehYg0l4JApJnMbBLBPPnPo1aBtGEKApFmMLMc4NfADwhmghwU34pEmk9BINI8PyW4fuwGFATSxikIRJoovEDMBIIro4GCQNo4XY9ARCTFqUUgIpLiFAQiIilOQSAikuIUBCIiKU5BICKS4hQEIiIpTkEgIpLi/j+1Tc1GWEUc7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_MSE_ridge(df_train.X,df_train.Y,df_test.X,df_test.Y,p=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lambda 100 seems to have lower error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### part(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [],
   "source": [
    "q4_train = sio.loadmat('PS2Q4.mat')\n",
    "q4_X = q4_train['x_tr']\n",
    "q4_Y = q4_train['y_tr']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=q4_X.reshape(1000,1)\n",
    "Y=q4_Y.reshape(1000,1)\n",
    "# Initialize beta\n",
    "beta = np.zeros((2,1))\n",
    "# Initialize mu\n",
    "mu = np.zeros((1000,1))\n",
    "\n",
    "# initialize gradient\n",
    "gradient = X.T@(mu-Y)\n",
    "\n",
    "# initialize epsilon\n",
    "epsilon = 0.0001\n",
    "# initialize eta\n",
    "learningrate = 10^(-5)\n",
    "\n",
    "while((gradient.T@gradient) > epsilon ):\n",
    "  # update beta\n",
    "  beta = beta - learningrate*(X.T@(mu-Y))\n",
    "  # update mu\n",
    "  for k in range(mu.shape[0]):\n",
    "    mu[k] = np.exp(beta[0] + beta[1] * (X[k]))\n",
    "  \n",
    "  # Calculte gradient\n",
    "  gradient = X.T@(mu-Y)\n",
    "  \n",
    "#   # Update \n",
    "#   maxIterations = maxIterations-1\n",
    "  \n",
    "  # Print \n",
    "  print(gradient.T@gradient)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}