\n",
"\n",
"__Given date:__ Monday September 27\n",
"\n",
"__Due date:__ Friday October 15"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Question 1 Least Absolute Shrinkage and Selection Operator \n",
"\n",
"__(13pts)__\n",
"\n",
"Learning a model through the OLS loss can be done very efficiently through either gradient descent or even through the Normal equations. The same is true for ridge regression. For the Lasso formulation however, the non differentiability of the absolute value at $0$ makes the learning more tricky.\n",
"\n",
"\n",
"\n",
"One approach, known as _ISTA_ (see Amir Beck and Marc Teboulle, _A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems_) consists in combining traditional gradient descent steps with a projection onto the $\\ell_1$ norm ball. Concretely, for the LASSO objective \n",
"\n",
"\\begin{align}\n",
"\\ell(\\boldsymbol \\beta) = \\|\\boldsymbol X\\boldsymbol \\beta - \\boldsymbol t\\|^2_2 + \\lambda \\|\\boldsymbol \\beta\\|_1\n",
"\\end{align}\n",
"\n",
"where $\\boldsymbol \\beta = (\\beta_1, \\beta_2,\\ldots, \\beta_D)$ (note that we don't include the bias) and the feature vectors $\\left\\{\\boldsymbol x_i\\right\\}_{i=1}^N$ (corresponding to the rows of the matrix $\\boldsymbol X$) as well as the targets $t_i$ are assumed to be centered, i.e.\n",
"\\begin{align}\n",
"\\boldsymbol x_{ij} \\leftarrow \\boldsymbol x_{ij}- \\frac{1}{N}\\sum_{i=1}^{N} x_{ij}\\\\\n",
"t_i \\leftarrow t_i - \\frac{1}{N}\\sum_{i=1}^N t_i\n",
"\\end{align}\n",
"\n",
"(Note that this is equivalent to taking $\\beta_0 = \\frac{1}{N}\\sum_{i=1}^N t_i$)\n",
"The ISTA update takes the form \n",
"\n",
"\\begin{align}\n",
"\\boldsymbol \\beta^{k+1} \\leftarrow \\mathcal{T}_{\\lambda \\eta} (\\boldsymbol \\beta^{k} - 2\\eta \\mathbf{X}^T(\\mathbf{X}\\mathbf{\\beta} - \\mathbf{t}))\n",
"\\end{align}\n",
"\n",
"where $\\mathcal{T}_{\\lambda \\eta}(\\mathbf{x})_i$ is the thresholding operator defined component-wise as\n",
"\n",
"\\begin{align}\n",
"\\mathcal{T}_{\\lambda t}(\\mathbf{\\beta})_i = (|\\beta_i| - \\lambda t)_+ \\text{sign}(\\beta_i)\n",
"\\end{align}\n",
"\n",
"In the equations above, $\\eta$ is an appropriate step size and $(x)_+ = \\max(x, 0)$ \n",
"\n",
"##### Question 1.1. (5pts)\n",
"\n",
"Complete the function 'ISTA' which must return a final estimate for the regression vector $\\mathbf{\\beta}$ given a feature matrix $\\mathbf{X}$, a target vector $\\mathbf{t}$ (the function should include the centering steps for $\\mathbf{x}_i$ and $t_i$) regularization weight $\\lambda$, and the choice for the learning rate $\\eta$. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def ISTA(beta_init, X, t, lbda, eta):\n",
" \n",
" '''The function takes as input an initial guess for beta, a set \n",
" of feature vectors stored in X and their corresponding \n",
" targets stored in t, a regularization weight lbda, \n",
" step size parameter eta and must return the \n",
" regression weights following from the minimization of \n",
" the LASSO objective'''\n",
" \n",
" beta_LASSO = np.zeros((np.shape(X)[0], 1))\n",
" \n",
" # add your code here\n",
" \n",
" \n",
" return beta_LASSO\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Question 1.2. (3pts)\n",
"\n",
"Apply your algorithm to the data (in red) given below for polynomial features up to degree 8-10 and for various values of $\\lambda$. Display the result on top of the true model (in blue)."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiQUlEQVR4nO3deXxU5b3H8c8vYZEAArJpCSEoO4IVAqioqBVlUamAdaFat3K1pe3VtgKCFatUsK27FKPXuqXSe2URVypu0AqytBrCDmExgLLJYiBAkuf+MYMOMSQHcjIzZ+b7fr18hZlzknmO4Jfjc873OeacQ0REgi8l1gMQERF/KNBFRBKEAl1EJEEo0EVEEoQCXUQkQdSI1Qc3adLEZWZmxurjRUQCafHixdudc03L2xazQM/MzGTRokWx+ngRkUAysw1H26YpFxGRBKFAFxFJEAp0EZEEoUAXEUkQCnQRkQShQBeR5JGTA5mZkJIS+pqTE+sR+arSQDez58xsq5nlVbJfDzMrMbOh/g1PRMQnOTkwfDhs2ADOhb4OH55Qoe7lDP15oF9FO5hZKjARmOXDmERE/DdmDOzbd+R7+/aF3o+iklJHaWn1LFteaaA75+YAOyvZ7RfAVGCrH4MSEfHdxo3H9n41WP3lXq6a/DF/X/R5tfz8Ks+hm1kL4Epgsod9h5vZIjNbtG3btqp+tIiIdxkZx/a+jw4Wl/L4e6sZ+OhH5K8uoN6NP66WOXw/qv+PAiOdcyVmVuGOzrlsIBsgKytLj0oSkegZPz40Zx457ZKWFnq/Gn32+S5GTs1lxRd7uXzVx9w7axJN9u0ObRw+PPR12DBfPsuPQM8CpoTDvAkwwMyKnXMzfPjZIiL+OByaY8aEplkyMkJh7lOYlrX/YAkPv7uS//nnOprVP4Fn5kym77w3jtzp8By+T2MwL88UNbNM4A3n3OmV7Pd8eL9XK/uZWVlZTotziUgi+njtdkZPW8KGHfu4rlcGo/p34MS02qG7a8oyg9JSzz/bzBY757LK21bpGbqZvQJcADQxswLgXqAmgHOu0nlzEZFksXv/ISa8vZxXFnxOZuM0XvnpWZx9WuPQxoyM0K2SZfk4h19poDvnrvX6w5xzN1ZpNCIi8Swn56hTNv9Y+gX3vJbHtr0H+K8+p3LHxe04oWbqt98bhTn8mK2HLiISKIeLSYcDOVxM2nbIGFezA28u2UKHk+vzzA1ZdE1v+N3vj8IcvgJdRMSLMsUkB0w99Szu/6wG+9O+5LeXtmf4+adSM7WCu8GHDau2i7CgQBcR8SaigFRwYlPuvnQEc07tTveCZUx8+HbaNKsXw8GFKNBFRLzIyKB0w0Ze7DaQh/r8BHOO+96dzPU78khp9ttYjw5QoIuIeLLmngcZ+dFmFn+vA+fnL+YPs54kvbgQsrNjPbRvKNBFRCpwsLiUpz9ayxP5DUlrWY+H//UCV855FcvIgPGPVOuc+LFSoIuIHEVuwS7uejVU27+s6ymMu6IzTepdHuthHZUCXUSkjP0HS3hk9iqenZtP0/q1eeaGLPp2ah7rYVVKgS4iEiGytn9tz5aMHtCRE0+oGetheaJAFxGhktp+QCjQRSTpVVrbDwgFuogkrW17DzDu9aW8mVtJbT8gFOgiknScc0z79yZ+/8Yy9h8s4dd923HbBadVXNsPAAW6iCSVgq/2cff0POas2kb3Vo2YOKQLbZrVj/WwfKFAF5GkUFrqeHHeeh6atRKAcZd34oazM0lJqfjRmUGiQBeRhLdm615GTl3C4g1f0addU8ZfeTrpjdJiPSzfKdBFJGEdKill8odreeL9NaTVTuXPV53B4G4tqOyB9kGlQBeRhBRZ2x/Y9RTGXd6ZpvVrx3pY1UqBLiIJpWxt/+nru3Np55NjPayoCPY9OiKSXHJyIDMTUlJCX3Nyjtg8b+0O+j82h+w5+VzdoyX/uKNP0oQ56AxdRILiKM/0BNgz5Ec8+NYKXlmwkVaN0/jbrb04p02TGA42NioNdDN7DrgM2OqcO72c7cOAkeGXXwO3O+c+83WUIiJlnukJwL59vPvUFMZuaM62vQcYfn6otl+nVvBq+37wcob+PPAk8OJRtq8D+jjnvjKz/kA20Muf4YmIhEU80xNge1oDxl38X7zR8Xw6pNUKfG3fD5UGunNujpllVrD944iX84F0H8YlInKkjAzYsAEHTO98Ib//wU/ZV7MOv8mdyX+N/0vga/t+8PvfwC3A20fbaGbDzWyRmS3atm2bzx8tIglt/HgKmmdw41XjuPOyX3Pqzk28NeUuRgw7X2Ee5ttFUTO7kFCgn3u0fZxz2YSmZMjKynJ+fbaIJLbSUsdLrXvz0M1P4Q4cZNzsp7l+xxJSJz4QV8/0jDVfAt3MugLPAv2dczv8+JkiIgBrtn7NqKm5LNrwFee1bcofruxCyz8PjvWw4lKVA93MMoBpwPXOuVVVH5KISKi2//RHa3n8vVBt/+EfncGVZyZubd8PlU48mdkrwDygvZkVmNktZnabmd0W3uV3QGNgkpl9amaLqnG8IhJklRSDDltSsJvLn/gnf/rHKvp2as67d/RhcLd0hXklvNzlcm0l228FbvVtRCJydDk5ofuxN24M3fUxfnxw5pArKAYdPob9B0t4dPYqnpmbT5N6yVXb94M5F5trk1lZWW7RIp3Mi3hWNhAB0tIgOzsYoZ6ZGQrxslq1gvXrmbd2B6On5bJ+xz6u6dGS0QM60qBOzagPM96Z2WLnXFa52xToIgFRSSDGvZQUKCdv9tSuy4N/m8crCzaScVIaEwZ3ScravlcVBbrWchEJijJNyUrfjzfhYlCkd9v0ZGz/X7Jt4UZ+el5r7uzbPmlr+35QoIsERTmB+M37QTB+/DdTRkfU9k8oIfuW3pzRsmGsRxh4CnSRoIgIxG+kpYXeD4Jhw3AOZjwzg/u+P5jCWmnc0Ww/t/9yMLVqqOnpBwW6SFAcvvAZ0LtcNu3az5hDbfnw7BvpltGQiUO60rZ5/VgPK6Eo0EWCZNiwwAT4YaWljpc/2cDEt1fggHsv78QNZ2eSmqJ7yv2mQBeRanNkbb9JqLZ/Ulqsh5WwFOgi4rtDJaVkz8nnsdmrqVMrlT9fdQaDu6m2X90U6CLiqyUFu7lrai7Lt+xhYJdTGHdFZ5rWrx3rYSUFBbqI+KLoUAmPzF7Fs3PX0bhuLSb/uDv9TldtP5oU6CJSZfPzdzB62hLWbS9UbT+GdPOnSDLxuNqhV3uKDnH39CVckz2fklLH327txYQhXRXmMaIzdJFk4WG1w2Mxe9mXjJ2Rx9a9RQw//1TuuLidavsxpsW5RJKFT4t7bf/6AONmLuWN3C10OLk+E4d0VW0/irQ4l4hUeXEv5xwzPt3E719fxtcHirmzbztu63OaavtxRIEukiyqsLjXpl37GTN9CR+u3KbafhxToIski+NY3Ku01JHzyQYmqLYfCAp0kWRxjIt7rd0Wqu0vXK/aflAo0EWSiYfFvb6p7b+3mjo1U/nTVWcwRLX9QFCgi8g3ytb2772iE83qnxDrYYlHlQa6mT0HXAZsdc6dXs52Ax4DBgD7gBudc//2e6AiUn3K1vafvr47l3ZWbT9ovJyhPw88Cbx4lO39gbbhf3oBfwl/FZEAiKztX53VkrsHdKRBmpqeQVRpoDvn5phZZgW7DAJedKGG0nwza2hmpzjntvg1SBHx396iQ0x4ewU5n2wk46Q0cm7tRe82TWI9LKkCP+bQWwCfR7wuCL+nQBeJU5G1/VvPbc2vL2mv2n4C8CPQy7v0Xe56AmY2HBgOkBGUJ5WLJJDtXx/gvteX8fpnm2nfvD6Tr+/O91XbTxh+dHYLgJYRr9OBzeXt6JzLds5lOeeymjZt6sNHi0SZz6sVRotzjhn/2UTfhz/inbwt3HFxO17/xbkK8wTjxxn6TGCEmU0hdDF0t+bPJSH5vFphtGwO1/Y/WLmNM8O1/Xaq7SekSldbNLNXgAuAJsCXwL1ATQDn3OTwbYtPAv0I3bZ4k3Ou0mUUtdqiBI5PqxVGS2mpI2fBRia+vYKSUsdvL23PT85RbT/oqrTaonPu2kq2O+Dnxzk2keCo4mqF0bR229eMnrqEBet3cm6bJjw4WLX9ZKCmqIhXVVitMFoia/sn1Ejhj0O7MrR7umr7SUKBLuLVcaxWGE15m3Zz16u5LNuyh/6nn8x9gzqrtp9kFOgiXh3jaoXRUnSohEdnr+aZufmcVLcWk3/cjX6nnxLTMUlsKNBFjoWH1Qqj6ZP8HYxSbV/CFOgiAaTavpRHgS4SMO8t/5Ix07+t7d95STvSauk/ZVGgiwTGjnBtf6Zq+3IUCnSROOecY+Znmxk3cylfHyjmvy9uy88uaEOtGn6s3CGJRH8iROLY5l37ufn5hfxqyqdkNqnLm788j/++uF3swjyga9kkC52hi8ShsrX9313WKfa1/YCuZZNMKl3LpbpoLReR8sVtbT9ga9kkqorWctGUi0g0VTBlcaiklEkfrqH/Y3NZ8cUeHhralZdu6RkfYQ6BWssmWWnKRSRaKpiyyLvgsviv7QdgLZtkp0AXiZYxY45cBwYoOljMY3//hOyljeK/th/na9mIAl0kespMTSxI78yofr8gv3E6V3dLj//afpyuZSPf0kVRkWgJX1TcW6sOE/vcyMvdBtJy1xdM+M//0vuTWbEenQRElR5wISI+GT+e9//wF8b0uZUv653ErQumc+fiaaRNeiLWI5MEoUAXiYIdXx/g96mdeO3y0bTftYlJOXdxZs0imPSEpizENwp0kWpUtrZ/x8XtuP2C/tSaPDzWQ5MEpEAXqSabd+1n7Iw83l+xle+3bMhDQ7vSrnn9WA9LEpgCXcRnZWv7Ywd25KberWNb25ekoEAX8VH+tq8ZFY+1fUkKngLdzPoBjwGpwLPOuQlltjcAXgYywj/zT865v/o8VpG4daiklGfm5vPo7NWcUCOFh4Z25aru6ZjprFyip9JAN7NU4CmgL1AALDSzmc65ZRG7/RxY5py73MyaAivNLMc5d7BaRi0SR/I27Wbk1FyWbg7X9q/oTLMT46y2L0nByxl6T2CNcy4fwMymAIOAyEB3QH0LnY7UA3YCxT6PVSSuFB0q4bH3VpM9Jz/+a/uSFLwEegvg84jXBUCvMvs8CcwENgP1gaudc6Vlf5CZDQeGA2RoQR8JsAXrdjJqai752wu5Oqtl/Nf2JSl4CfTyJgHLrhdwKfApcBFwGvCumc11zu054pucywayIVT9P+bRisTY3qJDPPTOSl6av4GWJ9Xh5Vt6cW7bJrEelgjgLdALgJYRr9MJnYlHugmY4EILw6wxs3VAB2CBL6MUiQPvr/iSMdPz+GJPETf3bs1vLm1HWi3dKCbxw8ufxoVAWzNrDWwCrgGuK7PPRuAHwFwzaw60B/L9HKhIrOwsPMjvX1/KjE830655PZ4adg7dMhrFelgi31FpoDvnis1sBDCL0G2LzznnlprZbeHtk4H7gefNbAmhKZqRzrnt1ThukWp3uLZ/3+vL2Ft0iP++uC0/u6BN7B7QLFIJT/+/6Jx7C3irzHuTI369GbjE36GJxM6W3fsZOz2P91Zs5YyWDXloSFfan6zavsQ3nWpIcqngmZ4Qqu2/PH8DfR+ew8drdzB2YEem3X6OwlwCQVd0JHlU8ExPhg0L1fanLWHBup30btOYB6/sSkZj1fYlOBTokjzKeaYn+/ZRPPYenmlxFo/MXqXavgSaAl2SR5lnegIsbdaakX1+Rd47K+jX+WR+P0i1fQkuBbokj4yM0DQLUJRak8d7X8vTvYbQ6EAhfxnWjf5dVNuXYNNFUUke48dDWhoL0jsz4KYnmHT2jxi8Yg6zuzmFuSQEnaFL0tg75Ec8tKkOL+2sTfquL3npgyc475fXw7CyPTmRYFKgS1L4YMVWxkxfwpY9tbm5d2t+fcml1K19c6yHJeIrBboktMjafttm9Zh6u2r7krgU6JKQnHO8nruFcTOXsrfoEL/6QVt+duFp1K6RGuuhiVQbBbokHNX2JVkp0CVhlJY6Xlm4kQffWkFJqWPswI7c1Ls1qSkqCElyUKBLQlBtX0SBLgFXXFLKM3PX8ejsVdRWbV+SnIpFEl2VrHZ4LJZu3s0PJ/2Lie+s4IL2TZl9Zx9+lNVSYS5JS2foEj2VrHboVdGhEh5/bzVPz8mnUVot1fZFwiz0GNDoy8rKcosWLYrJZ0uMZGZ+s5bKEVq1gvXrPf2IBet2MmpqLvnbC7mqezpjBnakYVotX4cpEs/MbLFzLqu8bZpykWNTlSmTclY7rPD9CHuLDnHPjDx+9PQ8Dq7fwEt/v4c//nYQDaf/n/fPF0lwmnIR76o6ZRKx2uF33q/AN7X93fu56dO3+M37f6XuoaLQxuOYshFJVDpDF++O8oAIxozx9v3h1Q6PkJYWer8cOwsPcsffP+Wm5xdSt3YNXn33z9w76y/fhvmxfr5IgtMZunhXhSkT4Nuz6DFjQt+TkREK8zJn10et7f/mo6p9vkiC8xToZtYPeAxIBZ51zk0oZ58LgEeBmsB251wf30Yp8eE4p0yOMGxYhdMjR9T20xswcWgvOpx8on+fL5LAKp1yMbNU4CmgP9AJuNbMOpXZpyEwCbjCOdcZuMr/oUrMHeOUybEoLXW8PH8DfR+ew7/WbmfswI5M+1nvb8O8mj9fJBF4OUPvCaxxzuUDmNkUYBCwLGKf64BpzrmNAM65rX4PVOKAxymTY7VueyEjp+ayYN1OzjmtMRMGH6W2X02fL5IoKr0P3cyGAv2cc7eGX18P9HLOjYjY51FCUy2dgfrAY865F8v5WcOB4QAZGRndN5T3v8+SNIpLSnn2n+t45N1V1KqRwtiBHdX0FKlERfehezlDL++/rrJ/C9QAugM/AOoA88xsvnNu1RHf5Fw2kA2hYpGHz5YEtWzzHu6a+hl5m/ZwSafm3P/D02l+4gmxHpZIoHkJ9AKgZcTrdGBzOftsd84VAoVmNgc4A1iFSISiQyU88f5qJn+UT6O0mkwa1o3+p5+ss3IRH3gJ9IVAWzNrDWwCriE0Zx7pNeBJM6sB1AJ6AY/4OVAJvoXrdzJyai752woZ2j2dsarti/iq0kB3zhWb2QhgFqHbFp9zzi01s9vC2yc755ab2TtALlBK6NbGvOocuATH3qJD/HHWSl6ct4H0RnV48eaenN+uaayHJZJwtDiXVKtvavt7irjxnEx+c0l76tZWn03keFX1oqjIMdtZeJD731jG9P9som2zerx62zl0b9Uo1sMSSWgKdPGVc443wrX93fsP8csftOXnF55G7RqpsR6aSMJToItvtuzezz0z8pi9PFTbz/lpryObniJSrRToUmWlpY4pCz/nwbeWc6i0lDEDOnLzua1JTdGtiCLRpECXKlm3vZBRU3P5JFzbf3BwF1o1rhvrYYkkJQW6HJfiklL+55/reDhc258wuAtX91BtXySWFOhyzJZt3sPIqbks2bSbvp2a84Bq+yJxQYEunhUdKuHJ99cw+aO1NFRtXyTu6BF0yeY4H/K8aP1OBjw+lyc/WMOg77dg9p19GNDlFIW5SBzRGXoyOY6HPH99oJg/vrOCF+dv4HsNVNsXiWeq/ieTzMzyH+HWqhWsX/+dtz9YuZUx00K1/Z+cnclvL1VtXyTWVP2XEI8Pef4qXNuf9p9NnNa0Lq/edjbdW50UhQGKSFUo0JNJJQ9Z/k5t/6I2/PyiNqrtiwSEAj2ZjB9/5Bw6fPOQ5S92FzF2Rh6zl39J1/QGvHxrLzqeotq+SJAo0JNJOQ9ZLn1gPFPanMeDD3/0TW3/pt6Z1EjVDVAiQaOLokls3fZCRk/LZX7+Ts4+tTEThqi2LxLvdFFUjqDavkhiUqAnmcja/iWdmnO/avsiCUOBniRU2xdJfAr0JLB4w07uejWXtdsKGdItnbEDO9Kobq1YD0tEfKZAT2CFB4r546yVvDBvPd9rUIcXbu5JH9X2RRKWp3vTzKyfma00szVmNqqC/XqYWYmZDfVviHI8Plq1jUsemcML89bzk7Mz+ccd5yvMRRJcpWfoZpYKPAX0BQqAhWY20zm3rJz9JgKzqmOg4s1XhQe5/81lTPu3avsiycbLlEtPYI1zLh/AzKYAg4BlZfb7BTAV6OHrCMUT5xxvLgnV9nftU21fJBl5CfQWwOcRrwuAXpE7mFkL4ErgIioIdDMbDgwHyAivHyJVV7a2/9Itqu2LJCMvgV7efW1l66WPAiOdcyUV3QbnnMsGsiHUFPU4RjkK5xxTFn7OH95crtq+iHgK9AKgZcTrdGBzmX2ygCnhMG8CDDCzYufcDD8GKd+1fnsho6ctYV7+Ds469SQmDO5KZhPV9kWSmZdAXwi0NbPWwCbgGuC6yB2cc60P/9rMngfeUJhXj+KSUp77V6i2XzMlhQcHd+Ea1fZFBA+B7pwrNrMRhO5eSQWec84tNbPbwtsnV/MYJWz5llBtP7dgN307Nef+QadzcgPV9kUkxFOxyDn3FvBWmffKDXLn3I1VH5ZEOlAcqu3/5cNQbf+p67oxoItq+yJyJDVF41xkbX9wtxbcM7CTavsiUi7dDhGnCg8UM27mUoZOnkfRoVJeuLknDx9aSqPO7SAlJfTA55ycWA9TROKIztDj0IcrtzJmeh6bd+/nhrNa8dt+Haj36t+PfHzchg2h1/Dtk4hEJKnpiUVxpGxt/6GhXb+t7Wdmlv+A51atYP36aA5TRGJITyyKc2Vr+yMubMOIi9pwQs2I2v7GjeV/89HeF5Gko0CPsS92F3HPa3m8u6yS2n5GRvln6FpCQUTCgnVRNCcnNPWQABcFS0sdf/tkI30f/oi5q7dx94AOTLv9nKOvwTJ+PKSlHfleWlrofRERgnSGnpOTMBcF128vZNS0XObn7/Re2z98jGPGhKZZMjJCYR6wYxeR6hOci6IJcFHwcG3/z/9YRa3UFEYP6Mi1PVXbFxHvEuOiaMAvCkbW9i/u2JwHfqjavoj4KziBHtCLgpG1/QZ1avLkdWcysMspOisXEd8FJ9DHjz9yDh3i/qLg4g07GTl1CWu2fs3gM1twz2Wq7YtI9QnOXS7DhkF2dmjO3Cz0NTs7uhcFPd5lE1nb33egmOdv6sHDV39fYS4i1So4Z+gQCu9Y3dXh8S6bOau2MXrakiNr+7WD9a9ZRIIpOHe5xFold9ns2neQ+99YztR/F3Bq07o8NKQrWZknRX2YIpLYEuMul1g7yt00buNG3srdwr0z845e2xcRiQIFulfl3GXzZb2TGDvo17z7t3/TpUUDXry5F52+d5Smp4hINVOgexVxl40D/t71EsZfdAsHT0jj7v4duLl3a2qkBucas4gkHgW6V+ELnxsefIRRXQYzr9UZ9KpbzMTbL6y8ti8iEgXJdUpZhcW9iktKeabl2Vx65f3ktevOH67switjrlCYi0jcSJ4z9Cos7nVkbb8ZD/ywi2r7IhJ3kue2xeNY3OtAcQlPvb+GSeHa/n2DOqu2LyIxVeXbFs2sH/AYkAo865ybUGb7MGBk+OXXwO3Ouc+Of8jV4BgX94qs7V95Zgt+p9q+iMS5SgPdzFKBp4C+QAGw0MxmOueWRey2DujjnPvKzPoD2UCv6hjwcfO4uFfhgWL+OGslL8xbz/ca1OH5m3pwQftmURqkiMjx83KG3hNY45zLBzCzKcAg4JtAd859HLH/fCDdz0H6wsPiXodr+5t27eeGs1txl2r7IhIgXtKqBfB5xOsCKj77vgV4u7wNZjYcGA6QEe1lbyt44k/Z2v7/3XY2PVTbF5GA8RLo5V0BLPdKqpldSCjQzy1vu3Mum9B0DFlZWdG/GltmcS/n3BG1/Z9feBq/uKitavsiEkheAr0AaBnxOh3YXHYnM+sKPAv0d87t8Gd41efLPUXcMyOPfyz7ktNbnMgLN/ek8/caxHpYIiLHzUugLwTamllrYBNwDXBd5A5mlgFMA653zq3yfZQ+cs7xv4s+54E3l3OwuJSR/Trw0/NU2xeR4Ks00J1zxWY2AphF6LbF55xzS83stvD2ycDvgMbApPA92sVHu08yljbu2Mfo6bn8a80OerU+iQlDutJaTU8RSRBJUSwqKXU89891/PndldRMSWHUgA5c2yODlBQVhEQkWJJ6PfQVX+xh5Ku5fBau7d//w9M5pUGdWA9LRMR3CRvoB4pLeOqDtUz6YA0N6tTk8WvP5PKuqu2LSOJKyEBfvOErRk3NZXW4tn/PZZ04SbV9EUlwCXVrR+GBYu57fSlDJ39M4YFi/npjDx65+vvxE+ZVWL5XRKQyCXOGPnd1qLZf8FWc1varsHyviIgXgb/LZde+gzzw5nJeXRyq7U8c0jU+a/vHsXyviEhZCXmXi3OOt/O+4HevLeWrfQcZcWEbRlzUJn5r+8e4fK+IyLEKZKBv3VPEPa/lMWvp4dp+j/iv7XtcvldE5HgFLtA/WLmVX77yHw4WlzK6fwduOTcgtX0Py/eKiFRF4AK9deO6dMtoxLgrOgertl/B8r0iIn4I/EVREZFkUtFF0QDMVYiIiBcKdBGRBKFAFxFJEAp0EZEEoUAXEUkQCnQRkQShQBcRSRAKdBGRBBGzYpGZbQPKWdzEkybAdh+HEwQ65uSgY04OVTnmVs65puVtiFmgV4WZLTpaUypR6ZiTg445OVTXMWvKRUQkQSjQRUQSRFADPTvWA4gBHXNy0DEnh2o55kDOoYuIyHcF9QxdRETKUKCLiCSIuA50M+tnZivNbI2ZjSpnu5nZ4+HtuWbWLRbj9JOHYx4WPtZcM/vYzM6IxTj9VNkxR+zXw8xKzGxoNMdXHbwcs5ldYGafmtlSM/so2mP0m4c/2w3M7HUz+yx8zDfFYpx+MbPnzGyrmeUdZbv/+eWci8t/gFRgLXAqUAv4DOhUZp8BwNuAAWcBn8R63FE45nOARuFf90+GY47Y733gLWBorMcdhd/nhsAyICP8ulmsxx2FY74bmBj+dVNgJ1Ar1mOvwjGfD3QD8o6y3ff8iucz9J7AGudcvnPuIDAFGFRmn0HAiy5kPtDQzE6J9kB9VOkxO+c+ds59FX45H0iP8hj95uX3GeAXwFRgazQHV028HPN1wDTn3EYA51zQj9vLMTugvpkZUI9QoBdHd5j+cc7NIXQMR+N7fsVzoLcAPo94XRB+71j3CZJjPZ5bCP0NH2SVHrOZtQCuBCZHcVzVycvvczugkZl9aGaLzeyGqI2ueng55ieBjsBmYAnwK+dcaXSGFxO+51eNKg2nelk575W9x9LLPkHi+XjM7EJCgX5utY6o+nk55keBkc65ktDJW+B5OeYaQHfgB0AdYJ6ZzXfOraruwVUTL8d8KfApcBFwGvCumc11zu2p5rHFiu/5Fc+BXgC0jHidTuhv7mPdJ0g8HY+ZdQWeBfo753ZEaWzVxcsxZwFTwmHeBBhgZsXOuRlRGaH/vP7Z3u6cKwQKzWwOcAYQ1ED3csw3ARNcaIJ5jZmtAzoAC6IzxKjzPb/iecplIdDWzFqbWS3gGmBmmX1mAjeErxafBex2zm2J9kB9VOkxm1kGMA24PsBna5EqPWbnXGvnXKZzLhN4FfhZgMMcvP3Zfg04z8xqmFka0AtYHuVx+snLMW8k9H8kmFlzoD2QH9VRRpfv+RW3Z+jOuWIzGwHMInSF/Dnn3FIzuy28fTKhOx4GAGuAfYT+hg8sj8f8O6AxMCl8xlrsArxSncdjTihejtk5t9zM3gFygVLgWedcube/BYHH3+f7gefNbAmh6YiRzrnALqtrZq8AFwBNzKwAuBeoCdWXX6r+i4gkiHiechERkWOgQBcRSRAKdBGRBKFAFxFJEAp0EZEEoUAXEUkQCnQRkQTx/2M32Zz70Ne9AAAAAElFTkSuQmCC\n",
"text/plain": [
"