RegularizedRegression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# math \n",
    "import numpy as np\n",
    "from scipy.linalg import norm\n",
    "\n",
    "# random\n",
    "from numpy.random import rand, randn\n",
    "np.random.seed(6)\n",
    "\n",
    "# plotting\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# robust regression \n",
    "from sklearn.linear_model import LinearRegression, Lasso, Ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_data(n, w):\n",
    "    X = randn(n,len(w))\n",
    "    y = X@w\n",
    "    return X, y\n",
    "\n",
    "def generate_noisy_data(n, w):\n",
    "    X = randn(n,len(w))\n",
    "    y = X@w + .1*randn(n)\n",
    "    return X, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OLS is invariant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 10\n",
    "n = 20\n",
    "w = randn(d)\n",
    "X, y = generate_data(n, w);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "LinearRegression?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.71420156, -1.63516395, -3.8608605 , -0.57894597,  6.59638624,\n",
       "       -4.79152845, -0.5736281 ,  5.18756905, -7.26414655,  5.32047682,\n",
       "        3.56641161, -0.245102  ,  1.61110249,  2.61211864, -3.38682977,\n",
       "        1.19532532,  4.86224655,  6.26788358,  4.48640507, -3.14514452])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin = LinearRegression().fit(X,y)\n",
    "yhat = lin.predict(X)\n",
    "yhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -3.57100781,  -8.17581975, -19.30430252,  -2.89472986,\n",
       "        32.98193119, -23.95764226,  -2.8681405 ,  25.93784523,\n",
       "       -36.32073277,  26.6023841 ,  17.83205803,  -1.22551002,\n",
       "         8.05551247,  13.06059321, -16.93414883,   5.9766266 ,\n",
       "        24.31123275,  31.33941792,  22.43202537, -15.72572261])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yscale = 5*y\n",
    "Xscale = 3*X\n",
    "lin = LinearRegression().fit(Xscale,yscale)\n",
    "yhatscale = lin.predict(Xscale)\n",
    "yhatscale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.4221204859038887e-14"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# predictions match up to scaling\n",
    "norm(yhat - yhatscale/5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ridge regression is not scaling invariant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ridge?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.73252921, -1.61138302, -3.871725  , -0.59222297,  6.55771702,\n",
       "       -4.72195274, -0.56508244,  5.18633252, -7.21591302,  5.28679447,\n",
       "        3.54741707, -0.21884335,  1.60849856,  2.58391179, -3.40339278,\n",
       "        1.19486892,  4.83695027,  6.26115137,  4.48779252, -3.108016  ])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha = .1 # regularization parameter\n",
    "\n",
    "lin = Ridge(alpha=alpha).fit(X,y)\n",
    "yhat = lin.predict(X)\n",
    "yhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -3.58216226,  -8.16285939, -19.31217038,  -2.90274427,\n",
       "        32.95995686, -23.91745638,  -2.8628421 ,  25.93815093,\n",
       "       -36.29302326,  26.58314257,  17.82078142,  -1.20993494,\n",
       "         8.0550225 ,  13.04380232, -16.94425996,   5.97684365,\n",
       "        24.29672683,  31.33609432,  22.43311342, -15.70431196])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yscale = 5*y\n",
    "Xscale = 3*X\n",
    "lin = Ridge(alpha=alpha).fit(Xscale,yscale)\n",
    "yhatscale = lin.predict(Xscale)\n",
    "yhatscale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10923658080907492"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# predictions do not match up to scaling\n",
    "norm(yhat - yhatscale/5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.1041804156888074e-15"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# standardize\n",
    "\n",
    "lin = Ridge(alpha=alpha, normalize=True).fit(X,y)\n",
    "yhat = lin.predict(X)\n",
    "yhat\n",
    "\n",
    "yscale = 5*y\n",
    "Xscale = 3*X\n",
    "lin = Ridge(alpha=alpha, normalize=True).fit(Xscale,yscale)\n",
    "yhatscale = lin.predict(Xscale)\n",
    "yhatscale\n",
    "\n",
    "# predictions match up to scaling\n",
    "norm(yhat - yhatscale/5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# let's compare different kinds of regularized regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate data\n",
    "\n",
    "d = 30\n",
    "w_randn = randn(d)\n",
    "w_sparse = w_randn*(abs(w_randn) > 1)\n",
    "w_pos = w_randn*(w_randn > 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.        ,  0.35483766,  0.2628088 , -0.        , -0.        ,\n",
       "       -0.        , -0.        ,  1.13694227,  1.18705819,  1.69089583,\n",
       "        0.49028137, -0.        , -0.        , -0.        , -0.        ,\n",
       "        1.37681404, -0.        , -0.        ,  0.52907247,  0.73336413,\n",
       "        0.96671585, -0.        , -0.        , -0.        ,  0.74307591,\n",
       "        0.19329784, -0.        ,  1.89687218, -0.        , -0.        ])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w_pos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.        ,  0.        ,  0.        , -0.        , -1.37105209,\n",
       "       -1.06862477, -0.        ,  1.13694227,  1.18705819,  1.69089583,\n",
       "        0.        , -1.44553018, -0.        , -1.15404442, -0.        ,\n",
       "        1.37681404, -0.        , -0.        ,  0.        ,  0.        ,\n",
       "        0.        , -0.        , -0.        , -0.        ,  0.        ,\n",
       "        0.        , -0.        ,  1.89687218, -0.        , -0.        ])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w_sparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find best model for each type of data\n",
    "w = w_randn\n",
    "\n",
    "alpha_ridge=.1\n",
    "alpha_lasso=.1\n",
    "\n",
    "X,y = generate_data(30, w)\n",
    "w_ridge = Ridge(alpha = alpha_ridge).fit(X,y).coef_\n",
    "w_lasso = Lasso(alpha = alpha_lasso).fit(X,y).coef_\n",
    "w_nonneg = LinearRegression(positive=True).fit(X,y).coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([3., 1., 2., 2., 4., 0., 1., 1., 3., 0., 2., 3., 0., 2., 0., 0., 2.,\n",
       "        1., 1., 2.]),\n",
       " array([-1.37606868, -1.23186246, -1.08765624, -0.94345001, -0.79924379,\n",
       "        -0.65503756, -0.51083134, -0.36662512, -0.22241889, -0.07821267,\n",
       "         0.06599355,  0.21019978,  0.354406  ,  0.49861223,  0.64281845,\n",
       "         0.78702467,  0.9312309 ,  1.07543712,  1.21964334,  1.36384957,\n",
       "         1.50805579]),\n",
       " <BarContainer object of 20 artists>)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQEUlEQVR4nO3df6jdd33H8efLLKLMQudyXbP8aJSVMStUyyW2K4wgbvQXywYdVJiVMggtFRSEkTmo+J/uDxltXEOYxZaJIui60Ca4brO0/SPVJEvaxtgZpaMhwcaKqaFFF/feH/erXE/Pued7cs/99dnzAYf7/fE53/P+3E/yyjff8/2RqkKStPa9aaULkCRNh4EuSY0w0CWpEQa6JDXCQJekRvzGSn3whg0batu2bSv18ZK0Jh05cuRHVTUzbN2KBfq2bds4fPjwSn28JK1JSf571DoPuUhSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RG9A70JOuS/GeSR4esS5L7kpxK8mySa6dbpiRpnEn20D8GnByx7ibgqu61C3hgkXVJkibUK9CTbAZuAf5xRJOdwMM15xBweZKNU6pRktRD3ytF/x74a+CyEes3AS/Nmz/dLTs7v1GSXcztwbN169ZJ6vx/b9vuxxb1/hc/c8uUKpG0Wo3dQ09yK/ByVR1ZqNmQZW94FFJV7auq2aqanZkZeisCSdIl6nPI5QbgT5O8CHwF+ECSfxpocxrYMm9+M3BmKhVKknoZG+hV9TdVtbmqtgG3A/9RVX850Gw/cEd3tst1wPmqOju4LUnS0rnkuy0muQugqvYCB4CbgVPAa8CdU6lOktTbRIFeVU8AT3TTe+ctL+CeaRYmSZqMV4pKUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhrR5yHRb0nyrSTHk5xI8ukhbXYkOZ/kWPe6d2nKlSSN0ueJRT8DPlBVF5KsB55OcrCqDg20e6qqbp1+iZKkPsYGevd4uQvd7PruVUtZlCRpcr2OoSdZl+QY8DLweFU9M6TZ9d1hmYNJrp5mkZKk8XoFelX9oqreC2wGtid5z0CTo8CVVXUNcD/wyLDtJNmV5HCSw+fOnbv0qiVJbzDRWS5V9RPgCeDGgeWvVtWFbvoAsD7JhiHv31dVs1U1OzMzc8lFS5LeqM9ZLjNJLu+m3wp8EPjuQJsrkqSb3t5t95WpVytJGqnPWS4bgYeSrGMuqL9aVY8muQugqvYCtwF3J7kIvA7c3n2ZKklaJn3OcnkWeN+Q5XvnTe8B9ky3NEnSJLxSVJIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhrR55mib0nyrSTHk5xI8ukhbZLkviSnkjyb5NqlKVeSNEqfZ4r+DPhAVV1Ish54OsnBqjo0r81NwFXd6/3AA91PSdIyGbuHXnMudLPru9fgA6B3Ag93bQ8BlyfZON1SJUkL6bOHTpJ1wBHg94DPV9UzA002AS/Nmz/dLTs7sJ1dwC6ArVu3XmLJsG33Y5f8XoAXP3PLot6v5bGYcV7JMV6rdWvt6/WlaFX9oqreC2wGtid5z0CTDHvbkO3sq6rZqpqdmZmZuFhJ0mgTneVSVT8BngBuHFh1Gtgyb34zcGYxhUmSJtPnLJeZJJd3028FPgh8d6DZfuCO7myX64DzVXUWSdKy6XMMfSPwUHcc/U3AV6vq0SR3AVTVXuAAcDNwCngNuHOJ6pUkjTA20KvqWeB9Q5bvnTddwD3TLU2SNAmvFJWkRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RG9Hmm6JYk30xyMsmJJB8b0mZHkvNJjnWve5emXEnSKH2eKXoR+ERVHU1yGXAkyeNV9Z2Bdk9V1a3TL1GS1MfYPfSqOltVR7vpnwIngU1LXZgkaTITHUNPso25B0Y/M2T19UmOJzmY5OoR79+V5HCSw+fOnZu8WknSSL0DPcnbgK8BH6+qVwdWHwWurKprgPuBR4Zto6r2VdVsVc3OzMxcYsmSpGF6BXqS9cyF+Zeq6uuD66vq1aq60E0fANYn2TDVSiVJC+pzlkuALwAnq+pzI9pc0bUjyfZuu69Ms1BJ0sL6nOVyA/Bh4Lkkx7plnwS2AlTVXuA24O4kF4HXgdurqqZfriRplLGBXlVPAxnTZg+wZ1pFSZIm55WiktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1Ig+zxTdkuSbSU4mOZHkY0PaJMl9SU4leTbJtUtTriRplD7PFL0IfKKqjia5DDiS5PGq+s68NjcBV3Wv9wMPdD8lSctk7B56VZ2tqqPd9E+Bk8CmgWY7gYdrziHg8iQbp16tJGmkPnvov5JkG/A+4JmBVZuAl+bNn+6WnR14/y5gF8DWrVsnLHV12Lb7sZUuQctgLY7zYmp+8TO3TLGSyVj39PT+UjTJ24CvAR+vqlcHVw95S71hQdW+qpqtqtmZmZnJKpUkLahXoCdZz1yYf6mqvj6kyWlgy7z5zcCZxZcnSeqrz1kuAb4AnKyqz41oth+4ozvb5TrgfFWdHdFWkrQE+hxDvwH4MPBckmPdsk8CWwGqai9wALgZOAW8Btw59UolSQsaG+hV9TTDj5HPb1PAPdMqSpI0Oa8UlaRGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEb0eabog0leTvL8iPU7kpxPcqx73Tv9MiVJ4/R5pugXgT3Awwu0eaqqbp1KRZKkSzJ2D72qngR+vAy1SJIWYVrH0K9PcjzJwSRXj2qUZFeSw0kOnzt3bkofLUmC6QT6UeDKqroGuB94ZFTDqtpXVbNVNTszMzOFj5Yk/dKiA72qXq2qC930AWB9kg2LrkySNJFFB3qSK5Kkm97ebfOVxW5XkjSZsWe5JPkysAPYkOQ08ClgPUBV7QVuA+5OchF4Hbi9qmrJKpYkDTU20KvqQ2PW72HutEZJ0grySlFJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqxNhAT/JgkpeTPD9ifZLcl+RUkmeTXDv9MiVJ4/TZQ/8icOMC628Crupeu4AHFl+WJGlSYwO9qp4EfrxAk53AwzXnEHB5ko3TKlCS1M/Yh0T3sAl4ad786W7Z2cGGSXYxtxfP1q1bp/DRl2bb7sdW7LPXIn9fWq38s/nrpvGlaIYsq2ENq2pfVc1W1ezMzMwUPlqS9EvTCPTTwJZ585uBM1PYriRpAtMI9P3AHd3ZLtcB56vqDYdbJElLa+wx9CRfBnYAG5KcBj4FrAeoqr3AAeBm4BTwGnDnUhUrSRptbKBX1YfGrC/gnqlVJEm6JF4pKkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY3oFehJbkzyQpJTSXYPWb8jyfkkx7rXvdMvVZK0kD7PFF0HfB74Y+A08O0k+6vqOwNNn6qqW5egRklSD3320LcDp6rqB1X1c+ArwM6lLUuSNKk+gb4JeGne/Olu2aDrkxxPcjDJ1cM2lGRXksNJDp87d+4SypUkjdIn0DNkWQ3MHwWurKprgPuBR4ZtqKr2VdVsVc3OzMxMVKgkaWF9Av00sGXe/GbgzPwGVfVqVV3opg8A65NsmFqVkqSx+gT6t4GrkrwzyZuB24H98xskuSJJuunt3XZfmXaxkqTRxp7lUlUXk3wU+AawDniwqk4kuatbvxe4Dbg7yUXgdeD2qho8LCNJWkJjAx1+dRjlwMCyvfOm9wB7pluaJGkSXikqSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjegV6EluTPJCklNJdg9ZnyT3deufTXLt9EuVJC1kbKAnWQd8HrgJeDfwoSTvHmh2E3BV99oFPDDlOiVJY/TZQ98OnKqqH1TVz4GvADsH2uwEHq45h4DLk2yccq2SpAX0eUj0JuClefOngff3aLMJODu/UZJdzO3BA1xI8gKwAfjRBDWvBauuT/nsojex6vq0kJ79XXV9WuQ4XXJ/pvDnY6msujGagg357KL6dOWoFX0CPUOW1SW0oar2Aft+7Y3J4aqa7VHHmmGf1obW+tRaf8A+TarPIZfTwJZ585uBM5fQRpK0hPoE+reBq5K8M8mbgduB/QNt9gN3dGe7XAecr6qzgxuSJC2dsYdcqupiko8C3wDWAQ9W1Ykkd3Xr9wIHgJuBU8BrwJ0T1LBvfJM1xz6tDa31qbX+gH2aSKrecKhbkrQGeaWoJDXCQJekRix7oCf5iyQnkvxvkpGn7iR5MclzSY4lObycNU5qgj4teAuF1STJ25M8nuR73c/fGtFuVY9Ti7et6NGnHUnOd2NyLMm9K1FnX0keTPJykudHrF+LYzSuT0szRlW1rC/gD4DfB54AZhdo9yKwYbnrW6o+MfeF8veBdwFvBo4D717p2hfo098Bu7vp3cBn19o49fmdM/dl/kHmrqW4DnhmpeueQp92AI+udK0T9OmPgGuB50esX1Nj1LNPSzJGy76HXlUnq+qF5f7cpdSzT31uobCa7AQe6qYfAv5s5Uq5ZC3etmKt/Tkaq6qeBH68QJO1NkZ9+rQkVvMx9AL+NcmR7pYBa92o2yOsVr9T3bUE3c93jGi3msepz+98rY1L33qvT3I8ycEkVy9PaUtmrY1RX1Mfoz6X/k8syb8BVwxZ9bdV9S89N3NDVZ1J8g7g8STf7f7VWxFT6FOv2yMsp4X6NMFmVtU4DZjabStWkT71HgWurKoLSW4GHmHuTqhr1Voboz6WZIyWJNCr6oNT2MaZ7ufLSf6Zuf9qrlhQTKFPq+72CAv1KckPk2ysqrPdf29fHrGNVTVOA1q8bcXYeqvq1XnTB5L8Q5INVbVWb3K11sZorKUao1V5yCXJbya57JfTwJ8AQ78tXkP63EJhNdkPfKSb/gjwhv+FrIFxavG2FWP7lOSKJOmmtzP39/yVZa90etbaGI21ZGO0At/+/jlz/+L+DPgh8I1u+e8CB7rpdzH37f1x4ARzhzVW/JvrxfSpm78Z+C/mzlJY7X36beDfge91P9++Fsdp2O8cuAu4q5sOcw9w+T7wHAucebVaXj369NFuPI4Dh4A/XOmax/Tny8zdavt/ur9Hf9XAGI3r05KMkZf+S1IjVuUhF0nS5Ax0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1Ij/A7yP/CZA2hh1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(w_ridge, label=\"ridge coefficients\", bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcceea10940>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYRklEQVR4nO3deXRV5b3/8fcXEpsyiiR6RSpoV52ahBijjEW4QbCo8EOwVsGB/pCyqsVrf1CgdgHLYVFvWUKhVWSBor0UR8BWvIq0UJEFFgJxABQcIjL8JEFFwIIg3/tHDudCzHByhpw89PNaKysnZ+99nm+ek3yy8+y9n23ujoiIhKdJugsQEZH4KMBFRAKlABcRCZQCXEQkUApwEZFAZTRkY9nZ2d6xY8eGbFJEJHglJSUV7p5T9fkGDfCOHTuybt26hmxSRCR4ZvZRdc9rCEVEJFAKcBGRQCnARUQC1aBj4CISn8OHD7N9+3YOHjyY7lIkhbKysmjfvj2ZmZkxra8AFwnA9u3badmyJR07dsTM0l2OpIC7s2fPHrZv384555wT0zYaQhEJwMGDB2nbtq3C+yRmZrRt27Ze/2UpwEUCofA++dX3PVaAi4gESmPgIgGa9sqWpL7eXVecV+9t+vfvz5/+9CdOPfXUE56fPHkyLVq0YMyYMUmqLjGHDh3iqquuoqKiggkTJtCuXTtGjRpFZmYmS5Ys4c477+TZZ5+tcfsRI0bwi1/8gosuuqjeba9YsYJTTjmFbt26JfIt1EgBLieP5VPi37b3hOTVcZJzd9ydF198Md2lxGTDhg0cPnyY0tJSAEaNGsWYMWMYPnw4QK3hDTBnzpy4216xYgUtWrRIWYBrCEVE6lRWVsaFF17Iz372MwoLC/n444/p2LEjFRUVANx///2cf/759OnTh3fffTe63dq1a8nPz6dr166MHTuW3NxcAL7++mvGjh3LpZdeSn5+Po888ki17T7xxBPk5+fTqVMnbrrpJgA++ugjiouLyc/Pp7i4mG3btgFQXl7O4MGDufTSS7n00ktZtWoVu3fvZtiwYZSWllJQUMAjjzzC008/zT333MPQoUMpKys7oaYxY8aQl5dHfn4+M2fOBKBXr17RKUCWLl1K165dKSws5LrrrmP//v1A5TQhkyZNorCwkLy8PN555x3KysqYNWsW06ZNo6CggJUrV/LMM8+Qm5tLp06d6NmzZ8Lvi/bARSQm7777Lo899hgPPfTQCc+XlJTw5JNPsmHDBo4cOUJhYSGXXHIJAMOHD2f27Nl069aN8ePHR7eZO3curVu3Zu3atRw6dIju3bvTt2/fE06f27hxI/fffz+rVq0iOzubTz/9FIA77riDm2++mVtuuYVHH32U0aNHs3jxYu68807uuusuevTowbZt2+jXrx+bN29mzpw5TJ06lRdeeAGA1atXc/XVVzNkyBDKysqi7c2ePZsPP/yQDRs2kJGREW3vmIqKCu677z6WLVtG8+bNeeCBB3jwwQeZOHEiANnZ2axfv56HHnqIqVOnMmfOHEaNGnXCcFJeXh4vv/wyZ511Fp9//nnC74kCXERi0qFDB7p06fKN51euXMmgQYNo1qwZAAMGDADg888/Z9++fdHhgxtvvDEaokuXLuXNN9+MDl/s3buXrVu3nhDgf/vb3xgyZAjZ2dkAnHbaaUBlAC9cuBCAm266iV/+8pcALFu2jE2bNkW3/+KLL9i3b1/M39+yZcsYNWoUGRkZJ7R3zJo1a9i0aRPdu3cH4KuvvqJr167R5ddeey0Al1xySbS+qrp3786tt97Kj370o+j6iVCAi0hMmjdvXuOy6k5/q+2G6e7OzJkz6devX63rxHJa3bF1jh49yurVq/n2t79d5zbxtOfuXHHFFSxYsKDa5d/61rcAaNq0KUeOHKl2nVmzZvH666+zZMkSCgoKKC0tpW3btnHVCxoDF5EE9ezZk0WLFvHPf/6Tffv28Ze//AWANm3a0LJlS9asWQPAk08+Gd2mX79+PPzwwxw+fBiALVu2cODAgRNet7i4mKeffpo9e/YARIc0unXrFn2t+fPn06NHDwD69u3L73//++j2xw5axqpv377MmjUrGr5Vh1C6dOnCqlWreO+99wD48ssv2bKl9rOBWrZsecJ/Ae+//z6dO3fmnnvuITs7m48//rheNValPXCRAMVz2l+qFBYWcv3111NQUECHDh34wQ9+EF02d+5cbrvtNpo3b06vXr1o3bo1UHlqXllZGYWFhbg7OTk5LF68+ITX/f73v8/dd9/N5ZdfTtOmTbn44ouZN28eM2bM4Cc/+Qm//e1vycnJ4bHHHgNgxowZ3H777eTn53PkyBF69uzJrFmzYv4+RowYwZYtW8jPzyczM5PbbruNO+64I7o8JyeHefPmccMNN3Do0CEA7rvvPs47r+b34pprrmHIkCE8//zzzJw5k2nTprF161bcneLiYjp16hRzfdWx2v7NSbaioiLXDR0kZU7i0wg3b97MhRdemO4y6m3//v20aNECgN/85jfs2rWL3/3ud2muqnGr7r02sxJ3L6q6rvbARSRllixZwpQpUzhy5AgdOnRg3rx56S7ppKIAF5GUuf7667n++uvTXcZJSwcxRUQCpQAXEQmUAlxEJFAKcBGRQOkgpkiIEjllsjoxnEbZokWL6ORNIVm5cmV0+tjVq1czceJEXnzxRfr37893v/tdmjVrxs0331zttjt37mT06NF1zlhYk+nTpzNy5MjoNAPJVmeAm9mjwNXAbnfPjTx3GvAU0BEoA37k7p+lpEIRkQTMnz//hOljH3nkEcrLy6OXvtemXbt2cYc3VAb4sGHDUhbgsQyhzAOurPLceOCv7v494K+Rr0XkX8D+/fspLi6OTp36/PPPA3DgwAGuuuoqOnXqRG5uLk899RQA48eP56KLLiI/Pz86K19NU8JWbWf48OHR6V2fe+45ABYsWEBeXh65ubmMGzcuun51U73OmTPnhOljBwwYwIEDB+jcuTNPPfUUkydPZurUqQC899579OnTh06dOlFYWMj777//jelmq5sCd8WKFfTq1YshQ4ZwwQUXMHToUNydGTNmsHPnTnr37k3v3r35+uuvufXWW8nNzSUvL49p06Yl/F7UuQfu7q+aWccqTw8EekUePw6sAMYhIie9rKwsFi1aRKtWraioqKBLly4MGDCAl156iXbt2rFkyRKgcobBTz/9lEWLFvHOO+9gZtEpVGuaEvZ49957L61bt+att94C4LPPPmPnzp2MGzeOkpIS2rRpQ9++fVm8eDE9evSocarX1157LTp9LFQOBR2bJ2Xy5MnR9oYOHcr48eMZNGgQBw8e5OjRo+zevTu6vKYpcKHyphEbN26kXbt2dO/enVWrVjF69GgefPBBli9fTnZ2NiUlJezYsYO3334bICnTycZ7EPMMd98FEPl8esKViEgQ3J1f/epX5Ofn06dPH3bs2MEnn3xCXl4ey5YtY9y4caxcuZLWrVvTqlUrsrKyGDFiBAsXLowOJaxevZobb7wRqJwS9rXXXvtGO8uWLeP222+Pft2mTRvWrl1Lr169yMnJISMjg6FDh/Lqq6+eMNVrQUEBjz/+OB999FHM39O+ffvYsWMHgwYNAir/SFUd9li6dClPPPEEBQUFdO7cmT179rB161YALrvsMtq3b0+TJk0oKCg4YZ7xY84991w++OADfv7zn/PSSy/RqlWrmOurScrPQjGzkWa2zszWlZeXp7o5EUmx+fPnU15eTklJCaWlpZxxxhkcPHiQ8847j5KSEvLy8pgwYQL33HMPGRkZ/OMf/2Dw4MEsXryYK6+sOhpbqabpaKs+X9PcTcemei0tLaW0tJRNmzYxd+7cmL+nWOaEOjYF7rE2Pvzww+ge+PHj6TVNJ9umTRveeOMNevXqxR/+8AdGjBgRc301iTfAPzGzMwEin3fXtKK7z3b3IncvysnJibM5EWks9u7dy+mnn05mZibLly+P7unu3LmTZs2aMWzYMMaMGcP69evZv38/e/fupX///kyfPj06dFHTlLDHqzo97GeffUbnzp35+9//TkVFBV9//TULFizg8ssvj2uq1+O1atWK9u3bR4dxDh06xJdffnnCOrFMgVvV8dPJVlRUcPToUQYPHsy9997L+vXrY66vJvGeRvhn4BbgN5HPzydciYjELo2zJw4dOpRrrrmGoqIiCgoKuOCCCwB46623GDt2LE2aNCEzM5OHH36Yffv2MXDgQA4ePIi7Rw/c1TQl7PF+/etfc/vtt5Obm0vTpk2ZNGkS1157LVOmTKF37964O/3792fgwIEA9Z7qtao//vGP/PSnP2XixIlkZmbyzDPP0KTJ/+7jxjIFblUjR47khz/8IWeeeSbTp09n+PDhHD16FIApUxI/FbTO6WTNbAGVByyzgU+AScBi4GngbGAbcJ27f1rDS0RpOllJKU0nKyeBpE4n6+431LCoOL7yREQkGXQpvYhIoBTgIoFoyLtnSXrU9z1WgIsEICsriz179ijET2Luzp49e8jKyop5G01mJRKA9u3bs337dnQtxcktKyuL9u3bx7y+AlwkAJmZmZxzzjnpLkMaGQ2hiIgESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigUoowM3sLjPbaGZvm9kCM8tKVmEiIlK7uAPczM4CRgNF7p4LNAV+nKzCRESkdokOoWQA3zazDKAZsDPxkkREJBZxB7i77wCmAtuAXcBed19adT0zG2lm68xsXXl5efyViojICRIZQmkDDATOAdoBzc1sWNX13H22uxe5e1FOTk78lYqIyAkSGULpA3zo7uXufhhYCHRLTlkiIlKXRAJ8G9DFzJqZmQHFwObklCUiInVJZAz8deBZYD3wVuS1ZiepLhERqUNGIhu7+yRgUpJqERGRetCVmCIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigUoowM3sVDN71szeMbPNZtY1WYWJiEjtMhLc/nfAS+4+xMxOAZoloSYREYlB3AFuZq2AnsCtAO7+FfBVcsoSEZG6JDKEci5QDjxmZhvMbI6ZNa+6kpmNNLN1ZrauvLw8geZEROR4iQR4BlAIPOzuFwMHgPFVV3L32e5e5O5FOTk5CTQnIiLHSyTAtwPb3f31yNfPUhnoIiLSAOIOcHf//8DHZnZ+5KliYFNSqhIRkTolehbKz4H5kTNQPgCGJ16SiIjEIqEAd/dSoCg5pYiISH3oSkwRkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCQe4mTU1sw1m9kIyChIRkdgkYw/8TmBzEl5HRETqIaEAN7P2wFXAnOSUIyIiscpIcPvpwC+BljWtYGYjgZEAZ599doLNSX1Ne2VLzOvedcV5lQ+WT0ms0d4TEtteRGIS9x64mV0N7Hb3ktrWc/fZ7l7k7kU5OTnxNiciIlUkMoTSHRhgZmXAk8C/m9l/JaUqERGpU9wB7u4T3L29u3cEfgz8zd2HJa0yERGplc4DFxEJVKIHMQFw9xXAimS8loiIxEZ74CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigUrKbIT/SuK6Rdm/mgRvybb6gz0xr9v13LYJtVUTvc8SAu2Bi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhIoBbiISKDiDnAz+46ZLTezzWa20czuTGZhIiJSu0Ru6HAE+H/uvt7MWgIlZvaKu29KUm0iIlKLuPfA3X2Xu6+PPN4HbAbOSlZhIiJSu6TcUs3MOgIXA69Xs2wkMBLg7LPPTkZzJ6XqbuHVZdvsuF9vzdkjEynnpFTbrdrWHIn9FmqJqvpe1/Y+13nLuN4TklJDbVJ1y7j61JDKOkKW8EFMM2sBPAf8h7t/UXW5u8929yJ3L8rJyUm0ORERiUgowM0sk8rwnu/uC5NTkoiIxCKRs1AMmAtsdvcHk1eSiIjEIpE98O7ATcC/m1lp5KN/kuoSEZE6xH0Q091fAyyJtYiISD3oSkwRkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUAl5ZZqDWL5lPi3jfO2U4mq7y2jQlLb7clClMjt61jeNm0/Yw0h2bf7gwRv+ZeuLEik3UTbroH2wEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCQW4mV1pZu+a2XtmNj5ZRYmISN3iDnAzawr8AfghcBFwg5ldlKzCRESkdonsgV8GvOfuH7j7V8CTwMDklCUiInUxd49vQ7MhwJXuPiLy9U1AZ3e/o8p6I4Fj90/KBd6Ov9yUygYq0l1EDRprbY21LlBt8WisdYFq6+DuOVWfTOSemFbNc9/4a+Dus4HZAGa2zt2LEmgzZVRb/TXWukC1xaOx1gWqrSaJDKFsB75z3NftgZ2JlSMiIrFKJMDXAt8zs3PM7BTgx8Cfk1OWiIjUJe4hFHc/YmZ3AC8DTYFH3X1jHZvNjre9BqDa6q+x1gWqLR6NtS5QbdWK+yCmiIikl67EFBEJlAJcRCRQKQ1wM/utmb1jZm+a2SIzO7WG9Rr8knwzu87MNprZUTOr8RQgMyszs7fMrNTM1jWy2hq038zsNDN7xcy2Rj63qWG9BuuzuvrAKs2ILH/TzApTWU896uplZnsjfVRqZhMbqK5HzWy3mVV7PUa6+ivG2tLVZ98xs+Vmtjnye3lnNeukp9/cPWUfQF8gI/L4AeCBatZpCrwPnAucArwBXJTKuiLtXgicD6wAimpZrwzITnU99a0tHf0G/CcwPvJ4fHXvZ0P2WSx9APQH/pvK6xa6AK83krp6AS805M9VpN2eQCHwdg3LG7y/6lFbuvrsTKAw8rglsKUx/Jy5e2r3wN19qbsfiXy5hspzxatKyyX57r7Z3d9NdTvxiLG2dPTbQODxyOPHgf+T4vbqEksfDASe8EprgFPN7MxGUFdauPurwKe1rJKO/oq1trRw913uvj7yeB+wGTirympp6beGHAP/CZV/oao6C/j4uK+3883OSScHlppZSWRagMYiHf12hrvvgsofauD0GtZrqD6LpQ/S0U+xttnVzN4ws/82s++nuKZYNfbfx7T2mZl1BC4GXq+yKC39lsil9ACY2TLg36pZdLe7Px9Z527gCDC/upeo5rmknNsYS20x6O7uO83sdOAVM3snsqeQ7tpS0m+11VWPl0lJn1Ujlj5I2c9XLWJpcz2V81vsN7P+wGLgeymuKxbp6K9YpbXPzKwF8BzwH+7+RdXF1WyS8n5LOMDdvU9ty83sFuBqoNgjg0VVpOyS/Lpqi/E1dkY+7zazRVT+e5xwGCWhtpT0W211mdknZnamu++K/Hu4u4bXSEmfVSOWPkjHlA91tnl8ALj7i2b2kJllu3u6J2xqtFNkpLPPzCyTyvCe7+4Lq1klLf2W6rNQrgTGAQPc/csaVmu0l+SbWXMza3nsMZUHZRvLbIrp6Lc/A7dEHt8CfOM/hQbus1j64M/AzZGzBLoAe48NA6VQnXWZ2b+ZmUUeX0bl7+KeFNcVi3T0V0zS1WeRNucCm939wRpWS0+/pfjo7XtUjguVRj5mRZ5vB7xY5QjuFiqP3N/dEEdvgUFU/tU8BHwCvFy1NirPIngj8rGxMdWWjn4D2gJ/BbZGPp+W7j6rrg+AUcCoyGOj8sYj7wNvUcsZRw1c1x2R/nmDygP83RqorgXALuBw5Gfs/zaG/oqxtnT1WQ8qh0PePC7L+jeGftOl9CIigdKVmCIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhKo/wFV99BRk7LM1gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(w_ridge, label=\"ridge coefficients\", bins=20, alpha=.5)\n",
    "plt.hist(w_lasso, label=\"lasso coefficients\", bins=20, alpha=.5)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which coefficients are more sparse? (More 0 coefficients)\n",
    "* A) ridge\n",
    "* B) lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcceeb3a1c0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZU0lEQVR4nO3de3RV5bX38e98IYgCckt8rVABO45UCiGGYAMIKChSpfbQguCVghqpKAgCrw47xAvWC1ROobWYAqI1FW+UWrFVQThYRI4gAcEIXhoh52AJWBCkApH5/pGwD4lJ9s7eO4Qn/D5jZJC9LvuZi8jPlbWfNZe5OyIiEp7/U9cFiIhIfBTgIiKBUoCLiARKAS4iEigFuIhIoBoey8FSU1O9ffv2x3JIEZHgrV27dqe7p1VcfkwDvH379qxZs+ZYDikiEjwz+7Sy5bqEIiISKAW4iEigFOAiIoE6ptfARaTmDh06RFFREV999VVdlyK1rHHjxrRt25aUlJSYtleAixznioqKaNasGe3bt8fM6rocqSXuzq5duygqKqJDhw4x7aNLKCLHua+++orWrVsrvOs5M6N169Y1+k1LAS4SAIX3iaGmP2cFuIhIoHQNXCQwM17fktT3G3/x2TXe59JLL+UPf/gDLVq0KLf8nnvuoWnTpkycODFJ1SXmwIEDXHbZZezcuZM777yTM844g9GjR5OSksLixYsZN24cL7zwQpX733DDDUyYMIFOnTrVeOzly5fTqFEjevbsmcghVEsBfqJY9mDs2154Z+3VIUFzd9ydV155pa5Licm6des4dOgQ+fn5AIwePZqJEycycuRIgGrDG2DOnDlxj718+XKaNm1aqwGuSygiUq3CwkLOOeccbr75ZjIzM9m2bRvt27dn586dADzwwAN07NiRiy66iM2bN0f2e+edd0hPT6dHjx5MmjSJzp07A/D1118zadIkunfvTnp6Oo8//nil4z711FOkp6fTtWtXrr32WgA+/fRT+vfvT3p6Ov3792fr1q0AFBcX85Of/ITu3bvTvXt3Vq5cyY4dO7jmmmvIz88nIyODxx9/nOeee4777ruPq6++msLCwnI1TZw4kS5dupCens6sWbMAuOCCCyLtP1577TV69OhBZmYmQ4cOZd++fUBpi5ApU6aQmZlJly5d+OCDDygsLGT27NnMmDGDjIwM3nzzTZ5//nk6d+5M165d6dOnT1J+NjoDF5GoNm/ezBNPPMFjjz1WbvnatWtZsGAB69ato6SkhMzMTLp16wbAyJEjyc3NpWfPntxxxx2RfebOnUvz5s155513OHDgAL169WLAgAHlps5t2rSJBx54gJUrV5Kamsrnn38OwC233MJ1113HiBEjmDdvHmPHjmXRokWMGzeO8ePHc/7557N161YuueQSCgoKmDNnDtOnT+fll18GYNWqVQwaNIghQ4ZQWFgYGS83N5e///3vrFu3joYNG0bGO2Lnzp1MnTqVJUuW0KRJEx5++GEeffRR7r77bgBSU1N59913eeyxx5g+fTpz5sxh9OjR5S4ndenShVdffZU2bdqwe/fupPxcFOAiElW7du3Izs7+xvI333yTwYMHc8oppwBw+eWXA7B792727t0buXxw1VVXRUL0tddeY8OGDZHLF3v27OHDDz8sF+BvvPEGQ4YMITU1FYBWrVoBpQG8cOFCAK699lomT54MwJIlS3j//fcj+3/xxRfs3bs35uNbsmQJo0ePpmHDhuXGO+Ltt9/m/fffp1evXgAcPHiQHj16RNb/+Mc/BqBbt26R+irq1asXP/3pT7niiisi2ydKAS4iUTVp0qTKdZVNfavuYenuzqxZs7jkkkuq3SaWKXVHtjl8+DCrVq3i5JNPjrpPPOO5OxdffDHPPPNMpetPOukkABo0aEBJSUml28yePZvVq1ezePFiMjIyyM/Pp3Xr1nHVe4SugYtI3Pr06cMf//hH/vWvf7F3717+/Oc/A9CyZUuaNWvG22+/DcCCBQsi+1xyySX89re/5dChQwBs2bKFL7/8stz79u/fn+eee45du3YBRC5p9OzZM/JeeXl5nH/++QAMGDCAX//615H9j3xoGasBAwYwe/bsSPhWvISSnZ3NypUr+eijjwDYv38/W7ZUPxuoWbNm5X4L+Pjjj/n+97/PfffdR2pqKtu2batRjZXRGbhIYOKZ9ldbMjMzGTZsGBkZGbRr147evXtH1s2dO5cbb7yRJk2acMEFF9C8eXOgdGpeYWEhmZmZuDtpaWksWrSo3Pt+73vf46677qJv3740aNCAc889l/nz5zNz5kxGjRrFtGnTSEtL44knngBg5syZjBkzhvT0dEpKSujTpw+zZ8+O+ThuuOEGtmzZQnp6OikpKdx4443ccsstkfVpaWnMnz+fK6+8kgMHDgAwdepUzj676p/FD3/4Q4YMGcKf/vQnZs2axYwZM/jwww9xd/r370/Xrl1jrq8qVt2vOsmWlZXleqBDHdE0wmAVFBRwzjnn1HUZNbZv3z6aNm0KwEMPPcT27dv51a9+VcdVHf8q+3mb2Vp3z6q4rc7ARaRWLF68mAcffJCSkhLatWvH/Pnz67qkeidqgJvZPGAQsMPdOx+1/FbgFqAEWOzuk2utShEJzrBhwxg2bFhdl1GvxfIh5nxg4NELzOxC4EdAurt/D5ie/NJERKQ6UQPc3VcAn1dY/DPgIXc/ULbNjlqoTUREqhHvNMKzgd5mttrM/tPMule1oZnlmNkaM1tTXFwc53AiIlJRvAHeEGgJZAOTgOesilnw7p7r7lnunpWWlhbncCIiUlG8s1CKgIVeOgfxv8zsMJAK6BRbpLbVZEpoLGKYNmpmTJgwgV/+8pcATJ8+nX379nHPPfdwzz338Mgjj1BYWMhpp50GQNOmTSPNno7+/ojNmzdz0003sXv3bg4cOEDv3r3Jzc1N7nGVKS4uZtCgQRw8eJCZM2fy2Wefcffdd3P66aczbdo0nnrqKWbOnFnl/lW1zo3FokWLOPvss+NqRxuLeM/AFwH9AMzsbKARsDNJNYnIceakk05i4cKFkQ6EFaWmpkbCPRZjx45l/Pjx5OfnU1BQwK233pqsUr9h6dKlfPe732XdunX07t2buXPn8thjj7Fs2TKysrKqDW+AV155Ja7whtIAP7pHS7JFDXAzewZYBXQ0syIzux6YB5xlZhuBBcAIP5Z3BInIMdWwYUNycnKYMWNGpetHjRrFs88++41b0Kuyfft22rZtG3ndpUuXSrd75JFH6NKlC127do10NMzPzyc7O5v09HQGDx7MP//5T6D0VvWBAwfSrVs3evfuzQcffEB+fj6TJ0/mlVdeISMjg3vvvZe//e1vjB49mkmTJrF8+XIGDRoElN54NHLkyEhL2RdffBGgXOvcp59+mvPOO4+MjAxuuukmvv76a6D0t4y77rqLrl27kp2dzT/+8Q/eeustXnrpJSZNmkRGRgYff/wxM2fOpFOnTqSnpzN8+PCY/q6qE8sslCvd/VvunuLubd19rrsfdPdr3L2zu2e6+xsJVyIix7UxY8aQl5fHnj17vrGuadOmjBo1KuY7LcePH0+/fv34wQ9+wIwZMyptr/qXv/yFRYsWsXr1atavXx/pPHjdddfx8MMPs2HDBrp06cK9994LQE5ODrNmzWLt2rVMnz6dm2++mYyMDO677z6GDRtGfn4+U6ZMISsri7y8PKZNm1ZuvPvvv5/mzZvz3nvvsWHDBvr161dufUFBAc8++ywrV64kPz+fBg0akJeXB8CXX35JdnY269evp0+fPvzud7+jZ8+eXH755UybNo38/Hy+853v8NBDD7Fu3To2bNhQo1v9q6JmViISk1NPPZXrrruuyksOY8eO5cknn+SLL76I+l4jR46koKCAoUOHsnz5crKzsyM9Ro5YsmQJI0eOjLSqbdWqFXv27GH37t307dsXgBEjRrBixQr27dvHW2+9xdChQyNnx9u3b6/R8S1ZsoQxY8ZEXrds2bLc+qVLl7J27Vq6d+9ORkYGS5cu5ZNPPgGgUaNGkTP5bt26les1frT09HSuvvpqnn766Ujr2kToVnoRidltt91GZmZm5JFkR2vRogVXXXXVNx76UJUzzjiDUaNGMWrUKDp37szGjRsjD4OA2FvKQmk72RYtWtS4C+HRYmkpO2LECB588JsfIqekpET2ra6l7OLFi1mxYgUvvfQS999/P5s2bUooyHUGLiIxa9WqFVdccQVz586tdP2ECRN4/PHHqwywI/76179G2sl+9tln7Nq1izZt2pTbZsCAAcybN4/9+/cDpS1emzdvTsuWLXnzzTcB+P3vf0/fvn059dRT6dChA88//zxQGrbr16+v0bFVbEl75Nr6Ef379+eFF15gx44dkXo+/fTTat/z6Jayhw8fZtu2bVx44YU88sgj7N69+xuzc2pKZ+AioanjbpG33357uaA7WmpqKoMHDy73Yef+/fvLfWA5YcIEioqKGDduHI0bNwZg2rRpnH766eXea+DAgeTn55OVlUWjRo249NJL+cUvfsGTTz7J6NGj2b9/P2eddVakpWxeXh4/+9nPmDp1KocOHWL48OE1atn685//nDFjxtC5c2caNGjAlClTyj05p1OnTkydOpUBAwZw+PBhUlJS+M1vfkO7du2qfM/hw4dz4403MnPmTBYsWMD111/Pnj17cHfGjx8f9+yWI9RO9kShdrLBCrWdrMSnJu1kdQlFRCRQCnARkUApwEUCoPvkTgw1/TkrwEWOc40bN2bXrl0K8XrO3dm1a1fkg91YaBaKyHGubdu2FBUVoXbM9V/jxo3LzdiJRgEucpxLSUmhQ4cOdV2GHId0CUVEJFAKcBGRQCnARUQCpQAXEQlULA90mGdmO8oe3lBx3UQzczNLrZ3yRESkKrGcgc8HBlZcaGbfBi4Gtia5JhERiUEsT+RZAVT2nKQZwGRAdxeIiNSBuK6Bm9nlwH+7e9SGu2aWY2ZrzGyNbkQQEUmeGge4mZ0C3AXcHcv27p7r7lnunpWWllbT4UREpArxnIF/B+gArDezQqAt8K6ZnV7tXiIiklQ1vpXe3d8DTjvyuizEs9x9ZxLrEhGRKGKZRvgMsAroaGZFZnZ97ZclIiLRRD0Dd/cro6xvn7RqREQkZroTU0QkUApwEZFAKcBFRAKlABcRCZQCXEQkUApwEZFAKcBFRAKlABcRCZQCXEQkUApwEZFAKcBFRAKlABcRCZQCXEQkUApwEZFAKcBFRAIVywMd5pnZDjPbeNSyaWb2gZltMLM/mlmLWq1SRES+IZYz8PnAwArLXgc6u3s6sAW4M8l1iYhIFFED3N1XAJ9XWPaau5eUvXyb0gcbi4jIMZSMa+CjgL8k4X1ERKQGEgpwM7sLKAHyqtkmx8zWmNma4uLiRIYTEZGjxB3gZjYCGARc7e5e1XbunuvuWe6elZaWFu9wIiJSQdSn0lfGzAYC/w/o6+77k1uSiIjEIpZphM8Aq4COZlZkZtcDvwaaAa+bWb6Zza7lOkVEpIKoZ+DufmUli+fWQi0iIlIDuhNTRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAqUAFxEJVCxP5JlnZjvMbONRy1qZ2etm9mHZny1rt0wREakoljPw+cDACsvuAJa6+78BS8tei4jIMRQ1wN19BfB5hcU/Ap4s+/5J4N+TW5aIiEQT7zXw/+vu2wHK/jytqg3NLMfM1pjZmuLi4jiHExGRimr9Q0x3z3X3LHfPSktLq+3hREROGPEG+D/M7FsAZX/uSF5JIiISi3gD/CVgRNn3I4A/JaccERGJVSzTCJ8BVgEdzazIzK4HHgIuNrMPgYvLXouIyDHUMNoG7n5lFav6J7kWERGpAd2JKSISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFKKMDNbLyZbTKzjWb2jJk1TlZhIiJSvbgD3MzaAGOBLHfvDDQAhierMBERqV6il1AaAiebWUPgFOB/Ei9JRERiEXeAu/t/A9OBrcB2YI+7v1ZxOzPLMbM1ZramuLg4/kpFRKScRC6htAR+BHQAzgCamNk1Fbdz91x3z3L3rLS0tPgrFRGRchK5hHIR8Hd3L3b3Q8BCoGdyyhIRkWgSCfCtQLaZnWJmRulT6guSU5aIiESTyDXw1cALwLvAe2XvlZukukREJIqGiezs7lOAKUmqRUREakB3YoqIBEoBLiISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFKKMDNrIWZvWBmH5hZgZn1SFZhIiJSvYSeyAP8Cviruw8xs0bAKUmoSUREYhB3gJvZqUAf4KcA7n4QOJicskREJJpELqGcBRQDT5jZOjObY2ZNKm5kZjlmtsbM1hQXFycwnIiIHC2RAG8IZAK/dfdzgS+BOypu5O657p7l7llpaWkJDCciIkdLJMCLgCJ3X132+gVKA11ERI6BuAPc3T8DtplZx7JF/YH3k1KViIhElegslFuBvLIZKJ8AIxMvSUREYpFQgLt7PpCVnFJERKQmdCemiEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoFSgIuIBCrhADezBmUPNX45GQWJiEhsknEGPg4oSML7iIhIDSQU4GbWFrgMmJOcckREJFaJPhPzP4DJQLOqNjCzHCAH4Mwzz0xwuLo34/Ut1a4ff/HZx6iSWrTswZptf+GdtVOHiFQr7jNwMxsE7HD3tdVt5+657p7l7llpaWnxDiciIhUkcgmlF3C5mRUCC4B+ZvZ0UqoSEZGo4g5wd7/T3du6e3tgOPCGu1+TtMpERKRamgcuIhKoRD/EBMDdlwPLk/FeIiISG52Bi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhKopEwjPBai9SCBetKHJAbH4u9i1Se7om7T46zWCY0hIonRGbiISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gESgEuIhKoRJ6J+W0zW2ZmBWa2yczGJbMwERGpXiJ3YpYAt7v7u2bWDFhrZq+7+/tJqk1ERKqRyDMxt7v7u2Xf7wUKgDbJKkxERKqXlF4oZtYeOBdYXcm6HCAH4Mwzz0zGcHGLpYfIiSKWXicxW/ZgzGO8fWZO5Ptk9K6J9jOtdIwq6q1yjJKfRN3mROnDI8eXhD/ENLOmwIvAbe7+RcX17p7r7lnunpWWlpbocCIiUiahADezFErDO8/dFyanJBERiUUis1AMmAsUuPujyStJRERikcgZeC/gWqCfmeWXfV2apLpERCSKuD/EdPe/AZbEWkREpAZ0J6aISKAU4CIigVKAi4gESgEuIhIoBbiISKAU4CIigVKAi4gEKinNrI4XITSrCqHG2pS9Nfd/XyxrXf3GF96Z8HiV/X1nb61ZI69scqNvFO1YKpOE4zuu1aRpWIh/FzVsilYbx6gzcBGRQCnARUQCpQAXEQmUAlxEJFAKcBGRQCnARUQCpQAXEQmUAlxEJFCJPhNzoJltNrOPzOyOZBUlIiLRJfJMzAbAb4AfAJ2AK82sU7IKExGR6iVyBn4e8JG7f+LuB4EFwI+SU5aIiERj7h7fjmZDgIHufkPZ62uB77v7LRW2ywFyyl52BDbHX26tSAV21nURx4COs37RcdYv0Y6znbunVVyYSDOryh5o/I3/G7h7LsTSDahumNkad8+q6zpqm46zftFx1i/xHmcil1CKgG8f9bot8D8JvJ+IiNRAIgH+DvBvZtbBzBoBw4GXklOWiIhEE/clFHcvMbNbgFeBBsA8d9+UtMqOneP28k6S6TjrFx1n/RLXccb9IaaIiNQt3YkpIhIoBbiISKAU4ICZDTWzTWZ22Mzq1ZSlE6XdgZnNM7MdZraxrmupTWb2bTNbZmYFZf/NjqvrmmqDmTU2s/8ys/Vlx3lvXddUm8ysgZmtM7OXa7KfArzURuDHwIq6LiSZTrB2B/OBgXVdxDFQAtzu7ucA2cCYevozPQD0c/euQAYw0Myy67akWjUOKKjpTgpwwN0L3P14u0M0GU6YdgfuvgL4vK7rqG3uvt3d3y37fi+l/+jb1G1Vyeel9pW9TCn7qpczLsysLXAZMKem+yrA67c2wLajXhdRD/+xn6jMrD1wLrC6jkupFWWXFfKBHcDr7l4vjxP4D2AycLimO54wAW5mS8xsYyVf9fKMtExM7Q4kPGbWFHgRuM3dv6jremqDu3/t7hmU3uV9npl1ruOSks7MBgE73H1tPPsn0gslKO5+UV3XUAfU7qAeMrMUSsM7z90X1nU9tc3dd5vZcko/46hvH1L3Ai43s0uBxsCpZva0u18Ty84nzBn4CUrtDuoZMzNgLlDg7o/WdT21xczSzKxF2fcnAxcBH9RpUbXA3e9097bu3p7Sf59vxBreoAAHwMwGm1kR0ANYbGav1nVNyeDuJcCRdgcFwHOBtjuIysyeAVYBHc2syMyur+uaakkv4Fqgn5nll31dWtdF1YJvAcvMbAOlJyKvu3uNptidCHQrvYhIoHQGLiISKAW4iEigFOAiIoFSgIuIBEoBLiISKAW4iEigFOAiIoH6/3dIukv9Pj1NAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(w_ridge, label=\"ridge coefficients\", bins=20, alpha=.5)\n",
    "plt.hist(w_nonneg, label=\"NNLS coefficients\", bins=20, alpha=.5)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which coefficients are more sparse? (More 0 coefficients)\n",
    "* A) ridge\n",
    "* B) nonnegative least squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcceec56250>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAywElEQVR4nO3deVxVdf748dcbxEBccCtFVLTMLEUlLLNlKk1KympmtGWaFqdsvm20UTpTDmM/R8rKoWZaLMtqbMwatQVTK9tX19RS0pAKxdwSFSEWP78/7gW53Asc7r3nbryfj4cP7v2cc+59e4D75nzO5/P+iDEGpZRSqq6oYAeglFIq9GhyUEop5UaTg1JKKTeaHJRSSrnR5KCUUspNq2AH4A9dunQxycnJwQ5DKaXCyqpVq3YbY7p62hYRySE5OZmVK1cGOwyllAorIvJDQ9u0W0kppZQbTQ5KKaXcaHJQSinlJiLuOXhSWVlJUVER5eXlwQ4loGJjY0lKSiImJibYoSilwljEJoeioiLatWtHcnIyIhLscALCGMOePXsoKiqiT58+wQ5HKRXGIrZbqby8nM6dO7eYxAAgInTu3LnFXS0pFUx5BXmMfm00KS+kMPq10eQV5AU7JL+I2CsHoEUlhhot8f+sVLDkFeSR/Vk25dWOP8iKS4vJ/iwbgIy+GUGMzHcRe+WglFJ2y12dW5sYapRXl5O7OjdIEfmPJocAGjNmDPv27XNrz87O5uGHHw58QEopn+wo3dGs9nCiySFAjDG89dZbJCQkBDsUpZSfdIvv1qz2cKLJwWnRmm2cnrOcPpPyOD1nOYvWbPP5NQsLCxkwYAA33XQTqampREdHs3v3bgCmTZtG//79GTVqFPn5+bXHrFixgpSUFE477TSysrIYOHAgANXV1WRlZTFs2DBSUlJ4+umnfY5PKeWbzNRMYqNjXdpio2PJTM0MUkT+o8kBR2KYvGA92/aVYYBt+8qYvGC9XxJEfn4+V199NWvWrKF3794ArFq1innz5rFmzRoWLFjAihUrave/7rrreOqpp/j888+Jjo6ubZ89ezYdOnRgxYoVrFixgmeeeYatW7f6HJ9SynsZfTPIHpFN9/juCEL3+O5kj8gO+5vREOGjlayasTSfsspql7ayympmLM3nkqE9fHrt3r17M3z4cJe2jz/+mEsvvZQ2bdoAMHbsWAD27dvHgQMHGDFiBABXXnklb731FgDLli1j3bp1vPbaawCUlJSwefNmnc+gVJBl9M2IiGRQnyYHYPu+sma1N0d8fLzHdk9DTo0xDb6OMYbHH3+c9PR0n2NSSqmmaLcSkJgQ16x2X5111lksXLiQsrIyDhw4wJtvvglAx44dadeuHV988QUA8+bNqz0mPT2dJ598ksrKSgC+++47SktLbYlPKaX0ygHISu/P5AXrXbqW4mKiyUrvb8v7paamctlllzFkyBB69+7NmWeeWbtt9uzZ3HDDDcTHx3P22WfToUMHAK6//noKCwtJTU3FGEPXrl1ZtGiRLfEppZQ01pURLtLS0kz9xX42btzIgAEDLL/GojXbmLE0n+37ykhMiCMrvb/P9xu8cfDgQdq2bQtATk4OxcXF5OY2b0JNc//vSqmWSURWGWPSPG3TKwenS4b2CEoyqC8vL4/p06dTVVVF7969mTNnTrBDUkq1QJocQsxll13GZZddFuwwlFItXFBvSIvIcyKyU0Q21GnrJCLviMhm59eOwYxRKaVaomCPVpoDnF+vbRLwnjGmH/Ce87lSSqkACmpyMMZ8BOyt13wx8ILz8QvAJYGMSSmlVPCvHDw5xhhTDOD8erSnnURkooisFJGVu3btCmiASikV6UIxOVhijJlljEkzxqR17do12OF4VDMkVSml7PD593vYffBXW147FJPDzyLSHcD5dWeQ41FKqZDy2qoikiflccUzX/DUB9/b8h6hmBzeAK5xPr4GeD0g77puPswcCNkJjq/r5vvtpQ8ePMjIkSNJTU1l0KBBvP66479UWlpKRkYGgwcPZuDAgbzyyisATJo0iRNPPJGUlBTuvvtuAH744QdGjhxJSkoKI0eO5Mcff/RbfEqp8LB8088kT8rj7le/BqBdbCtuPuc4W94rqPMcROS/wNlAFxEpAv4G5ADzReRPwI/AONsDWTcf3rwNKp2F9kp+cjwHSBnv88vHxsaycOFC2rdvz+7duxk+fDhjx45lyZIlJCYmkpfnWJC8pKSEvXv3snDhQjZt2oSI1K4cd8stt3D11VdzzTXX8Nxzz3Hbbbdp+QylWohVP+zld09+7tL28T3n0LNTG9veU8tngONKoeQn9/YOPeGODe7tFrVt25aDBw9SWVnJHXfcwUcffURUVBT5+fls3bqV/fv3k56ezvjx47nwwgs588wzqaqq4uSTTyYtLY2MjAwuvPBCWrduTZcuXSguLiYmJobKykq6d+9eu3BQfVo+Q6nIsPnnA5w38yOXtiW3n8kJ3dr7peRPY+UzQrFbKfBKiprX3kxz585l165drFq1irVr13LMMcdQXl7O8ccfz6pVqxg0aBCTJ09m6tSptGrViq+++orf/e53LFq0iPPPrz8NxMFTyW+lVGTYvq+M5El5Lolh/o2nUZiTUZsY7FqgrIaWzwDokNTAlUOSX16+pKSEo48+mpiYGN5//31++OEHALZv306nTp246qqraNu2LXPmzOHgwYMcOnSIMWPGMHz4cI47ztGfOGLECObNm8cf//hH5s6dyxlnnOGX2JRSoeOX0gpGPfohe0orattm/fFkRp/kuia1nQuU1dDkADByius9B4CYOEe7H/zhD3/goosuIi0tjSFDhnDCCScAsH79erKysoiKiiImJoYnn3ySAwcOcPHFF1NeXo4xhpkzZwLw2GOPMWHCBGbMmEHXrl15/vnn/RKbUir4DlVU8dsnPmPTjgO1bQ/+bhCXDevlcX87FyirockBjtx0fm+qoyupQ5IjMfh4M/rgwYMAdOnShc8//9xte3JysseV3b766iuP+y5fvtyneJRSoaWsopqx//qEzTsP1rZlpfdvcgRSYkIcJ+9/h3tazSdRdrPddOGhqvGsan+e32LT5FAjZbxfRiYppVRTqg8bjv3LYpe2a0ck87eLTrR0P/GfJ25m4KpniRNH91OS7ObBmGfZcGIycK5fYtTkoJRSAWKMoc/kxW7tiR1iGdIzwfJAk2HfPw5S4dIWJxWOdm70R6iaHJRSKhDOeHA5Rb80cK+gpJzJC9YDWLuhbPMIS9ChrEopZasJc1aQPCnPJTF0ax/rtl/NaCNLGhpJ6acRlqDJQSmlbPH/3vqW5El5LN90pDzc25lnUpiTwc/7yz0eY3m00cgpjhGVdflxhCVot5JSSvnVS58Xcv/r37i2/ekUzux3pHp0YkIc2zwkgsSEOLc2j2waYVmXXjnYSES46667ap8//PDDZGdnA5CdnU2bNm3YufPIXxV1S3x7Kvedn5/P2WefzZAhQxgwYAATJ060L3ilVLO8862jKF7dxDDj9ykU5mS4JAZwDFeNi4l2aYuLiSYrvb/1N0wZ7yjvk73P8dXPoy01OdjoqKOOYsGCBQ3WQOrSpQuPPPKI5de77bbbuOOOO1i7di0bN27k1ltv9VeoSikvLf1mB8mT8rjhxSP13W4b2Y/CnAzGpfX0eMwlQ3sw/beD6JEQhwA9EuKY/ttBfpvd7A/areSUV5BH7upcdpTuoFt8NzJTM8nom+HTa7Zq1YqJEycyc+ZMpk2b5rZ9woQJzJkzh3vvvZdOnTo1+XrFxcUkJR254TRo0CCf4lNKeW99UQkX/esTl7axgxN57Iqhlo6/ZGiPkEoG9emVA47EkP1ZNsWlxRgMxaXFZH+WTV5Bns+vffPNNzN37lxKSkrctrVt25YJEyaQm5tr6bXuuOMOzj33XC644AJmzpxZW85bKRU4O0rKSZ6U55YYCnMyLCeGcKDJAchdnUt5tevogfLqcnJXW/vQbkz79u25+uqreeyxxzxuv+2223jhhRfYv39/k6913XXXsXHjRsaNG8cHH3zA8OHD+fVXe5YIVEq5Kv21iuRJeQyf/p5Le2FOBoU5vvUyhCLtVgJ2lO5oVntz3X777aSmpnLddde5bUtISODKK6/kiSeesPRaiYmJTJgwgQkTJjBw4EA2bNjAySef7Jc4lVLuPJW6ACj4xxiionwonb9uvq2jjXylyQHoFt+N4tJij+3+0KlTJ8aPH8/s2bOZMGGC2/Y777yTYcOGUVVV1ejrLFmyhJEjRxITE8OOHTvYs2cPPXqEbp+lUuEueZJ71/LGqecT1zraw97NYPPqk/6g3UpAZmomsdGuMxZjo2PJTM3023vcddddjY5auvTSS126iA4dOkRSUlLtv0cffZRly5YxcOBABg8eTHp6OjNmzKBbN/8kMBUYeQV5jH5tNCkvpDD6tdF+ua+l/C95Up5bYvjyLyMpzMnwPTGA44qhst48h8oyR3uI0GVCnewYrRQsukxoaKoZ+FD3/lZsdCzZI7LD9mct0vSZnEf9j8S8287gpMQODR7j1XKd2QmAp89eccxbCJDGlgnVbiWnjL4Z+guqbNXYwAf92Quu3z35Gat++MWlbfY1aYwccEyjx9Us11mzKlvNcp3QRAE9m1ef9AftVlIqQOwe+KCa7/5FG0ielOeSGO7LGEBhTkaTiQEaX66zUQGojeSriL5yMMZYro8eKSKhmzBS2T3wQVn34ueFTKlX/+iytJ48+PuUZr2O18t1BqA2kq8iNjnExsayZ88eOnfu3GIShDGGPXv2EBvrXg5YBV9maqbHew7+HPigGvdB/k6ufX6FS9uJ3duzOPNMr17PpwJ6Ib76ZMQmh6SkJIqKiti1a1ewQwmo2NhYlxIbKnTU3FeIlIEP4SR/xwHS//mRW7uvk9ey0vu73HMALwrohaiIHa2klFK7DvzKsGnvurX7c0azV6OVQoSOVlJKtSjlldWccP8St3Y7ylyEegE9b2lyUEpFjMOHDX09lLrYMu0CWkXr4Mzm0OSglIoInkpdrM8eTbvYmCBEE/40OSilwpqnpPDJveeQ1LFNEKKJHJoclFJhKSV7KfvLXYtVLrr5dIb0TAhOQBFGk4NSKqz8cfaXfLzZtYjlv69MJSOle5AiikyaHJRSYeEfizcy66MCl7a7Rx/PLef2C1JEkS1kk4OIFAIHgGqgqqGxuEqpyDZ/xU/c8791Lm0XpnTnX1emBimiliFkk4PTOcYYz4sgKKUi2uff7+GKZ75waevduQ0fZp0TpIhallBPDkqpFqZg10HOfeRDt/ZIXKc5lIVycjDAMhExwNPGmFnBDkgpZZ9fSisY+sA7bu2aFIIjlJPD6caY7SJyNPCOiGwyxtRWzhKRicBEgF69egUrRqWUjyqqDnP8fW+7tW+dPqbFVFQORSGbHIwx251fd4rIQuAU4KM622cBs8BReC8oQSqlvGaMoc9k91IXm6ddQIyWugi6kEwOIhIPRBljDjgfjwZCZ+VtpZRPPM1q/nrKaDq00VIXoSIkkwNwDLDQeUnZCnjZGONeYlEpFVY8JYUPs86md+f4IESjGhOSycEYUwAMDnYcSqlGrJtveZlLT0nhqatO5vyBukRqqArJ5KCUL/IK8nS1Nbutmw9v3gaVziUyS35yPAeXBDFs2rvsOvCry6F3jDqezFE6qznUWUoOItIb6GeMeVdE4oBWxpgD9oamVPPlFeS5rNNcXFpM9mfZAJog/Om9qUcSQ43KMkd7ynhumruKxet3uGweNeBonr1mWACDDJBmXEGFkyaTg4jcgGPIaCfgWCAJeAoYaW9oSjVf7urc2sRQo7y6nNzVuZoc/KmkyGPzU3uGkFOvC6ldbCvWZ6cHIqrAs3gFFY6sXDncjGMY6ZcAxpjNzrkHSoWcHaU7mtWuvNQhyfFB6PRudSrXV97ttpunCWzhvOaymyauoMKZleTwqzGmomYyioi0wjF7WamQ0y2+G8WlxR7blR+NnAJv3sZ3v3ZidMUMt80NzWpetGYbkxesp6yyGoBt+8qYvGA9QHgmiAauoBpsDyNWZpp8KCJ/AeJE5DzgVeBNe8NSyjuZqZnERse6tMVGx5KZmhmkiCLT3mMvIfnAbLfEUJiT0Wi5ixlL82sTQ42yympmLM23JU7bdUhqXnsYsZIcJgG7gPXAjcBi4D47g1LKWxl9M8gekU33+O4IQvf47mSPyNb7DX5SUXWY5El5pHqogRQXE82iNdsaPX77vrJmtYe8kVMgJs61LSbO0R7mmuxWMsYcBp5x/lMq5GX0zdBk4GcNlbqoq+YKoLHuocSEOLZ5SASJCXEe9g4DNfcVWuhopa14uMdgjOlrS0RKqZDiaQJbQ5q6AshK7+9yzwEcVxxZ6f29ji/oUsZHRDKoz8oN6borsMUC43AMa1VKRTBPSWHhTSMY2qsjp+cs9+oKoOaqImJGK0UwK91Ke+o1/VNEPgHCv1NNKeXGU1L452VDXD7AfbkCuGRoD00GYcBKt1LdhVqjcFxJtLMtIqVUUIx85AO+31Xq0nbjb/oy+YIBbvvqFUDks9Kt9Eidx1VAIRB5HWxKtVD3vraOV1b+5NJ2ap9OvHLjaY0ep1cAkc1Kt5Ku5q1UBHrx80KmvP6NW7suy6mgkeQgInc2dqAx5lH/h6OUstsnm3dz1ewv3do1Kai6Grty0PsKSkWQrbtLOefhD9zaNSkoTxpMDsaYvwcyEKXs0tLXdygpq2Tw35e5tWtSUI2xMlopFvgTcBKOeQ4AGGMm2BiXUn7Rktd3qD5sOPYv7rOav//HGKKjJAgRqXBiZbTSS8AmIB2YCvwB2GhnUEr5S0td38HTXIUNf0+n7VG6+KOyxspPynHGmHEicrEx5gUReRlYandgSvlDS1vfwVNS+Piec+jZqU0QolHhzEpyqHR+3SciA4EdQLJtESnlRy1lfQdPSWHWH09m9EmR9f9UgWMlOcwSkY7A/cAbQFvnY6VCXmZqpss9B7B/fYdA3gD3lBRuH9WP20cdb8v7+UWErrkcaawkh+eNMdXAh4BWYlVhpeZDOVAf1pZugPvhw7HfXxdTWe1aLPmM47rwn+tP9fn/YKsIXnM50ogxja/4KSI/AkuAV4DlpqkDgiAtLc2sXLky2GEoxejXRnvsxuoe351lv1/m/uEIjsVhLnrM0ofjhDkrWL5pp1t72AxLnTnQZe3pWh16wh0bAh9PCyciq4wxaZ62Wbly6A9cBNwMPCcibwLzjDGf+DFGpULKojXbvCoq1+QNcC8XpP/3+1s8LqUZNkmhRgSvuRxprNRWKgPmA/Od9x5ycXQxRdscm1JBsWjNNpdy1Nv2lTF5wXqAJhNEkzfAm/nh+P6mnVw3Z4Vbe9glhRodkhq4cgj/NZeby9d7U3bf27KyhjQi8hsReQJYjWMinHYOqog1Y2m+yzoFcGQJzKZkpmYSGx3r0uZyA9zigvQFuw6SPCnPLTEU5mRYSwzr5ju6cLITHF/XzW/6mECI4DWXm6Pm3lRxaTEGU3tvKq/A2qp7vh5vhdVlQtfiuHrIMsaUNn6EUuGtoaUum1oCE47cdJ7+xaOUVOzkcGUClF5EZckQxw4jp3i+5+D8cDxQXsmgbB9LXYTyTd8IXnO5OXydnBmIyZ1W7jkMNsbs98u7KUXo1zpKTIjzagnMGkdvKOKV/EK6s5vt5hAPVf1Sp1vK84fj4YHj6OthWOrmaRcQE23pAv8IL+9rBEyErrncHL5OzgzE5E4r9xw0MSi/CYdaR1np/flk4RPczjwSZTfbTRf+yeWckX5T0wevm8/A1fcRJxUAJMlucmKehUqYsbQ1l0R/6pYYkl+OB1xrIK28bxRd2h7lVfympAhPlZMaaleB5+vkzEBM7mzmnyRK+aaxy+FQcUn0p+TEPEtS1G6iBJKiHB/wl0R/2uSxh96eQhwVLm1tpIJ7Ws0nbf87ju6dkp8AQ/LPDzoTwxFv3XoGhTkZXicGgJ/p0qx2FXhN3puy+XgrtAqXCqiwqHX03lSWxkaR2zGRHa2i6VZVTeYv+8iw0C0TW+b5/5Eoe5jc+lWoLCO5/GW37bmXD+HiIf5ZcnN6xTimxzxLGzmSpA6Z1kyvHEfopOCWzdfJmYGY3KkrwamACodaR3lVe8nu0onyKMeFdXFMK7K7dILde2nqV2/74c4kRe0mL74NuR0TapPLVXurmbJ3ptv+10fncV/MyzBkn9/iX9n+PCbth3tazSdR9rDddOahqvGsan+e395D+S6jb4ZPH+a+Ht8UKyvB9QeG4airBI4JcR/ZFpGTiJyPY05FNPCsMSbH7vdUNqlTLiKzaxLZ7VpTbiprN9td66i55SpyO3eivN56B+VRUeR27tRkcni29VWktH6O6V3a1yaX77Y8TP2BmvGU8U3snwA4FNcdf9ZMzUrvz+QFFbxRcUZtW1xMNNPT+/vxXVSka3IlOBFZBqQaYw44n2cDr9oZlIhEA/8GzgOKgBUi8oYx5ls731fZoN6wyoxdP8GvCeR268mOyv32j1byYljnjmjPt20baq9rSMZEpq16m/KoMkoL/8zhsmS3fQpjr6x9fMi05qHKy8hu8pWtq5mo580Mb6VqWLnn0Atc7rBVYH/J7lOALcaYAgARmQdcDGhyCDcehlVm7N9HhrQLTC0dL4Z1dovv3kDXV/cm3+6SoT3IeudcKvee6bat7QmTOPe7Cyk63MWlu+fNX0/xa3KoiUOTgfKF1ZXgvhKRhYABLgVetDUq6AHUnWNfBLiUmxSRicBEgF69etkcjvJasGvpePH+3pb5fnt9Mf83dzXgmhjaDZgEgFR15I3DZ7h09wB0bBPT6OsqFQxW5jlME5G3OfITf50xZo29YXkepl0vrlnALHBUZbU5HuWtYNfS8eL9mzsSZMvOA4x61P02XE1SAEdyubDnDfynQNxKbR8sr2LRmm36l74KKVaHsrYB9htjnheRriLSxxiz1ca4ioCedZ4nAdttfD9llybKRYTq+1sZCdJQqYut08eweOticld3d0suCz9exr6ySpf9Kw8bZizND7vkEOoz3ZVvrNRW+huQhmPU0vNADPAf4HQb41oB9BORPsA24HLgysYPUSEp2LV0bHh/Ywx9Ji92a9849XziWjuKFTeUXErqJYYaVuo2hZJwmOmufGPlyuFSYCiOiqwYY7aLSLvGD/GNMaZKRG4BluIYyvqcMeYbO99TNVNzhocGu5aOH9/f07KcH9x9Nsld4j3s7c7Xuk2hwpfCb3rFER6sJIcKY4wREQMgItZ+C3xkjFlM/YIzKjSEctVPm3hKCs9dm8a5Jxzj+YAGkqdjDsJ6l5LgcTHRZIXZHARvZ7rrFUf4sJIc5ovI00CCiNwATACetTcsFdJCveqnH3lKCreeexx3jW7kw7yR5FlTlTXc5yB4O9M9EKWmlX9YGa30sIicB+zHcd9hijHmHdsjU6Er2MNTmznj2RueksLgngm8frOFW21NJM9ImIPg7XDfsKitpQBrN6QfNMbcC7zjoU21RHEdoWyv53a72dyl5SkpQDMX2wl28gwAbwu/hUNtLeVgpVvpPKB+IrjAQ5tS9rOpS2v8U5/zVaF7wvNqreZgz+0IEG8Kv3l7xaECr7GqrP8H3AQcKyLr6mxqB3xmd2AqhJX90rx2f/LzX+W5725m5rvfubV7lRRqBHtuRwgLRKlp5R+NXTm8DLwNTAcm1Wk/YIzx0KegWoxg/mXsp/f+8LtdXPPcV27tPiWFGsGe2xHi7C41rfyjsaqsJUCJiOQCe+tUZW0nIqcaY74MVJAqxATzL2Mf3/unvYc486H33dr9khTqCvbcDqV8ZOWew5NAap3npR7aVEsSzL+MvXzv8spqTrh/iVt7wT/GEBWlKysrVZ+V5CDGmNpKYcaYwyKiy4u2dMH8y7gZ791QqYt12aNpH6vVUJVqiJUP+QIRuQ3H1QI4blIX2BeSUv7haVjq0tvPon83W6u/KBURrCSHPwOPAffhKJv9Hs51FJQKRZ6SwmNXDGXs4MQgRKNUeLIyQ3onjqqoSoU0T0nhmtN6M7RXRx58exOZ/10TtuUqlAq0xuY53GOMeUhEHqfeQjsAxpjbbI1MKYs8JYWjWkVRUXWY19du5+WvfqxdYGfbvjImL1gPoAlCqUY0duWw0fl1ZSACUZFr0ZptthSaO+fhD9i6u9StPS4murbqaf2FdQDKKqvDcnEdpQKpsXkObzq/vhC4cFSkWbRmm0uJan/85f7v97cwY2m+W3thTgan5yz3uF5CfeG2uI5SgdZYt9KbeOhOqmGMGWtLRCqizFia77J2AXj/l/u73/7M9S+6X8jWncC2fV8ZY6M+4Z5W80mU3Ww3XXioajxvHD7D5ZhwW1xHqUBrrFvpYefX3wLdcCwNCnAFUGhjTCqCNPQXenP+ct/88wHOm/mRW7unWc3XtP2KeyqfpY1UAJAku8mJeRYqqU0Q4bi4jlKB1li30ocAIvKAMeasOpveFBH331SlPEhMiOPnw59xVNelSMw+TGUCv+5K55ioEU0eu+9QBUOmui8d0lipi3tiXqFNVYVLWxupYFLr+bxZfoaOVlLKIivzHLqKSF9jTAGAiPQButoblooUo0/Zxqs/LECiHDeGpfU+YrsvYHTv3g0eU1V9mOP++rZb+/f/GEN0E6Uu2pR5XjQmkT1s9Xf9JKUimJXkcAfwgYjUzIpOBm60LSIVUT7d+1JtYqghUZV8uvcl4I9u+3salro+ezTtrJa6aCFrKShlNyuT4JaISD/gBGfTJmPMr/aGpSKF1WUhPSWF9+8+mz5d4pv3hrqWglJ+YWWZ0DbAnUBvY8wNItJPRPobY96yPzwV7ppaFtJTUnhxwimcdbyXPZe6loJSfiF1Cq563kHkFWAVcLUxZqCIxAGfG2OGBCA+S9LS0szKlTpXLxTlFeR5XBZy14Zst33vyxjA9Wf2DWB0SrVsIrLKGJPmaZuVew7HGmMuE5ErAIwxZSKiBfCVJfWXhawuvp5d+4512Wfs4EQeu2JoMMJTSjXASnKocF4tGAARORbQew7Ksoy+GWwpOJ6ZK13Xaj7u6La8e+dvghSVUqoxVpLD34AlQE8RmQucDlxrZ1AqcnxZsIfLZn3h0pbQJoa1U0YHKSKllBWNJgcRiQI64pglPRwQINMYszsAsakwVrDrIOc+8qFbe2FOBovWbOP0nOVeF+Kzq5CfUuqIRpODc0nQW4wx8wH3YSVK1bO3tILUBxqe1exrIT47CvkppdxZ6VZ6R0TuBl4BausjG2P22haVCjvlldWccP8St/b6pS58LcTnz0J+SqmGWUkOE5xfb67TZoCwH3OYV5BXO4qmW3w3MlMza0fXqMb9fflL/G/rM1RH/UJpfo7b9oZKXfhaiM8fhfyUUk2zMkO6TyACCbT64++LS4vJ/iwbQBNEE/6+/CVe/WEmB7fcB+Yol23fTk2nTeuGf6wSE+I8rrdgtYS2r8crpayJamoHEYkVkTtFZIGI/E9EbheR2EAEZ6fc1bkuE7MAyqvLyV2dG6SIwsd/vviRg/kPuCSG+OOm0b7fg40mBoCs9P7ExUS7tDWnhLavxyulrLHSrfQicAB43Pn8CuAlYJwdAYlINnADsMvZ9BdjzGJ/v4/Vmj/qiCMrsA2pbYs/7h9ExewH4HDjk+0Bx03jHj+9Rc/VMzja7GKndOWn1CyGDT3fUgw19xV0tJJS9rKSHPobYwbXef6+iHxtV0BOM40xDze9m/eaqvmjjnhtVRF3v+r6LW/T91Gij9rp0ibVHZt+sXXzGbb+b0AZCHRjF93W/w2SO1qufxTTYS3xx+XSrnQH8fHdiOmQCWhyUMqfmuxWAtaIyPCaJyJyKvCpfSEFRmZqJrHRrr1jsdGxZKZmBimi0PPhd7tInpTnkhjm33ga7busIyrmF5d9zeEY5JcLmn7R96a6VkwFx/P3plqKqeZeUXFpMQZTe68or0BHWivlT1auHE4FrhaRH53PewEbRWQ9YIwxKTbEdYuIXA2sBO4yxvxSfwcRmQhMBOjVq1ez3yCjbwb8+AW5BQvZEQXdDkNm7/P1ZjSwYVsJFz7+iUvbk39I5YJB3QE4sDuF6Ipqt9Xdqvdb+FEoKWpeez2N3SvS751S/mMlOVjrDG4GEXkXx7rU9f0VeBJ4AMdw2QeARzgynLaWMWYWMAscVVmbHcS6+WR8+gwZdf+KLX4GOg1qseWdf9xziLNmvO/Sln3RiVx7uuuANceIoaFU7XctltfDyoghHxfj0XtFSgWGlaGsP/j7TY0xo6zsJyLPAPasG9FY90YLSw6eZjXfeFZfJo8Z4HH/rPT+LrOUoRkjhnxcjEfvFSkVGFauHAJKRLobY2p++y8FNtjyRj52b0SCsopqBkxxndWcMag7//5DaqPH+TRiyMfFeDJTMz2uD6H3ipTyr5BLDsBDIjIER7dSIXatV92C1xquqj7M8fe97TL09KTE9rx16xlYXarjkqE9vB8+mjLe66uz+utD6Mx2pewRcsnBGOO+6rwdWuBaw8YYRj76IQW7aktkEd86mnXZ6R5LXYSqjL4ZmgyUslnIJYeAaWFrDV//wgre3eg6L2HTA+cTW2+2sVJKQUtODuBT90a4yH7jG+Z8VujS9vWU0XRoExOcgJRSYaFlJ4cI9sxHBUxbvNGl7dNJ51obbqqUavE0OUSY19duI3PeWpe2tzPPZED39sEJSCkVljQ5BEAglrX8bMturnz2S5e2l68/lRHHdfHr+yilWgZNDjaze1nLjcX7uSD3Y5e23MuHcPEQLUSnlPKeJgeb2bWs5bZ9ZZyes9yl7S9jTmDiWcd6/ZpKKVVDk4PN/L2sZcmhSgZPXebSdu2IZLLHnuTV6ymllCeaHGyWmBDHyfvf4Z5W80mU3Ww3XXioajyr2p/XrNcpr6zmhPtdS12MPOFoZl87zJ/hKqUUoMnBdv88cTMDVz1LnFQAkCS7eTDmWTacmAyc2+Tx1YcNKdlLKa040jV1bNd43r3zN5ZLXSilVHNpcrDZsO8fB2diqBEnFY72RspGGWPIeOwTvi3eX9sWHSXkP3A+raKtrNGklFLe0+RgNy+qv940dxWL17uuT7Bx6vnEtdZSF0qpwNDkYLdmVH+dvngjT39U4NK2+v7z6BTf2q7olFLKI00OdrNQ/XXOp1vJfvNbl8M+yjqHXp3bBCpKpZRyocnBbo1Uf317fTH/N3e1y+5v3nIGg5I6BCFQpZQ6QpNDINSr/rqicC/jJuW57DLnumGc3f/oQEemlFIeaXIIoM0/H+C8mR+5tD08bjC/PznyV59TSoUXTQ4BsO9QBUOmvuPSlpXen5vPOS5IESmlVOM0OdiorKKa3z/1Gd9sPzJX4YpTejH9t4OCGJVSSjVNk4MNKqsPc8OLK/kgf1dt2x2jjidzVL8gRqWUUta16OTg73UWjDHc+791zF95ZILbVcN78cDFA7XUhVIqrLTY5ODvdRYeXprPv97fUvs8/aRj+PeVqVrqQikVllpscvDXOgv1J7AN6ZnAvInDiY3RUhdKqfDVYpODr+ssvPn1dm7975ra5z0S4liceSYd4mLc9s0ryCN3dS47SnfQLb4bmamZZPTN8C5wpSzQnznlqxabHBIT4tjmIREkJsQ1etwnm3dz1ewjazXHRAuf3Hsux7SP9bh/XkEe2Z9lU15dDkBxaTHZn2UD6C+rsoX+zCl/aLEd4lnp/Ymr1/UTFxNNVnp/j/uvK9pH8qQ8l8Tw/t1ns3namAYTA0Du6tzaX9Ia5dXl5K7O9SF6pRqmP3PKH1rslUPNfYWmRisV7DrIuY986NL2xi2nk5KUYOl9dpTuaFa7Ur7SnznlDy02OYAjQTR083nn/nJOf3A5ldWmtu0/fzqVM/p1adZ7dIvvRnFpscd2peygP3PKH1pst1JDSsoqGTH9PU75x3u1ieHxK4ZSmJPR7MQAkJmaSWy0a7dTbHQsmamZfolXqfr0Z075Q4u+cqirvLKay2d9wdqf9tW2/X3sSVwzItmn1625AagjR1Sg6M+c8gcxxjS9V4hLS0szK1eu9OrYqurD/N/c1bzz7c+1bbeccxx3N3BjWimlIoWIrDLGpHna1qKvHMoqqhkwZUnt88vSepLzu0Fa6kIp1eIF5Z6DiIwTkW9E5LCIpNXbNllEtohIvoik2x1Lv6Pbck7/rmyedgEP/j5FE4NSShG8K4cNwG+Bp+s2isiJwOXASUAi8K6IHG+MqXZ/Cd/FtY7mnTt/Y8dLK6VUWAvKlYMxZqMxJt/DpouBecaYX40xW4EtwCmBjU4ppVSoDWXtAfxU53mRs82NiEwUkZUisnLXrl2edlFKKeUl27qVRORdwNOsm78aY15v6DAPbR6HUxljZgGzwDFayasglVJKeWRbcjDGjPLisCKgZ53nScB2/0Sk/EmrfioV2UKtW+kN4HIROUpE+gD9gK+CHJOqp6bqZ3FpMQZTW/UzryAv2KEppfwkWENZLxWRIuA0IE9ElgIYY74B5gPfAkuAm+0aqaS8p1U/lYp8QRnKaoxZCCxsYNs0YFpgI1LNoVU/lYp8odatpMJAQ9U9teqnUpFDk4NqNq36qVTka9G1lZR3tOqnUpFPk4PySkbfDE0GSkUw7VZSSinlRpODUkopN5oclFJKudHkoJRSyo0mB6WUUm40OSillHKjyUGpAMsryGP0a6NJeSGF0a+N1oKFKiTpPAelAqimom1N4cKairaAzhtRIUWvHJQKIK1oq8KFJgelAkgr2qpwoclBqQDSirYqXGhyUCqAtKKtChd6Q1qpANKKtipcaHJQKsC0oq0KB9qtpJRSyo0mB6WUUm40OSillHKjyUEppZQbTQ5KKaXciDEm2DH4TER2AT80sVsXYHcAwvFWKMensXknlGOD0I5PY/Nec+LrbYzp6mlDRCQHK0RkpTEmLdhxNCSU49PYvBPKsUFox6exec9f8Wm3klJKKTeaHJRSSrlpSclhVrADaEIox6exeSeUY4PQjk9j855f4msx9xyUUkpZ15KuHJRSSlmkyUEppZSbiEoOIjJORL4RkcMiklZv22QR2SIi+SKS3sDxnUTkHRHZ7Pza0cZYXxGRtc5/hSKytoH9CkVkvXO/lXbFU+89s0VkW534xjSw3/nO87lFRCYFKLYZIrJJRNaJyEIRSWhgv4Cdt6bOgzg85ty+TkRS7Yyn3nv3FJH3RWSj83fDbeEIETlbRErqfL+nBDC+Rr9PwTp3ItK/zvlYKyL7ReT2evsE9LyJyHMislNENtRps/SZ5dXvqjEmYv4BA4D+wAdAWp32E4GvgaOAPsD3QLSH4x8CJjkfTwIeDFDcjwBTGthWCHQJ8HnMBu5uYp9o53nsC7R2nt8TAxDbaKCV8/GDDX2PAnXerJwHYAzwNiDAcODLAH4vuwOpzsftgO88xHc28FYgf8asfp+Cee7qfY934JgwFrTzBpwFpAIb6rQ1+Znl7e9qRF05GGM2GmPyPWy6GJhnjPnVGLMV2AKc0sB+LzgfvwBcYkugdYiIAOOB/9r9Xn52CrDFGFNgjKkA5uE4f7YyxiwzxlQ5n34BJNn9nk2wch4uBl40Dl8ACSLSPRDBGWOKjTGrnY8PABuBHoF4bz8J2rmrYyTwvTGmqSoMtjLGfATsrdds5TPLq9/ViEoOjegB/FTneRGef0GOMcYUg+OXCjg6ALGdCfxsjNncwHYDLBORVSIyMQDx1LjFeRn/XAOXqlbPqZ0m4Pir0pNAnTcr5yEUzhUikgwMBb70sPk0EflaRN4WkZMCGFZT36dQOHeX0/Afb8E6bzWsfGZ5dQ7DbiU4EXkX8LQa+1+NMa83dJiHNtvH8FqM9Qoav2o43RizXUSOBt4RkU3OvyBsiw14EngAxzl6AEe314T6L+HhWL+cUyvnTUT+ClQBcxt4GVvOm6dwPbTVPw9B+flzCUCkLfA/4HZjzP56m1fj6DI56Ly/tAjoF6DQmvo+BfXciUhrYCww2cPmYJ635vDqHIZdcjDGjPLisCKgZ53nScB2D/v9LCLdjTHFzkvXnd7EWKOpWEWkFfBb4ORGXmO78+tOEVmI4xLR5w85q+dRRJ4B3vKwyeo5bTYL5+0a4EJgpHF2qnp4DVvOmwdWzoNt58oKEYnBkRjmGmMW1N9eN1kYYxaLyBMi0sUYY3txOQvfp6CeO+ACYLUx5uf6G4J53uqw8pnl1TlsKd1KbwCXi8hRItIHR3b/qoH9rnE+vgZo6ErEX0YBm4wxRZ42iki8iLSreYzjZuwGT/v6U70+3UsbeM8VQD8R6eP86+pyHOfP7tjOB+4FxhpjDjWwTyDPm5Xz8AZwtXPkzXCgpKYrwG7Oe1qzgY3GmEcb2Kebcz9E5BQcnwt7AhCble9T0M6dU4NX9sE6b/VY+czy7nc1UHfaA/EPxwdZEfAr8DOwtM62v+K4Y58PXFCn/VmcI5uAzsB7wGbn1042xzsH+HO9tkRgsfNxXxwjC74GvsHRrRKI8/gSsB5Y5/wh6l4/NufzMThGv3wfwNi24Og/Xev891Swz5un8wD8ueZ7i+Oy/t/O7eupM5IuAOfrDBxdCOvqnLMx9eK7xXmevsZxk39EgGLz+H0KoXPXBseHfYc6bUE7bziSVDFQ6fyc+1NDn1n++F3V8hlKKaXctJRuJaWUUs2gyUEppZQbTQ5KKaXcaHJQSinlRpODUkopN5ocVIslIgkiclOw46ghIteKyL+CHYdSoMlBtWwJgMfkICLRgQ1FqdCiyUG1ZDnAsc5a/DOc9fnfF5GXgfUiklyvdv7dIpLtfHysiCxxFoz7WEROqPvCIhIljrUKEuq0bRGRY0TkIhH5UkTWiMi7InJM/cBEZI6I/L7O84N1HmeJyApnYcS/+/OEKFVDk4NqySbhKMU8xBiT5Ww7BccM0hObOHYWcKsx5mTgbuCJuhuNMYdxlDK4FEBETgUKjaNGzyfAcGPMUBzlk++xGrCIjMZR/uUUYAhwsoicZfV4pawKu8J7StnsK+NY86NBzgqnI4BXnaV1wLGQVH2vAFOA53HUs3nF2Z4EvOKsYdUaaPT96hnt/LfG+bwtjmRhR1FB1YJpclDKVWmdx1W4Xl3HOr9GAfuMMUOaeK3PgeNEpCuORVj+n7P9ceBRY8wbInI2jpX36qt9b2dxt9bOdgGmG2Oebvq/opT3tFtJtWQHcCyd2ZCfgaNFpLOIHIWjTDjGUap5q4iMg9p1jgfXP9g4CpctBB7FURW1pmJnB2Cb8/E19Y9zKuRIKfeLgRjn46XABOfVCyLSw7kWglJ+pclBtVjOD+tPRWSDiMzwsL0SmIpj5bS3gE11Nv8B+JOI1FQUbWjZxVeAqzjSpQSOK4VXReRjoKHa/88AvxGRr4BTcV7RGGOWAS8Dn4vIeuA1Gk9wSnlFq7IqpZRyo1cOSiml3GhyUEop5UaTg1JKKTeaHJRSSrnR5KCUUsqNJgellFJuNDkopZRy8/8BgRrnjVGTm7cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# which fits data best?\n",
    "X,y = generate_data(30, w_randn)\n",
    "Xtest,ytest = generate_data(20,w_randn) # also try w_sparse and w_nonneg\n",
    "\n",
    "yhat_ridge = Ridge(alpha = alpha_ridge).fit(X,y).predict(Xtest)\n",
    "yhat_lasso = Lasso(alpha = alpha_lasso).fit(X,y).predict(Xtest)\n",
    "yhat_nonneg = LinearRegression(positive=True).fit(X,y).predict(Xtest)\n",
    "\n",
    "plt.plot(ytest,ytest)\n",
    "plt.scatter(ytest, yhat_ridge, label=\"ridge\")\n",
    "plt.scatter(ytest, yhat_lasso, label=\"lasso\")\n",
    "plt.scatter(ytest, yhat_nonneg, label=\"NNLS\")\n",
    "plt.xlabel(\"true value\")\n",
    "plt.ylabel(\"predicted value\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcceed6d670>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyFklEQVR4nO3deXiU1fXA8e9JSEjYkR0CBFwQC2ExIrIoihIlitjWpVixoqWtC7jBD1BppFWxuGFbFyruWKEUUYwCiituLIKgIi4QNSHI1rCERLKc3x8zGTLJJJkkM3lnOZ/n8THvnXfmPUa4Z9773nuuqCrGGGOiT4zTARhjjHGGJQBjjIlSlgCMMSZKWQIwxpgoZQnAGGOiVCOnA6iNtm3banJystNhGGNMWFm/fv0eVW1XsT2sEkBycjLr1q1zOgxjjAkrIvK9r3YbAjLGmChlCcAYY6KUJQBjjIlSYfUMwJeioiKys7MpLCx0OpQGlZCQQFJSEnFxcU6HYowJU2GfALKzs2nevDnJycmIiNPhNAhVZe/evWRnZ9OjRw+nwzHGhKmwHwIqLCykTZs2UdP5A4gIbdq0ibq7HmNMYIV9AgCiqvMvE43/zcaYwIqIBGCMMZFq2+5D/OOtbygqKQ34Z1sCCILRo0eTl5dXqT0jI4P77ruv4QMyxoQdVeUPz63jrPvf5b6VX/PTgcAP+Yb9Q+BQo6q8+uqrxMRYbjXG1M2GH/7HRY986Dl+6NL+JLVuEvDrRF0vtXRDDkNnv0WPaZkMnf0WSzfk1Pszs7Ky6N27N9deey0DBw4kNjaWPXv2AHDXXXfRq1cvzj77bLZu3ep5z9q1a0lJSeG0005jypQp9OnTB4CSkhKmTJnCKaecQkpKCo8//ni94zPGhIcjxaUkT8v0dP4dWjRm61/PZeyALkG5XlQlgKUbcpi+ZDM5eQUokJNXwPQlmwOSBLZu3cr48ePZsGED3bt3B2D9+vW8+OKLbNiwgSVLlrB27VrP+VdddRWPPfYYH330EbGxsZ72+fPn07JlS9auXcvatWv517/+xfbt2+sdnzEmtN2+dDMn3P665/i5qwfxyYyzadwotpp31U9UDQHNWbGVgqISr7aCohLmrNha7wzbvXt3Bg8e7NX2/vvvc9FFF9GkievWbcyYMQDk5eVx8OBBhgwZAsC4ceN49dVXAVi5ciWbNm1i8eLFAOzfv59vvvnG5vsbE6H2FxTR786VXm3f3nUejWKD//08qhLAjryCWrXXRtOmTX22+5quqapVfo6q8ve//520tLR6x2SMCW2/fvRD1n3/P8/xpaldWf3tHo6/7XU6t0pkSlov4lpuZO6nc9mZv5OOTTsyeeBk0numB+T6UTUE1LlVYq3a6+v000/npZdeoqCggIMHD7Js2TIAWrduTfPmzfn4448BePHFFz3vSUtL49FHH6WoqAiAr7/+mvz8/KDEZ0JX5rZMRi0eRcozKYxaPIrMbZlOh2QC6JufDpI8LdOr83/o0v688tkOryHqGSuf4Y7VfyY3PxdFyc3PJePDjID9eYiqO4Apab2YvmSz1zBQYlwsU9J6BeV6AwcO5NJLL6V///50796d4cOHe16bP38+v//972natCkjRoygZcuWAFxzzTVkZWUxcOBAVJV27dqxdOnSoMRnQlPmtkwyPsygsMQ17a/sLz0QsG9+xjnJ07w77+evPpVhx7dl6Oy3Kg1RyzGvU6Q/e7UVlhQy99O5AfmzINUNR4Sa1NRUrbghzJYtW+jdu7ffn7F0Qw5zVmxlR16B5xYrWE/Yq3Po0CGaNWsGwOzZs8nNzWXu3Lm1+oza/reb8DBq8Shy83MrtXdq2omVv17p4x0mHLz39W7GP7nGqy1r9tFOvMe0TCr2xs1PnAY+Fv0LsOnKzX5fW0TWq2pqxfaougMAGDugiyMdfkWZmZncc889FBcX0717d55++mmnQzIhYmf+zlq1m9BX8Vv/y9cNpV/XVl5tnVslklPheWT74lJ2xVUeqe9YEpgv7lH1DCCUXHrppWzcuJHPP/+czMxM2rWrtF2niVIdm3asVbsJXc9+lFWp88+anV6p8wfXEHVinPeUzxv3/Y+EUu8SEAmlpUzeuy8g8UXdHYAxoW7ywMlezwAAEmITmDxwsoNRmdooLVV6znjNq+2j6WfRqWXVE07KRibKD1GP1BbE7NnH3Nat2Nkolo7FJUz+Xx7pjdoEJE5LAMaEmLKHe8Ga+meC67aXNrPgkx88x706NGfFTad7nVPVs8hKQ9SbZpH28g2kZ+/wNBXHJsCFMwMSqyUAY0JQes906/DDTMGREnrPXO7V9sWdaTRt7N3NllUkKJvxU1aRAKj0fHJpyVBWF13DjbxIZ9nLDm3DQ6WXMaxkKGMDELOjCUBEWgFPAH0ABSao6kdOxmSMMbU15h+r2ZS933N88clJzLm4n89zq61IEPsBrJoF+7OhZRIb83/F4iNDWMwQr/M/CkD1AnD+DmAusFxVfy0i8UDgy901gGbNmnHo0CGnwzDGNLCfDhRy6t2rvNq23T2amJiqN2yqqvJA6oE3YNlTUOR+ff+PTNVH2BdzhFdKh/n1GbXlWAIQkRbA6cDvAFT1CHDEqXiMMaY2Ks7umXn+SUwYVnPNLl/TPQGmx//naOfv1kSOMLXRIl454p0AAlW9wMlpoD2B3cBTIrJBRJ4QkUoFdURkooisE5F1u3fvrv9VNy2CB/tARivXvzctqv9nuh06dIiRI0cycOBA+vbty8svvwxAfn4+6enp9OvXjz59+rBw4UIApk2bxkknnURKSgq33norAN9//z0jR44kJSWFkSNH8sMPP1R5PWNMw/s8Z7/PqZ3+dP7ge7pnYlwsHdjj8/zOsrfSuYGqXuDkEFAjYCBwg6p+IiJzgWnAHeVPUtV5wDxwrQSu1xU3LYJlk7xusVg2yfVzyiX1+miAhIQEXnrpJVq0aMGePXsYPHgwY8aMYfny5XTu3JnMTNcfmv3797Nv3z5eeuklvvrqK0TEs4PY9ddfz/jx47nyyit58sknmTRpkpWCMCZEVOz4n/xdKmed2KFWnzF2QBe6/PgqXT+dQ3vdzS5px48DpyDfJbn6pAoKm3SkS2JiUKoXOFYKQkQ6Ah+rarL7eDgwTVWrnPpQ71IQD/bx+QumZVe46XN/Q6+k7BlAUVERN910E++99x4xMTFs3bqV7du3c+DAAdLS0rjkkks4//zzGT58OMXFxZx88smkpqaSnp7O+eefT3x8PG3btiU3N5e4uDiKioro1KmTZ3OZiqwUhDENY/nnO/nj8+u92sqXcaiVil9EAeISod84+OyFyu0XPFzvL6hVlYJwbAhIVXcCP4pI2b3MSODLoF50f3bt2mtpwYIF7N69m/Xr17Nx40Y6dOhAYWEhJ5xwAuvXr6dv375Mnz6dWbNm0ahRI9asWcOvfvUrli5dyrnnnuvzM32VkzbGNJzkaZlenf/yG4fXvfMH1yyfCmP9FBXANytdnX3LroC4/h2Azr86Ts8CugFY4J4BtA24KqhXa+n7FouWSQH5+P3799O+fXvi4uJ4++23+f777wHYsWMHxxxzDL/97W9p1qwZTz/9NIcOHeLw4cOMHj2awYMHc9xxxwEwZMgQXnzxRa644goWLFjAsGHDqrukMSZIHn3nO+5d/pXnODZG+O7u0fX/4Oq+iKZcEtQOvyJHE4CqbgQq3ZYEzciZvm+9RgZmVd3ll1/OBRdcQGpqKv379+fEE08EYPPmzUyZMoWYmBji4uJ49NFHOXjwIBdeeCGFhYWoKg8++CAADz/8MBMmTGDOnDm0a9eOp556KiCxGWP8U1KqHFuhjMPa286mXfPGgblAkL+I1kbUlYNm0yKvhRaMnNmgGTeQ7BmAMYE16d8beOWzo2UXBiUfw6I/nhbYi1T1DCCIwz1WDrpMA99iGROqMrdlWr0ht4OFRfTN8N5r4au/nEtCXBA2ZC/rf0Lgi2j0JQBjjO06Vs6IOW+Ttfew53jC0B7MvOCk4F40RL6IWgIwJgrN/XSuV7lpCOxWg+Hgx32HGf63t73att8zOqpm3lkCMCYKRfuuYxUXdM3+ZV8uG9TNoWicYwnAmCjUsWlHn/sOR/quY+u//x+/evRDr7Z6zekPc5YAjIlC0bjrWMVv/S/8/lSGHNvWoWhCg+0JHAAiwi233OI5vu+++8jIyAAgIyODJk2asGvXLs/rzZo18/lzma1btzJixAj69+9P7969mThxYvCCN1EpvWc6GUMy6NS0E4LQqWknMoZkROT4/9INOT6Lt0V75w92BxAQjRs3ZsmSJUyfPp22bSv/oWrbti33338/9957r1+fN2nSJG666SYuvPBCwLWQzJhAi/Rdx1SVHtO9F3S9fesIerStVHS4wYXKFNyouwPI3JbJqMWjSHkmhVGLR5G5LbPmN9WgUaNGTJw40bOat6IJEyawcOFC9u3b59fn5ebmkpR0dFVg37596x2jMdHkvhVbvTr/Nk3jyZqdHjKdf8aHGeTm56KoZwpuIPqi2oqqBBDMX/x1113HggUL2L9/f6XXmjVrxoQJE5g7d65fn3XTTTdx1llncd555/Hggw96SkUbY6p3pLiU5GmZ/OPtbz1tn80cxfo7znEwKm/VTcFtaFGVAIL5i2/RogXjx4/n4Ycf9vn6pEmTeOaZZzhw4ECNn3XVVVexZcsWLr74Yt555x0GDx7Mzz//XO8YjYlkE55eywm3v+45Prt3B7Jmp9OySZyDUVUWSlNwoyoBBPsXf+ONNzJ//nzy8/MrvdaqVSvGjRvHI4884tdnde7cmQkTJvDyyy/TqFEjPv+87vsVGBPJ8g4fIXlaJm99dXSixTd3nccTVzZcncnaqGqqrRNTcKMqAQT7F3/MMcdwySWXMH/+fJ+v33zzzTz++OMUFxdX+znLly+nqKgIgJ07d7J37166dAnMDkDGRJL+s1bSf9YbnuNJI48na3Y6cbGh27VNHjiZhNgErzanpuCG7m8pCBriF3/LLbdUuYNX27Ztueiii7yGcw4fPkxSUpLnnwceeICVK1fSp08f+vXrR1paGnPmzKFjx8heoGNMbXy3+xDJ0zLJO1zkadt+z2huPucEB6PyTyhNwY26ctChMv0qEKwctIlGFef0z72sPxf2tzvk6lg5aLdIn/tsTKT68Ns9jHviE6+2aC7jEAhRlwBMZImkOzpTtYrf+pdcO4SB3Vo7FE3kiIgEoKpRVcIVXP/N0c5q2ke+Fz75gRkvea+Et2/9geN4AhCRWGAdkKOq59f2/QkJCezdu5c2bdpETRJQVfbu3UtCQkLNJ0cwq2kfuXyVcVj9f2eS1LqJQxFFJscTADAZ2AK0qMubk5KSyM7OZvfu3YGNKsQlJCR4lYuIRqG0oMYETsYrX/D0h1me42PbNWXVLSMciyeSOZoARCQJSAfuAm6uy2fExcXRo0ePgMZlwkO01rSPVIVFJZx4x3Kvts/vTKNZ41D4nhqZnF4H8BAwFSit6gQRmSgi60RkXbR9yzfVC6UFNaZ+fv3oh16d/y8HdCFrdrp1/kHm2G9XRM4HdqnqehEZUdV5qjoPmAeudQANE50JB2Xj/IGYBWSziZyx62Ahg+5a5dX23d2jiY2Jjud5TnMyvQ4FxojIaCABaCEiz6vqbx2MyYSZQKzrsNlEzqg4tXPG6BOZePqxDkUTnRwbAlLV6aqapKrJwGXAW9b5GyeEUnneaLAl94DPHbqs8294NsBmop7NJmo4FTv+J8ancvZJHRyKxoREAlDVd4B3HA7DRCmbTRR8q7b8xNXPeNfxsgVdlTX0s6iQSADGOGnywMlezwDAZhMFUsVv/ZmThvGLzi0diiZ0OfEsyulpoMY4LpTK80aSJ97f5nOs3zp/35x4FmV3AMZgVWIDqaRUOXaGdxmHNTNG0r5FdJcuqYkTz6IsARhTga0JqLtbFn3Gfz/N9hz379qKpdcNdTCi8OHEsyhLAMaUE01rApZuyGHOiq3syCugc6tEpqT1YuyAum2skv9zMb/48wqvtq/+ci4JcbGBCDUqOPEsyhKAMeVES4XRpRtymL5kMwVFJQDk5BUwfYmr7HJtk8CoB9/l658OeY6vPK07d17YJ3DBRolArmz3lyUAY8qJljUBc1Zs9XT+ZQqKSpizYqvfCWBHXgFDZr/l1bbt7tHEWBmHOmvoZ1GWAIwpJ1rWBOzIK6hVe0UVZ/f8ZWwfrhjcvd5xmYZl00CNKSdaKox2bpVYq/YyG3/M8zm10zr/8GR3AMaU48Q4rBOmpPXyegYAkBgXy5S0XlW+p2LH/9zVgxh+fLugxWiCzxKACYpwnkoZDWsCysb5/ZkF9OqmHVz/wgavNivjEBn8SgAi0h04XlXfFJFEoJGqHgxuaCZcRdNUynA2dkCXah/4+tqX982bz+C49s2CHZppIDU+AxCR3wOLgcfdTUnA0iDGZMKclVcOf3Pf/Mar82+e0Iis2enW+UcYf+4ArgMGAZ8AqOo3ItI+qFGZsBYtUykjUVFJKcff9rpX24Y7zqF103iHIjLB5E8C+FlVj4i45vaKSCPAtmY0VYqWqZSR5g/PrWPFFz95jk8/oR3PThjkYEQm2PxJAO+KyAwgUUTOAa4FlgU3LBPOrLxyeNlfUES/O1d6tX391/OIb2SzxCOdPwlgGnA1sBn4A/Aa8EQwgzLhLVqmUoa9V28mefVIr6ZrRxzL1HNPdCgg09BENXxGc1JTU3XdunU1n2iMqdbHC2Zx2eaTvdq2Nx6HnHI1nP+AQ1GZYBGR9aqaWrHdn1lA20VkW8V/ghOmMSbYkqdlenX+Y2I+ICthHCLA+qcdi8s0PH+GgMpnjQTgYuCY+l5YRLoCzwIdgVJgnqraPEFjgmT+6u385dUvvdqyEsZ5n6TeBeJMZKsxAajq3gpND4nIamBmPa9dDNyiqp+KSHNgvYi8oapf1vRGY0ztVCreFvc0V8SurHyiWP3+aFJjAhCRgeUOY3DdETSv74VVNRfIdf98UES2AF0ASwDGBEjFqZ3gLuPw6irw9Tjt5N81SFwmNPgzBHR/uZ+LgSzgkkAGISLJwADci80qvDYRmAjQrVu3QF7WmJBV3926fJVxWPzH00hNdo/elj3oXf+0a9hHYl2dvz0Ajir+DAGdGcwARKQZ8F/gRlU94OP684B54JoFFMxYjAkF9d2t67gZr1Fc6v1XxVfxtqVdbmHO5+cfTTJdejG2/uGbMFJlAhCRm6t7o6rW+6uCiMTh6vwXqOqS+n6eMZGgrrt1Hfq5mD4V9uX9ZMZIOrRIqHTu0g05zFj5DNLmdZp2zCOvqBUzVp4HXFnnfYFN+KnuDqDe4/zVEVdtifnAlkAkE2MiRV1266r4kBeqL9l817sLiGm/GIkpAkDi89D2i7nr3UaMHTC1lhGbcFVlAlDVO4N87aHAFcBmEdnobpuhqq9V/RZjIl/nVonkuDv7MTGrmdpoEZ1lD7ukHWzKh5Sjj+C278nnzPve8Xq/P2UcDjddRoy78y8jMUUcbroMsAQQLfyZBZSAqxTEL3CtAwBAVSfU58Kquhqw3aONqaBst65zSt5lRKvnmXBMM3Y26krH4hImvzmFdICUSyp96+/Rtilv3zrCr2vExOXVqt1EJn9mAT0HfAWkAbOAy4EtwQzKmGhWNgaf/+Zk7mnbgsIY17f53LhGZLRuxuaV/+GxF5p6vae2O3S1jG/P/qJdPttN9PCn3N9xqnoHkK+qzwDpQN/ghmVMdBsb+wHPtI7xdP5ldm/9G4/t+Z3n+MrTutdpe8bpg28mThp7tcVJY6YPrnbuh4kw/twBlA0U5olIH2AnkBy0iIyJdpsWwbJJ7OzS1tP08+6zObLnbK/T6rMvr1VsNeBfApgnIq2BO4BXgGbun40xwbBqFhQV0LG4hNy4RhzcMtvr5WYdlvD5TfPrfZn0nunW4Uc5fxLAU6paArwL9AxyPMaY/dkAZH93J4fVe6y/Za/biN8/zte7jKk1fxLAdhFZDiwE3tJw2kDAmDBU3KIbx+26x6utSbdH6BSfRe6uy7lz1OUORWYijT8JoBdwAa7N4Z8UkWXAi+5pnMaYAHJN7fTu/LMSxnH4p3j+Fnct/dNtpa4JHH9qARUAi4BF7mcBc3ENB1ndWGMCYdMidq18kEF7bvdqXtduFm0PboWWXWkyciYZKQGtwWiMX3cAiMgZwKXAecBaAlwN1JiotWkRyS80Bbw7/6xx+ZBSqTiuMQHlz0rg7cBGXHcBU1Q1P9hBGRMNPvpuL7+psKDr28a/pZGUwqquXiUfjAkGf+4A+vkq02yMqTufxdvKb8/onglkTDD58wzAOn9j6sDXpi45eQXMWbHV67xK+/ICtExqoChNNPPrGYAxpnaWbshh9UuPsJAX6dx4DzsOt2XYwoe9zhndtyOP9P0OliVCUblSz3GJMLK+W24bUzNLAMYEwcbMeQxp9iwTjmnGdzl3U1KQ7PX60TIOJ7v+tWqWa9inZZKr87fxf9MAHN0RzJhI1TP+39zTrgW7t/7Nq31Mk0U8PPMZ75NTLrEO3zjCnx3BegGn4KoDBK5FYe8FMyhjwt3MfQ/CPu+25r2n8VVRsTMBGeNDjTuCichKYKCqHnQfZwD/aZDojAkzBUdK6D1zuVdbk+6PENvkBwB2NrL1kyZ0+PMMoBtwpNzxEawctDGV+Jra2bz3NK/jjvGtGigaY2rmz4YwzwFrRCRDRP4MfAI8G4iLi8i5IrJVRL4VkWk1v8OY0PPtrkOVOv97x0O7k7yrpidIHJMHT2/I0Iyplj/rAO4SkdeB4e6mq1R1Q30vLCKxwD+Bc4BsYK2IvKKqX9b3s41pKD4XdLln+DRLsA1XTGjzdxpoE+CAqj4lIu1EpIeqbq/ntQcB36rqNgAReRG4ELAEYELe0g053Lhwo1fbtrtHExMjnmPbcMWEOn9qAf0ZSMU1G+gpIA54Hhhaz2t3AX4sd5wNnOrj+hOBiQDdunWr5yWNqb/qvvUbE078uQO4CBgAfAqgqjtEpHn1b/GL+GirtNmMqs4D5gGkpqbaZjQm+DYt8rkwa/KLG3h54w6vU63jN+HMnwRwRFVVRBRARJrW9AY/ZQNdyx0nATuqONeYhuHekN1TmmH/j7Bskrtk81GjTurAvPGpDgRoTOD4kwAWicjjQCsR+T0wAXgiANdeCxwvIj2AHOAywDY7Nc7wfOv/0as5ufAFKPQ+1b71m0jhzyyg+0TkHOAArucAM1X1jfpeWFWLReR6YAWu3cWeVNUv6vu5xtTapkWw9FooLfI0qUKPn1/wOu3eX/Xl0lMC/xzKV9VQ2/bRNAR/HgLfq6r/B7zho61eVPU14LX6fo4x9fL6/3l1/smFL1Q6JavD/8Epnwf80ks35DB9yWYKikoAyMkrYPqSzQCWBEzQ+bMQ7BwfbecFOhBjnKIFrqI9BzWxUue/In4qWc2vDlp55jkrtno6/zIFRSWV9gwwJhiqqwb6J+Ba4FgR2VTupebAh8EOzJiG5Otb//aEcUjLrjDy4aBV69yRV1CrdmMCqbohoBeA14F7gPJlGg6q6j7fbzF1lbkt01aNOiAnr4ChFTr/LY1/R6Ic4ZAm0OymwA/7lNe5VSI5Pjr7zq0Sg3pdY6CaISBV3a+qWcBcYJ+qfq+q3wNFIlJpwZapu8xtmWR8mEFufi6KkpufS8aHGWRuq7zgyARO8rRMhs5+y3PcT74jK2EcieKqfVgscUGPYUpaLxLjvCuEJsbFMiWtV9CvbYw/00AfBQaWO8730WbqYe6ncyks8Z5rWFhSyNxP54blXUBD3s3UZQbNR9/t5Tf/+tirbXvjcUiFpYktOBTocCspi9VmARkn+JMARFU9K3BVtVREwmYryXAYWtmZv7NW7aGs7G6mLKGV3c0AAf+912UGTcUyDn84oyfjPz6/UucPkKttaIhueOyALtbhG0f4Mwtom4hMEpE49z+TgW3BDiwQwmVopWPTjrVqD2VV3c3MWD2DlGdSGLV4VMB+/7WZQfPUB9srdf5Zs9OZfl5v7i26hMMa7/XaYY3n3iLbptFENn8SwB+BIbhW65YVbJsYzKACpbqhlVAyeeBkEmITvNoSYhOYPHCyQxHVXVV3LaVaGvAkvCOvgDExq1kdP4ltjcexOn4SY2JWV5pBkzwtkzuXHS0y+49xA7xW82bqcKYVXUN2aVtKVcgubcu0omt4pXQYQ2e/xdINOfWO1ZhQ5M9K4F24yjSEnXAZWikbGgn1oSp/dGzakdz83GrPqc3zjerG+K9stoapRU/QxP3QNkn2MDvuCY6JiwfSmfTvDbzyWc3F235zalee/3gYrxwZVuk1W5hlIll16wCmqurfROTv+K7SOSmokQVAVZ1RKA6tRErt+MkDJ3s9A6iKP0m4pjH+qXELaVJ8xOs9TeQItzZaRPK0QV7tr08eTu9OLXxe569j+wLw709+pEQrF5wtG1ayBGAiTXVDQFvc/14HrPfxT8iLpKGVcJHeM52MIRl0atoJQYgR33/E/EnCNY3xNymonET6Fz5On7z7vNqyZqdX2fmX+evYvnx3z2ifNcohuAuzMrdlMmrxqIA/IzGmJlXeAajqMve/n2m4cAIrkoZWwkn5u5mKs4LA/yRc4yrZlkme6p0HNZG+P8/3Om/DHefQuml8xbdXq6EXZjXkrCljKqpuCGgZPoZ+yqjqmKBEFGCRMrQSruqThGvsjEfOdNXqP+jd8QvK9tnn1yneKWm9vIadILgLsyJtDYgJL9U9BC67j/4l0BHXNpAAvwGyghiTiTB1TcI1dcY7ul3AkIPeG7V8e1k+jfrXffpmQy/MCpeJCiYyVTcE9C6AiPxFVU8v99IyEXkv6JGZqFddZ1xxTv+Zvdrx1FWDfH1Mna7bUA98w2migok8/qzobSciPVV1G4B7B692wQ3LGJeKnfH67//nc0FXuPI1a8omKpiG4k8CuAl4R0TKVv8mA38IWkTGlFdug/bkwgVeLz10af+wn5ppExWMk0R9zHuudJJIY+BE9+FXqvpzUKOqQmpqqq5bt86JSxsnuDdo/29hKrcU/cnrpXD+1m9MQxOR9aqaWrHdny0hmwA3A91V9fcicryI9FLVV4MRqDEeq2ZVmuGzNP4O+rc+AlgCMKa+/BkCegrXwq/T3MfZwH+AOicAEZkDXAAcAb4DrlLVvLp+ngltdSnZfFfml/zrp3u92rISxrl+2F/Vci1jTG34kwCOVdVLReQ3AKpaIOKreG6tvAFMV9ViEbkXmA7Ue5N5E3pqW7K5pFQ5dsZrXm0fNb6eTlJuE7qWScEL2Jgo4k810CMikoh7UZiIHAvU6xmAqq5U1WL34ceA/Y2OULUp2XzxYx96df5tGiv/bH8tV3ZNICW5K6OSOpPZolVQNmi3cgwmGvlzB/BnYDnQVUQWAEOB3wUwhgnAwqpeFJGJuMtPd+vWLYCXNQ2hrGTz1EaL6Cx72KFt+VvxJSzLO1p580BhESkZK73et2XWubyVvZyM1W0o1CIAcuMakdG2DTRrGtAnAFaOwUSrahOAiMQArXGtBh4MCDBZVffU9MEi8iauFcQV3aaqL7vPuQ0oBhb4OA8AVZ0HzAPXLKCarmtCS00lmyvO6U/7RQcev8I1WWHup3M9nX+ZQi0KeJkEK8dgolW1CcC9/eP1qroIqNU9saqeXd3rInIlcD4wUv2Zi2rCUlUlmy+NebNSyebt94ym/OOlhiqTYOUYTLTyZwjoDRG5FdcwTX5Zo6ruq/ot1RORc3E99D1DVQ/X9XNMiCq3eKuJj3qCyYUvQLkv3FPSenHdmcdVOq+hyiRYOQYTrfx5CDwBuA54j6N7AdR3NdY/gOa4kstGEXmsnp9nQoV78ZarTLN35/9ByS9cnX85WbPTfXb+0HD7Odi+ESZa+bMlZI9AX1RVff+NN+Fv1SwoqlzCuWLHP+/0QkaN/lW1H9VQZRKsHIOJVjWWghCRBOBaYBiur3TvA4+pavV7/gWBlYIIfZrRCin3zf/Z4nOYWXyV1zlZ4/Ihpe4lm40xtVPnUhDAs8BB4O/u498AzwEXBy48Eyl+oi0d2Y0q9PjZ+1v/6/HT6D1rs0ORGWMq8icB9FLVfuWO3xaRz4IVkAlv9xy5mFNivuL24mu82rMSxkHLrg5FZYzxxZ8EsEFEBqvqxwAicirwQXDDMuGouKSUl0uH8XLp0UVemxpfTQspoIB4EoOwgtcYU3f+JIBTgfEi8oP7uBuwRUQ2A6qqKUGLzoSNv776JU+s3u457kkOz8bfSzMKyS5ty/2ll3JGyVDGOheiMaYCfxLAuUGPwoStg4VF9K1QxqFlQiO2FXZh2JGHvdrXrNga9hu4GBNJ/JkG+n1DBGLCzxXzP+H9b45WBfnL2D5cMbg7Pab5XjS+I6/y9NBQkbkt06aBmqjjzx2AMV5y8goYOvstr7byZRw6t0okx0dn37lVYoPEV1tWDM5EK39WAhvj0X/WSq/O/6mrTiFrdrpXDZ8pab1IjIv1el9iXCxT0no1WJy1UV0xOGMimd0BGL9sys5jzD+8J39VtS9v2Th/bXcBc4oVgzPRyhKAqVHFks2vTx5O704tqn3P2AFdQrbDr8iKwZloZUNApkorv9jp1fl3aZVI1uz0Gjv/cGPF4Ey0sjsAU4mq0mO69768a2aMpH2LhCreEd6sGJyJVpYAjJcn3t/GXzO3eI5HndSBeeMr1ZCKOOk9063DN1HHEoAB4EhxKSfc/rpX25ez0mgSb39EjIlU9rfbMH3JJv695kfP8fVnHsetITpl0xgTOJYAolje4SP0n/WGV9t3d48mNkaqeIcxJpJYAohSY//5ARt/zPMc33dxP359cpJzARljGpwlgCjz/d58zpjzjldbVQu6jDGRzdEEICK3AnOAdqq6p6bzTf30nJ5JabkdQF+cOJjBPds4F5AxxlGOJQAR6QqcA/xQ07mmftZl7ePXj33k1Wbf+o0xTt4BPAhMBV52MIaIV7GMw6pbzuDYds0cisYYE0ocSQAiMgbIUdXPyleRrOLcicBEgG7dujVAdJHh5Y05TH5xo+f4xI7NWX7j6c4FZIwJOUFLACLyJuCrmtZtwAxglD+fo6rzgHkAqampWsPpUa+0VOk5w7uMw6d3nMMxTeMdisgYE6qClgBU9Wxf7SLSF+gBlH37TwI+FZFBqmr1d+vh4VXf8MAbX3uOfzmwCw9c0t+5gIwxIa3Bh4BUdTPQvuxYRLKAVJsFVHeFRSWceMdyr7av/nIuCRU2ZTHGmPJsHUCYu+HfG1j22Q7P8dRze3HtiOMcjMgYEy4cTwCqmux0DOFo98GfOeWuN73ayu/La4wxNXE8ARi3TYtg1SzYnw0tk2DkTEi5xOepI+9/h+9253uO/zluIOkpnRoqUmNMhLAEEAo2LYJlk6CowHW8/0fXMXglgW9+Osg5D77n9VZb0GWMqStLAKFg1ayjnX+ZogJXuzsBVFzQ9dK1QxjQrXVDRWiMiUCWAELB/uwq2z/4dg+XP/GJpykxLpYtfzm3gQIzxkQySwChoGWSa9inguTCBVCu839/6pl0PaZJQ0ZmjIlgMU4HYGDtsTdQoEdX6i4sHkFy4Que40E9jiFrdrp1/saYgLI7gBBw45fHc3LRNdwS+x/OKJrr9dpnfx5Fy8Q4hyIzxkQySwAhYEdeATkM45XSYV7tAtb5G2OCxoaAHHakuJSq1m51bpXYsMEYY6KKJQAHLd2Qwwm3v+61S1eZxLhYpqT1avigjDFRw4aAHHCwsIi+GSs9x+f16ciokzpw38qv2ZFXQOdWiUxJ68XYAV0cjNIYE+ksATSwf723jbte2+I5fuuWM+jp3qHrooFJToVljIlClgAayK4DhQy6e5XneMLQHsy84CQHIzLGRDtLAA3gzmVf8NQHWZ7jNbeNpH3zBOcCMsYYLAEE1bbdhzjr/nc9x7en9+aa4T0djMgYY46yBBAEqsofn1/Pii9+8rRtzhhF8wSb02+MCR2WAAJsww//46JHPvQcz72sPxf2t9k8xpjQYwkgQEpKlQv+vpovcw8A0KllAu9OOZP4RrbUwhgTmiwBBMBbX/3EhKfXeY6fv/pUhh3f1sGIjDGmZo4lABG5AbgeKAYyVXWqU7HUVWFRCYPuepMDhcUAnJLcmoUTTyMmxvblNcaEPkcSgIicCVwIpKjqzyLS3ok46mPR2h+Z+t9NnuNXbxhGny4tHYzIGGNqx6k7gD8Bs1X1ZwBV3eVQHLW2/3AR/WYdLeMwtn9nHrpsgIMRGWNM3TiVAE4AhovIXUAhcKuqrvV1oohMBCYCdOvWreEi9OHvq77h/je+9hy/N+VMurWxTVqMMeEpaAlARN4EOvp46Tb3dVsDg4FTgEUi0lNVK9XFVNV5wDyA1NRUH3Uzgy93fwGn3fOW5/i6M49lStqJToRijDEBE7QEoKpnV/WaiPwJWOLu8NeISCnQFtgdrHjqavqSzfx7zQ+e4/W3n02bZo0djMgYYwLDqSGgpcBZwDsicgIQD+xxKBafvv7pIKMefM9zPOvCXzD+tGTnAjLGmABzKgE8CTwpIp8DR4ArfQ3/BMSmRbBqFuzPhpZJMHImpFxS5emqyvgn1/D+N658FB8bw8Y/n0OTeFsyYYyJLI70aqp6BPht0C+0aREsmwRFBa7j/T+6jsFnElizfR+XPP6R5/jRywdyXt9OQQ/TGGOcENlfa1fNOtr5lykqcLWXSwDFJaWMeug9tu3OB6Bn26asuOl04mKtjIMxJnJFdgLYn11j+/LPc/nj8596jhf94TQG9Tgm2JEZY4zjIjsBtExyDfv4aD98pJj+s97gSHEpAMOPb8uzEwYhYmUcjDHRIbLHOEbOhLhE77a4RJ7rNouTZq7wdP4rbjyd564+1Tp/Y0xUiew7gLJxfvcsoH3NT2Dg7j+De83xbwZ15Z5fpjgXnzHGOCiyEwCwtGQoc35+mJzCAlfRCbcPp51F51aJVb/RGGMiXEQngKUbcpi+ZDMFRSWetkYxwn0X97PO3xgT9SL6GcCcFVu9On+A4lJlzoqtDkVkjDGhI6ITwI68glq1G2NMNInoBFDVMI8N/xhjTIQngClpvUiMi/VqS4yLZUpaL4ciMsaY0BHRD4HHDugCuJ4F7MgroHOrRKak9fK0G2NMNIvoBACuJGAdvjHGVBbRQ0DGGGOqZgnAGGOilCUAY4yJUpYAjDEmSlkCMMaYKCXB2oo3GERkN/C903GU05YQ28y+CuESJ1iswWKxBke4xNpdVdtVbAyrBBBqRGSdqqY6HUdNwiVOsFiDxWINjnCK1RcbAjLGmChlCcAYY6KUJYD6med0AH4KlzjBYg0WizU4winWSuwZgDHGRCm7AzDGmChlCcAYY6KUJYB6EpEbRGSriHwhIn9zOp6aiMitIqIi0tbpWKoiInNE5CsR2SQiL4lIK6djqkhEznX/f/9WRKY5HU9VRKSriLwtIlvcf0YnOx1TTUQkVkQ2iMirTsdSHRFpJSKL3X9Wt4jIaU7HVFuWAOpBRM4ELgRSVPUXwH0Oh1QtEekKnAP84HQsNXgD6KOqKcDXwHSH4/EiIrHAP4HzgJOA34jISc5GVaVi4BZV7Q0MBq4L4VjLTAa2OB2EH+YCy1X1RKAf4RGzF0sA9fMnYLaq/gygqrscjqcmDwJTgZB+8q+qK1W12H34MZDkZDw+DAK+VdVtqnoEeBHXF4GQo6q5qvqp++eDuDqpkN0gQ0SSgHTgCadjqY6ItABOB+YDqOoRVc1zNKg6sARQPycAw0XkExF5V0ROcTqgqojIGCBHVT9zOpZamgC87nQQFXQBfix3nE0Id6plRCQZGAB84nAo1XkI15eUUofjqElPYDfwlHu46gkRaep0ULUV8TuC1ZeIvAl09PHSbbh+f61x3VqfAiwSkZ7q0NzaGmKdAYxq2IiqVl2sqvqy+5zbcA1hLGjI2PwgPtpC+q5KRJoB/wVuVNUDTsfji4icD+xS1fUiMsLhcGrSCBgI3KCqn4jIXGAacIezYdWOJYAaqOrZVb0mIn8Clrg7/DUiUoqrONTuhoqvvKpiFZG+QA/gMxEB15DKpyIySFV3NmCIHtX9XgFE5ErgfGCkUwm1GtlA13LHScAOh2KpkYjE4er8F6jqEqfjqcZQYIyIjAYSgBYi8ryq/tbhuHzJBrJVtexuajGuBBBWbAiofpYCZwGIyAlAPCFYGVBVN6tqe1VNVtVkXH94BzrV+ddERM4F/g8Yo6qHnY7Hh7XA8SLSQ0TigcuAVxyOySdxZfz5wBZVfcDpeKqjqtNVNcn9Z/Qy4K0Q7fxx/935UUR6uZtGAl86GFKd2B1A/TwJPCkinwNHgCtD8NtqOPoH0Bh4w33H8rGq/tHZkI5S1WIRuR5YAcQCT6rqFw6HVZWhwBXAZhHZ6G6boaqvORdSxLgBWOD+ErANuMrheGrNSkEYY0yUsiEgY4yJUpYAjDEmSlkCMMaYKGUJwBhjopQlAGOMiVKWAEzEc1dtvNbpOMqIyO9E5B9Ox2GMJQATDVoBPhOAu7KnMVHJEoCJBrOBY0Vko3uvgRHuGvkv4FoglexezAd49kzIcP98rIgsF5H1IvK+iJxY/oNFJEZEssrvWeDeI6CDiFzgLhS4QUTeFJEOFQMTkadF5Nfljg+V+3mKiKx174twZyB/IcaAJQATHaYB36lqf1Wd4m4bhKvwXE218efhKvh1MnAr8Ej5F1W1FHgZuAhARE4FslT1J2A1MFhVB+AqGT3V34BFZBRwvDvO/sDJInK6v+83xh9WCsJEqzWqur26E9wVNIcA/3GXpABXiYqKFgIzgadw1bBZ6G5PAhaKSCdcdaKqvV4Fo9z/bHAfN8OVEN6rxWcYUy1LACZa5Zf7uRjvu+EE979jgDxV7V/DZ30EHCci7YCxwF/d7X8HHlDVV9zljTN8vNdzbXfhtnh3uwD3qOrjNf+nGFM3NgRkosFBoHk1r/8EtBeRNiLSGFcZatx187eLyMXg6qBFpF/FN7sLAL4EPICr6uZe90stgRz3z1dWce0s4GT3zxcCce6fVwAT3HchiEgXEWlf03+oMbVhCcBEPHeH/IGIfC4ic3y8XgTMwrVT1qvAV+Vevhy4WkQ+A76g6q0fFwK/5ejwD7i+8f9HRN6n6jLh/wLOEJE1wKm470xUdSXwAvCRiGzGVW++uiRmTK1ZNVBjjIlSdgdgjDFRyhKAMcZEKUsAxhgTpSwBGGNMlLIEYIwxUcoSgDHGRClLAMYYE6X+Hxw8M7k5WxmNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# which fits data best?\n",
    "X,y = generate_data(30, w_sparse)\n",
    "Xtest,ytest = generate_data(20,w_sparse) # also try w_sparse and w_nonneg\n",
    "\n",
    "yhat_ridge = Ridge(alpha = alpha_ridge).fit(X,y).predict(Xtest)\n",
    "yhat_lasso = Lasso(alpha = alpha_lasso).fit(X,y).predict(Xtest)\n",
    "yhat_nonneg = LinearRegression(positive=True).fit(X,y).predict(Xtest)\n",
    "\n",
    "plt.plot(ytest,ytest)\n",
    "plt.scatter(ytest, yhat_ridge, label=\"ridge\")\n",
    "plt.scatter(ytest, yhat_lasso, label=\"lasso\")\n",
    "plt.scatter(ytest, yhat_nonneg, label=\"NNLS\")\n",
    "plt.xlabel(\"true value\")\n",
    "plt.ylabel(\"predicted value\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcceeea99a0>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAw9klEQVR4nO3deXhU5dn48e+dBbJBwh5ChAAioKwRFXFDUaLEBWoVtUUrKm9/aknVYgNViGBLFBSjttQFxfpKgfoCLrGs4grKXqBlkULQhAQIkBBC9jy/P2YyyYQkTGBmziRzf66LK/M8Z5l7gDz3nOeccx8xxqCUUsr/BFgdgFJKKWtoAlBKKT+lCUAppfyUJgCllPJTmgCUUspPBVkdQGO0b9/exMXFWR2GUko1KZs3b841xnSo3d+kEkBcXBybNm2yOgyllGpSRORgXf06BaSUUn5KE4BSSvkpTQBKKeWnmtQ5gLqUlZWRmZlJcXGx1aF4VUhICLGxsQQHB1sdilKqiWryCSAzM5NWrVoRFxeHiFgdjlcYYzh27BiZmZl0797d6nCUUk1Uk58CKi4upl27dn4z+AOICO3atfO7ox6l/FH6/nRGfjiSAe8NYOSHI0nfn+62fTf5IwDArwb/Kv74mZXyN+n700lZl0Jxhe3LXnZhNinrUgBI7JF43vtv8kcASinVXKVtSeP06QhKcm/AGNtwXVxRTNqWNLfsXxOAB4waNYq8vLwz+lNSUpg9e7b3A1JKNTnGGPbvuYHC/b+j9OhITHlrx7Kcwhy3vEezmALyJcYYPv30UwICNLcqpc7Njsx8bnv9G6A/ACExCwkIznMsjw6Pdsv7+N0otWxrFlelfk735HSuSv2cZVuzznufGRkZ9O3bl0cffZT4+HgCAwPJzc0F4I9//CO9e/fmxhtvZM+ePY5tNm7cyIABA7jyyiuZNGkS/fr1A6CiooJJkyZx2WWXMWDAAN54443zjk8p1TRUVhrunLvOPvhDq1Bof/EMgiO3OdYJCQwhKT7JLe/nVwlg2dYsJi/ZQVZeEQbIyiti8pIdbkkCe/bs4f7772fr1q1069YNgM2bN7Nw4UK2bt3KkiVL2Lhxo2P9Bx98kL/+9a+sX7+ewMBAR/+8efOIjIxk48aNbNy4kbfeeosDBw6cd3xKKd/27b5cekz5jM0HTwAw/8HL2DEtkeeufpbO4Z0RhM7hnUkZluKWE8DgZ1NAs1bsoaiswqmvqKyCWSv2MHpwl/Pad7du3Rg6dKhT39dff82YMWMICwsD4PbbbwcgLy+PgoIChg0bBsB9993Hp59+CsDKlSvZvn07H374IQD5+fn88MMPer2/Us1UWUUlw2d9QVZeEQCXxLTm48evJjDAdqVfYo9Etw34tflVAjhk/wt2tb8xwsPD6+yv63JNY0y9+zHG8Nprr5GQkHDeMSmlfFv69mweW7DF0V7y6DDiu7bx2vv71RRQTFRoo/rP17XXXsvSpUspKiqioKCATz75BIA2bdrQqlUrvvvuOwAWLlzo2CYhIYG5c+dSVlYGwN69eyksLPRIfEopa5wuLafXHz5zDP439OnIgZmj6hz8PXHesopfHQFMSujN5CU7nKaBQoMDmZTQ2yPvFx8fz9ixYxk0aBDdunXjmmuucSybN28ejzzyCOHh4QwfPpzIyEgAHn74YTIyMoiPj8cYQ4cOHVi2bJlH4lNKed/73x3k2WU7He1VT1xLr06t6ly36rxl1ZhVdd4SOO9pawBpaDrC1wwZMsTUfiDMrl276Nu3r8v7WLY1i1kr9nAor4iYqFAmJfR2y19kY506dYqIiAgAUlNTyc7OJi2tcTd3NPazK6Wsc6KwlMEzVjna915+ATN/NqDBba5K/dxxbqCmLlGhfJt8g8vvLSKbjTFDavd7/AhARN4BbgWOGGP62ftSgEeAo/bVphhjPvN0LGDLmlYM+LWlp6czc+ZMysvL6datG/Pnz7c6JKWUh7yyei+vrP7B0f42+Qa6uDD17MnzluCdKaD5wOvA32r1zzHG+O1tsWPHjmXs2LFWh6GU8qBDeUUMS/3c0Z44ohdP3nSRy9vHRIXWeQTgrvOWHj8JbIz5Cjju6fdRSilf8syyHU6D/5Znb2rU4A+285ahwYFOfe48b2nlSeDHReR+YBPwlDHmRF0ricgEYAJA165dvRieUko13r4jBdz48leO9nO3X8IDw+LOaV9V09WeOm9pVQKYC8wAjP3nS8D4ulY0xrwJvAm2k8DeClAppRrDGMOE9zez6j+HARCBnSkJhLc8v2HWk+ctLUkAxpjDVa9F5C3gUyviUEopd9j2Ux6j//yto/3avYO5bWCMhRG5xpIEICKdjTHZ9uYYYGdD6/u6iIgITp06ZXUYSikvq6g0jP7zt+zIygcgJjKELyZdT4sgN55e3b4Y1kyH/EyIjIURU2HA3W7ZtTcuA/07MBxoLyKZwDRguIgMwjYFlAH8j6fjUEopd/py71EeeGeDo/3+Q5dzTa8O7n2T7Yvhk4lQZr8SKP8nWxvckgS8cRXQvcaYzsaYYGNMrDFmnjFmnDGmvzFmgDHm9hpHA563fTHM6QcpUbaf2xe7bdenTp1ixIgRxMfH079/fz766CMACgsLSUxMZODAgfTr149FixYBkJyczMUXX8yAAQP43e9+B8DBgwcZMWIEAwYMYMSIEfz4449ui08pdf5Kyyu5/I+rHYP/oAui2P+nUe4f/MH2zb+s1mWgZUW2fjfwq1IQns6mISEhLF26lNatW5Obm8vQoUO5/fbbWb58OTExMaSn2x7mnJ+fz/Hjx1m6dCm7d+9GRBxPEHv88ce5//77eeCBB3jnnXeYOHGiloJQykd8/K9DTPz7Vkd72WNXMeiCKM+9YX5m4/obya+KwXk6mxpjmDJlCgMGDODGG28kKyuLw4cP079/f1avXs3vf/97vv76ayIjI2ndujUhISE8/PDDLFmyxFEyev369dx3330AjBs3jm+++cYtsSmlzl1hSTlxyemOwT/hkk4cmDnKs4M/2Ob8G9PfSP6VADycTT/44AOOHj3K5s2b2bZtG506daK4uJiLLrqIzZs3079/fyZPnsz06dMJCgpiw4YN3HnnnSxbtoybb765zn3WVU5aKeU98789wCXTVjjaq5+8jjfGDfHO7+aIqRBc667f4FBbvxv41xRQZKxt2qeufjfIz8+nY8eOBAcHs3btWg4ePAjAoUOHaNu2Lb/85S+JiIhg/vz5nDp1itOnTzNq1CiGDh3KhRdeCMCwYcNYuHAh48aN44MPPuDqq692S2xKqcY5dqqES59f7WiPG9qNGaP7eTeIqqnppnoVkE8ZMdX5HAC4NZv+4he/4LbbbmPIkCEMGjSIPn36ALBjxw4mTZpEQEAAwcHBzJ07l4KCAu644w6Ki4sxxjBnzhwAXn31VcaPH8+sWbPo0KED7777rltiU0q5bvaKPby+dp+jvX7yDXSO9MxzQ85qwN1uG/Br87ty0J68ptbbtBy0Uu6VeeI0V7+w1tF+8qaLmDiil4URuYdl5aB9jgezqVKq6fr9h9tZtKl6injb1JuICmthYUSe538JQCmlath7uICRc6qLt/1xTD9+cUU3CyPyHk0ASim/ZIzhV+9u5Mu9tudStQgKYNvUmwhr4T/Dov98UqWUstt88Dh3zl3vaM/9RTy39O9sYUTW0ASglPIbFZWGxFe/ZndOAQBd24ax5qnrCA70r1uiqmgCUEr5hbW7j/Dg/I2O9oKHr2DYhe0tjMh6/pn23ExEeOqppxzt2bNnk5KSAkBKSgphYWEcOXLEsTwiIqLO11X27NnD8OHDGTRoEH379mXChAmeC16pZq6kvIL4Gascg/9lcW3Y/6dRfj/4gyYAt2jZsiVLliwhNze3zuXt27fnpZdecnl/EydO5IknnmDbtm3s2rWL3/zmN+4KVSm/smRLJr2fWc7xwlIAPnn8av7x62EEBGiJFfDDBJC+P52RH45kwHsDGPnhSNL3p5/3PoOCgpgwYYLjbt7axo8fz6JFizh+/LhL+8vOziY2tro8Rf/+/c87RqX8SUFxGXHJ6Ty5+F8A3DqgMwdmjqJ/bKTFkfkWv0oA6fvTSVmXQnZhNgZDdmE2KetS3JIEHnvsMT744APy8/PPWBYREcH48eNJS0tzaV9PPPEEN9xwA7fccgtz5sxxlIpWSp3d21/vp3/KSkd77e+G8/p98VpYsQ5+lQDStqRRXFHs1FdcUUzaFtcG5oa0bt2a+++/n1dffbXO5RMnTuS9997j5MmTZ93Xgw8+yK5du7jrrrv44osvGDp0KCUlJecdo1LN2dGCEuKS03k+fRcAD14VR0ZqIt3bh1scme/yqwSQU5jTqP7G+u1vf8u8efMoLCw8Y1lUVBT33Xcff/nLX1zaV0xMDOPHj+ejjz4iKCiInTub9GOTlfKomf/cxWV/rK7cuWHKCKbddomFETUNfpUAosOjG9XfWG3btuXuu+9m3rx5dS5/8skneeONNygvL29wP8uXL6esrAyAnJwcjh07RpcuXdwSo1LNyU/HTxOXnM4bX+4HYFJCbzJSE+nYOsTiyJoGv0oASfFJhAQ6/8cICQwhKT7Jbe/x1FNPNXg10JgxY5ymc06fPk1sbKzjz8svv8zKlSvp168fAwcOJCEhgVmzZhEd7Z4kpVRz8eSibVzzYnXlzn9NG8lj119oYURNj9+Vg07fn07aljRyCnOIDo8mKT6JxB6J7g7VK7QctPJHu7JPckva1472C3f2Z+xlXS2MyPdpOWi7xB6JTXbAV8qfGWP45bzv+XbfMQAiWgax6ZkbCQkOtDiypsvvEoBSqunZcOA4d79RXbztjXGXknCJTouer2aRAIwxfneNb1OaulPqXJVXVHJz2tfsO3IKgB4dwln522sJ8tPibe7W5BNASEgIx44do127dn6TBIwxHDt2jJAQvdJBNV+r/nOYR/5Wfc5v4YShDO3RzsKImp8mnwBiY2PJzMzk6NGjVofiVSEhIU7lIpRqLorLKrjs+dUUlNgul76yRzsWPHKF33zB86YmnwCCg4Pp3r271WEopdxg8aafePrD7Y72ZxOv4eKY1hZG1Lw1+QSglGr68ovKGPhcdf2e0YNieOWewRZG5B80ASilLPXXL/9L6j93O9pfTbqeru3CLIzIf2gCUEpZ4sjJYi7/0xpHe8K1PZgySm9s9CZNAEopr3v+0//w9jcHHO0NfxhBx1Z6VZu3aQJQSnlNRm4hw2d/4WhPGdWHCdf2tC4gP6cJQCnlFb/5+1Y++dchR3t7ykhahwRbGJHyeAIQkXeAW4Ejxph+9r62wCIgDsgA7jbGnPB0LEop79uZlc+tr33jaM++ayA/v1TvYfEF3rifej5wc62+ZGCNMaYXsMbeVko1I5WVhrvfWO8Y/KPCgtk942Yd/H2Ix48AjDFfiUhcre47gOH21+8BXwC/93QsSinvWP/fY9z71neO9rwHhjCibycLI1J1seocQCdjTDaAMSZbRDrWt6KITAAmAHTtqjW/lfJlZRWV3Pjylxw8dhqA3p1akT7xai3e5qN8/iSwMeZN4E2wPRDG4nCUUvVYvjObX//vFkf7H7++ksvi2loYkTobqxLAYRHpbP/23xk4YlEcSqnzVFRaweAZKykuqwTgml7t+dv4y7V4WxNgVQL4GHgASLX//MiiOJRS5+HvG35k8pIdjvby315Dn2gt3tZUeOMy0L9jO+HbXkQygWnYBv7FIvIQ8CNwl6fjUEq5T/7pMgZOry7e9vNLY5l910ALI1LnwhtXAd1bz6IRnn5vpZT7vf75D8xeudfR/vrp67mgrRZva4p8/iSwUso35OQXM3RmdfG2R4f35Omb+1gYkTpfmgCUUmeV8vG/mb8uw9He9MyNtI9oaV1Ayi1cSgAi0g3oZYxZLSKhQJAxpsCzoSmlrLb/6ClueOlLR3vqrRcz/mp9Al9zcdYEICKPYLsRqy3QE4gF/orO4SvVbBljePSDLfxzZ46jb+dzCUS01EmD5sSVf83HgMuB7wGMMT80dOeuUqpp256Zx+2vf+tovzJ2EKMHd7EwIuUpriSAEmNMadVNHSISBOgduUo1M5WVhjv/uo6tP+YB0D6iJd8mX0/LoEBrA1Me40oC+FJEpgChInIT8CjwiWfDUkp50zc/5PLLed872vMfvIzhvfVAv7lzJQEkAw8BO4D/AT4D3vZkUEop7ygtr+T62V+QlVcEQL8urfnosasJDNAyDv7grAnAGFMJvGX/o5RqJtK3Z/PYguribUseHUZ81zYWRqS8zZWrgA5Qx5y/MaaHRyJSSnnU6dJy+qespKLS9ms9ok9H3n5giBZv80OuTAENqfE6BFvdHq3xqlQT9P76DJ796N+O9qonrqVXp1YWRqSs5MoU0LFaXa+IyDfAVM+EpJRytxOFpQyescrRvvfyrsz8WX8LI1K+wJUpoPgazQBsRwT6lUGpJmLOqr2krfnB0f42+Qa6RIVaGJHyFa5MAb1U43U5kAHc7ZFolFJucyiviGGpnzvaE0f04smbLrIwIuVrXJkCut4bgSil3OcPS3fwwfc/Otpbnr2JtuEtLIxI+aJ6E4CIPNnQhsaYl90fjlLqfOw7UsCNL3/laE+/4xLuvzLOuoCUT2voCEDn+ZVqIowxPPK3TazeZXu8dmCAsH3aSMK1eJtqQL3/O4wxz3kzEKXUudn64wnG/GWdo/3avYO5bWCMhRGppsKVq4BCsJWCuATbfQAAGGPGezAupdRZVFQa7vjzN+zMOglATGQIX0y6nhZBARZHppoKV44P3wd2AwnAdOAXwC5PBqWUatiXe4/ywDsbHO33H7qca3p1sDAi1RS5kgAuNMbcJSJ3GGPeE5EFwApPB6aUOlNJeQVXv7CWowUlAAy6IIol/28YAVq8TZ0DVxJAmf1nnoj0A3KAOI9FpJSq00fbskhauK26/dhVDLwgyrJ4VNPnSgJ4U0TaAM8CHwMR9tdKKS84VVJOv2nVB903XxLN3F/Ga/E2dd5cSQDvGmMqgC8BrQCqlAct25rFrBV7OJRXRExUKJd3b8vSrVmO5Wueuo6eHSIsjFA1J64kgAMishxYBHxujNHHQSrlAcu2ZjFl5XtIu38S1r6MrH3POgb/+6/sxvQ7+lkcoWpuXLlerDewGtvD4TNE5HURudqzYSnlf/745QcEdPyQ0twRFO6rnmVt22GPDv7KI1ypBVQELAYW288FpGGbDtInRSvVWNsXw5rpkJ8JkbEwYioMsNVWLAjYRvGeGY5VW7RfScsOn1NSGgU0WJlFqXPi0n3iInIdMBa4BdiIVgNVqvG2L4ZPJkKZ7fm75P9kawNxC8Kx3WJjE97reQKCTgEQEJzn5UCVv3D1kZDbsB0FTDLGFHo6KKWapTXTqwd/u0+KB/KbBeGOdkBIJuHdX3daJ7JFR6+Ep/yPK0cAA40xJz0eiVLNXX4m6eFhpLWJIjswkII9LzgtfuF+eGHTW5TVuMwiWFoyeahO/yjPcOUcgA7+StnVvkxzUkJvRg/u4tK26R1iSQkz5OeOovTYcEf/kLCtfDj1GQAiQiBtSxo5hTlEh0eTFJ9EYo9ET3wUpVw7B6CUsg3+k5fsoKisAoCsvCImL9kB4FISmBPZjqM7n3bqi+j9DCdbhgO2BJDYI1EHfOU1lpYNFJEMEdkhIttEZJOVsSh1NrNW7HEM/lWKyiqYtWLPWbf95dvfs6/G4B8UuYVWfZORgHKyy2ocZG9fDHP6QUqU7ef2xe4KX6kz+MITwa43xuS6aV9KecyhvKJG9QOcKCxl8IxVTn0RfaYgUuloS3mU7UUDVwlVXSqqlDu58kSw3sBl2OoAAdwGfFXnFko1YzFRoWTVMdjHRIXWuf6lM1ZxrLDU0Q5oeYiwuLlOg7+pDKb48Ehbo46rhCgrsvVrAlAeUO8UkDHmOftTwdoD8caYp4wxTwGXArFuen8DrBSRzSIywU37VMojJiX05uct1vFNi4nsb3kf37SYyM9brGNSQm+n9Q4eKyQuOd1p8D8wcxSdQ3tSnP0zKkujMAYqS6Mozv4ZHQOG2VbKz6z7jevrV+o8uXISuCtQWqNdivvKQV9ljDkkIh2BVSKy2xjjdHRhTwwTALp27eqmt1XKdVVX/gw5uYrUFm8Tav91iJVcUgPfJihwIFX3RsYlpztt+8rYQY4TxJMSejN5SSmF/x3sWB4aHMikn9kTSGSsbdqntkh3fd9SypmrTwTbICJLsX1jHwP8zR1vbow5ZP95xL7/y6k1vWSMeRN4E2DIkCFaiE55Vc0rfxa1WOwY/KsEVRTDmulsjryJO+euc1qWkep8NU9VIqj3MtIRU53PAQAEh9r6lfIAcaW4p4jEA9fYm18ZY7ae9xuLhAMBxpgC++tVwHRjzPL6thkyZIjZtEkvFlLec1Xq5455//0t76OuB2/FFS9wai+cMJShPdqd2xs2UCtIqXMlIpuNMUNq97t6H0AYcNIY866IdBCR7saYA+cZUydgqf2hFkHAgoYGf6WsUPMKn0OmPf+KOE1amyhyggIJz7uE7Jz7ndav/a2/0QbcrQO+8hpXagFNA4ZguxroXSAY+F/gqvN5Y2PMfmDg+exDKU+reeXPk6HDONB+C8UBARTsSqXmLfKrn7yOCzvqg1pU0+LKjWBjgNuBQnDM27dqcAulmolJCb0Ja/Mvwnum8p/O2ziRcw8Fu1KrV5AyLhqSpoO/apJcmQIqNcYYETHgmLtXyi8ER24jpPMSSitLObU71WlZeK8ZBAQVklOoz+ZVTZMrCWCxiLwBRInII8B44G3PhqWUb0jbkkbevkepLO3k1N+qb7LjdXR4tLfDUsotXKkGOltEbgJOYjsPMNUYs+osmynV5J0qKWfvpiSnvoiLpiKB1ZeChgSGkBSfVHtTpZoEV04Cv2CM+T22yzRr9ynlM86nVHNttW/okhZHiej5UnWHMXSOiNFyzapJc2UK6Cag9mB/Sx19SlnmfEs1V8nKK+Kq1M+d+tpf9HtKAqvn+UMqK0nJPUHir3a6IXKlrNNQNdD/BzwK9BSR7TUWtQLW1b2VUtZoqFSzqwmg9rf+0YNieCV7HOnHjjmu/Y8uryDpRB6JQed4o5dSPqShI4AFwD+BmUByjf4CY8xxj0alVCOdS6nmKlt/PMGYv9RTxmH7VBI/mUhi5qHqhcGhkKDlGVTTV28CMMbkA/kikgYcN8YUAIhIKxG5whjzvbeCVOpsGluquUrtb/3Jt/Th19f1rO6ouitXyzOoZsiVcwBzgfga7cI6+pSylK3S5g6naaDQ4MAzSjVX+XT7IR5f4FzSqt4yDlqeQTVTriQAMTUqxhljKkVEnyWsfMpZK23WUPtb/xvjLiXhEr2WX/kfVwby/SIyEdu3frCdGN7vuZCUOjejB3c5Y8BP359O2pY0cgpzCM6/nWOHrnRaft7F25RqwlxJAL8GXgWewfY8gDXYH9CilC9L359OyroUisqLzyjj8OlvrqZfl0iLIlPKN7hyJ/AR4B4vxKKUW6VtSePEwTspL+jv1H/RkDT6ddFv/ko1dB/A08aYF0XkNWzf/J0YYyZ6NDKlGin9i2dJ27+UnADoVBHADz/MdFoefuGfCAg+SbYWb1MKaPgIYJf9pz6CS/m89C+eJeXAUooDhVM/JHOyPMppec3ibVJrmVL+qqH7AD6x/3zPe+Eo1ThV3/qzA8CYUE7tSnFaHnHRM0hguaMdUllJjyPdvBylUr6poSmgT6hj6qeKMeZ2j0SklItqfut3ekgLEBCSSXj3121F28oqnMo4XFqy3qKIlfItDU0Bzbb//BkQje0xkAD3AhkejEkpl6TtX8rpinYU/vdpp/6IPpOxP7+IzuUVrKxZxgEwnL08hFL+oKEpoC8BRGSGMebaGos+EZGvPB6ZUmexd+8LTu3gqO8J6bzU0Q6prCTpRN4Z20lkrKdDU6pJcOWZwB1EpEdVQ0S6Ax08F5JSDdtw4PgZd/O26ptcPfgbQ+cKQ0qrASSW1prFDA611fJRSrl0I9gTwBciUnX3bxzwPx6LSKkG1B74wzt+REC76jn9kMpKxkUkMPGul20d2xdrITel6uHKjWDLRaQX0MfetdsYU+LZsJS/q1nCITo8mqERSfztC+d1ukSF0qWkiGNllRwNEjqUG9odjWeR3IrjJhUt5KZUvVx5JGQY8CTQzRjziIj0EpHexphPPR+e8kdVJRyKK4oB2Lspib01lr/7q8u4vk9Huienk8U9kG+7Uf0k8F9A9CSvUi5xZQroXWAzUFVFKxP4B6AJQHlE2pY0iiuKKTmSQOmx652W1Szedq7PAFBK2bhyErinMeZFoAzAGFME6L30ymOyT+VQsCvVafAP6z6H1n0nO603KaE3ocGBTn0NPQNAKeXMlSOAUhEJxX5TmIj0BPQcgPKIX727gYI9zjV8qso4RId3dupvzDMAlFJnciUBTAOWAxeIyAfAVcCvPBmU8j+l5ZVc9Mw/nfrCez1PQNApAExlMFe1HXfGdnU9A0Ap5ZoGE4CIBABtsN0NPBTb1E+SMSbXC7EpPzFo+kryTpc59QW13gqVQRgDpiyKkqMJrDzchWk3WBSkUs1QgwnA/vjHx40xi4H0htZVqrFOFJYyeMYqp77dM26m77PLKT85mPKTg52WHdKre5RyK1emgFaJyO+ARdgeCA+AMea4x6JSzV7tG7ouj2vL4l/bLjTTq3uU8g5XEsB4+8/HavQZoEcd6yp1huc+f5//O/AWlYEnMEU9KTz4iNPyAzNHIVJ9YdmkhN58s/Qv/JaFxEguh0x7XuEerk541NuhK9WsuXIncHdvBKKap+c+f59/HJyDBJVxqlbJ5gevimPabZecsc3owG+5Nfhtguw3gsVKLqmBbxMUOBDQu3qVcpez3gcgIiEi8qSILBGR/xOR34pIiDveXERuFpE9IrJPRJLPvoVqav7vwFtUFHU9o15/614v1Dn4A7BmumPwrxJUUWyr6aOUchtXpoD+BhQAr9nb9wLvA3edzxuLSCDwZ+AmbHcXbxSRj40x/zmf/SrfcvKH3zu1W0YvoUWbDVTW+6ghbIXbGtOvlDonriSA3saYgTXaa0XkX25478uBfcaY/QAishC4A9AE0ATVLt42JDSJBV87r1PzubwBFW3q31lkLOT/VHe/UsptXCkFsVVEhlY1ROQK4Fs3vHcXoOZveaa9z4mITBCRTSKy6ejRo254W+VuVcXbsguzMRj2bnIe/ENi33Ua/E1lMHd2f6SOPdmNmGqr21+T1vFXyu1cSQBXAOtEJENEMoD1wHUiskNEtp/He9dVT+iMiQFjzJvGmCHGmCEdOuhzaHxRVfG24sO3njHXn5GayKOdutKxrBIxho5llUyIGM60G868q9dhwN1w26sQeQEgtp+3vaplnZVyM1emgG720HtnAhfUaMcCh+pZV/mw7FM5FOx2HvjDe8wmsOUx2F7IxL0LmFhW47r+wwtg+9CGB3St46+Ux7lyGehBD733RqCX/RGTWcA9wH0eei/lIXe/sZ6CA3UXb2sd3NF25U5ZrZu6yops/TrAK2UpV44APMIYUy4ijwMrgEDgHWPMv62KR9WhgccpFpdV0OfZ5U6rR/SajgSdBmzz/CVHEiD/T3XvW6/oUcpyliUAAGPMZ8BnVsag6rF9MXwysfrbe/5PtjbQa1EEZRXOp2uCWm/FVLYAc9pRvK3w5CXQSa/oUcpXWZoAlA+rY+omtzSYIQvCqXmufu/zt3D97C/IyjuzeFuXKPuVOzUTCegVPUr5CE0Aqm75maSHh5HWJoqcoEBO7n7BafHw3h2Y/+DlwFlq9wyw12+uZypJKWUdTQD+ooH5/Lqkd4glJcxQWNqZ0/uecFpWu3jbWWv36BU9SvkkTQD+oIH5/PoG5lmtW3P0339w6mvRbi0x0esQSXReuaHaPTrwK+WzXLkRTDV1DV2KWYe1e46QUWvwb9U3mZYdV3CsouDMDbR2j1JNkh4B+AGTn1n3bdd19Nd+UEtIzCKCI7c62pVlUWfuSGv3KNUk6RGAHzhM+7P2z//2wBmDf0TvZ50Gf1MZTFjhbWfuSGv3KNUkaQLwAzNL7+K0aeHUd9q0YGapraJ3XHI6KZ9UF2H9x6+v5JWxg6g88nMqS6MwBipLo6g88nP+cN0vznwDrd2jVJOkU0B+YFPrm0g+CU8HLSZGjnHItOPF8rtZHXTdGd/6M1JtJ3gviwN4gFkrhnIor4iYqFAmJfRm9OAzCrba6JU+SjU5mgD8wKSE3kxZ+R1r2sYgwWFUlkZxev/VUFrhWOerSdfTtV2Y03ajB3epf8BXSjV5mgD8wM6SdwmMXgRA4YHHqSx2Pjlb9a1fKeVfNAE0c+n701m0ZxGmMphTe2Y4Les5aC5r7vnUosiUUlbTBNDMpW1J4/TBCVSc7uHok6A8InqlcrSkrotDlVL+QhNAE7BsaxazVuxx7WRsDScKS9m7KcmpL6LPFEQqAYgOj/ZIvEqppkETgI9btjWLyUt2UFRmO2GblVfE5CU7ABpMAvEzVnG8sNTRbtFhBS3br3VaJyk+qfZmSik/ovcB+LhZK/Y4Bv8qRWUVzFqxp871M3ILiUtOdxr8O1z87BmD/9jeY0nsoSd/lfJnegTg4w7lFbncX/ua/rR7BnFHwLekrz5GWuswcoICiS6vIOnkaRLDLvJIvEqppkMTgI+LiQrl0pOr7Ddx2Wrtv1h+N5tb3+RYZ/PB49w5d73Tdo5LO+dMJ/FkHokn85x3rJU6lfJ7mgB83CsX/0C/zW8TKrYpnVjJ5YXgt9l5cRxwwxnf+hdNGMoVPdpVd2ilTqVUPTQB+LjL/vsaSKlTX6iUcuTfXxC3zoUburRSp1KqHpoAfF2tRzNGl1ewd99sqPH8lTVPXUfPDhF1bz9iKuUf/cbpgS3lgSEEaaVOpfyeXgXk49I7xJLSvi3ZwUGUHL/ONvjbhbUIJCM1sf7BH1hWcRXJZQ+TWdmeSiNkVrYnuexhllVc5Y3wlVI+TI8AfEkdz+1NaxNFUelJTu2a6bRqj4tf4fP7V511l7NW7CGrdBgfMsypf/2KPVroTSk/pwnAV9R4bm96eBhprSrI2TKdopyfUZZ3hWO1wNAMwuL+Sq5xrYxDYy4jVUr5F00AvmLNdNJbCKmdu5AXEIAxLTi127l4W0TvZ5CAcsD1Mg4xUaFk1THYx0SF1rG2Usqf6DkAH/Fp+XFS2rclLzCQgh+mOlXuDI76jlZ9kx2Df0hgiMtlHCYl9CY0ONCpLzQ4kEkJvd0XvFKqSdIjAB8xp00bikw4p3ZNc+qP6DMZEQOAIESHR5MUn+RyGYeqef5zKSanlGreNAH4iANZj1JZFOdoB7f9kpBO/3S0pbwN2x/66pz2rU/2UkrVRROAxX46fpprXlwLxDn6WvVNdlrHVAZTfHikdwNTSjV7mgAsNCBlBSeLyx3tkNh3CW5VXeXTGDAVoZQcvp1OAcPq2oVSSp0zPQlsgZ1Z+cQlpzsN/hmpidx78R1Q3gYMdCyrZOaRXFYcPMbo4iI9aauUcjs9AvCy2sXb/pl0DX07twZg2g3jmLa9pXPphoDTpAa+TVDgQECrdyql3EePALxk7Z4jToN/58gQMlITHYO/w5rpTnV7AFt7zXRvhKmU8iOWHAGISArwCHDU3jXFGPOZFbF4mjGG7pOdP9r6yTfQObKeG7G0fLNSykusnAKaY4yZffbVmq6FG34k2f78XoBrerXn/YeuaGALtHyzUspr9ByAB1RUGnpOcf7Wvz1lJK1Dgs++8YipjppADsGhtn6llHIjK88BPC4i20XkHRFpU99KIjJBRDaJyKajR4/Wt5p3bF8Mc/pBSpTt5/bFZ6zy8so9ToP/uKHdyEhNdG3wB9tjGm97FSIvAMT287ZX9fGNSim3E2OMZ3Ysshqoq2LZH4DvgFzAADOAzsaY8Wfb55AhQ8ymTZvcGqfLalTrdAgOdQzORaUV9J263GmTvc/fQosgPc+ulLKWiGw2xgyp3e+xKSBjzI2urCcibwGfeioOt1kz3XnwB1t7zXSS/tOLj7YdcnRPGdWHCdf29HKASinVOFZdBdTZGJNtb44BdloRR6PUcRXOcdOK+MMvwOHqwf/AzFGIuFarXymlrGTVSeAXRWQQtimgDOB/LIrDdZGxpJcfczybt+TA45SUXOBY/Nq9g7ltYIyFASqlVONYkgCMMeOseN/zkT54DCkHlnK6vB2Fu592WpaR6lppZqWU8iV6GaiL0nK/5+je58FUX80T2u0NLuhQDGgCUEo1PZoAXLDtpzz2bnJ+AldVyeacQp3vV0o1TZoAzqJ28bawHi8T2PKIo906uIO3Q1JKKbfwuwSwbGuWS49HXP2fwzz8t+p7DgICSwi78HkkoMzRZyqDKTmS4JW4lVLK3fwqASzbmsU3S//CIhYS0zKXQ6fb88rSe4BHHUmgruJtG6aM4Io/raE4+2e07LACCc7DlEVRcjSBwpOXWPBJlFLq/DX7BFDzG/8dgd/yp6C3CJNSAGIll+nmTV5MD2L04Od4/7uDPLus+paEEX06Mu9XlwEQExVKVt5gyk8Odtp/l6h6qnoqpZSPa9YJoPY3/koCCJJKp3XCpJRflSwgLvlyp/6dzyUQ0bL6r2dSQm8mL9lBUVmFoy80OFCf1KWUarKadQLYlv4m0+VNxzf+ACpJDw9z3MwVXV5B18wbWV12s2Obh67uzrO3XnzGvqqmiFw5f6CUUk1Bs04AD5f+L2EBpY52engYKe3bUhwQgKkMZu++Geytsf4Pf7yF4MD6i7eNHtxFB3ylVLPRrBNATMAxp3ZamyiKAwIoPXEFJTljHP0dYlaxceIrXo5OKaWs1awTQHFoNGsD8h1TPpWVLTi16zlqPgYhok8yJVq8TSnlh5p1Alh7xV2kHFhKcYBQenwoJYdHO5aF9XiJwJa2B8xEh9f12AKllGremnUCSMv9nuIAoShrrOPyzeCo7wnpvNSxTkhgCEnxSfXtQimlmq1mnQByCnMACAz9ifKTgwm/cCYBwfkACEJ0eDRJ8Ukk9tBibkop/9OsE0B0eDTZhdm0aLuOFm3XOfo7h3dm5c9XWhiZUkpZr1k/sDYpPomQwBCnPp3yUUopm2Z9BFA1tZO2JY2cwhyd8lFKqRqadQIAWxLQAV8ppc7UrKeAlFJK1U8TgFJK+SlNAEop5ac0ASillJ/SBKCUUn5KjDFWx+AyETkKHLQ6jkZoD+RaHYSF/P3zg/4d6Of3jc/fzRjToXZnk0oATY2IbDLGDLE6Dqv4++cH/TvQz+/bn1+ngJRSyk9pAlBKKT+lCcCz3rQ6AIv5++cH/TvQz+/D9ByAUkr5KT0CUEopP6UJQCml/JQmAA8TkRQRyRKRbfY/o6yOyRtE5GYR2SMi+0Qk2ep4vE1EMkRkh/3ffJPV8XiDiLwjIkdEZGeNvrYiskpEfrD/bGNljJ5Uz+f36d9/TQDeMccYM8j+5zOrg/E0EQkE/gzcAlwM3CsiF1sblSWut/+b++x14G42H7i5Vl8ysMYY0wtYY283V/M58/ODD//+awJQnnA5sM8Ys98YUwosBO6wOCblYcaYr4DjtbrvAN6zv34PGO3NmLypns/v0zQBeMfjIrLdfojYbA+Ba+gC/FSjnWnv8ycGWCkim0VkgtXBWKiTMSYbwP6zo8XxWMFnf/81AbiBiKwWkZ11/LkDmAv0BAYB2cBLVsbqJVJHn79db3yVMSYe2zTYYyJyrdUBKUv49O9/s38kpDcYY250ZT0ReQv41MPh+IJM4IIa7VjgkEWxWMIYc8j+84iILMU2LfaVtVFZ4rCIdDbGZItIZ+CI1QF5kzHmcNVrX/z91yMAD7P/p68yBthZ37rNyEagl4h0F5EWwD3AxxbH5DUiEi4irapeAyPxj3/3unwMPGB//QDwkYWxeJ2v//7rEYDnvSgig7BNgWQA/2NpNF5gjCkXkceBFUAg8I4x5t8Wh+VNnYClIgK237EFxpjl1obkeSLyd2A40F5EMoFpQCqwWEQeAn4E7rIuQs+q5/MP9+Xffy0FoZRSfkqngJRSyk9pAlBKKT+lCUAppfyUJgCllPJTmgCUUspPaQJQzZ6IRInIo1bHUUVEfiUir1sdh1KaAJQ/iALqTAD2yqVK+SVNAMofpAI97fXYZ4nIcBFZKyILgB0iElerhvvvRCTF/rqniCy3F3X7WkT61NyxiATYa/9H1ejbJyKdROQ2EfleRLba60V1qh2YiMwXkZ/XaJ+q8XqSiGy0FxJ7zp1/IUqBJgDlH5KB/9rrsU+y910O/MEYc7bnFLwJ/MYYcynwO+AvNRcaYyqxlTcYAyAiVwAZ9how3wBDjTGDsZXEftrVgEVkJNDLHucg4FItKKfcTUtBKH+1wRhzoKEVRCQCGAb8w17WAaBlHasuAqYC72Kre7TI3h8LLLLXg2kBNPh+tYy0/9lqb0dgSwj+WFBOeYgmAOWvCmu8Lsf5aDjE/jMAyDPGDDrLvtYDF4pIB2wPPHne3v8a8LIx5mMRGQ6k1LGt473FlmVa2PsFmGmMeePsH0Wpc6NTQMofFACtGlh+GOgoIu1EpCVwK4Ax5iRwQETuAtsALSIDa29sbAW1lgIvA7uMMcfsiyKBLPvrB2pvZ5cBXGp/fQcQbH+9AhhvPwpBRLqIiD8+TEV5kCYA1ezZB+Rv7Q/pmVXH8jJgOvA9tnrtu2ss/gXwkIj8C/g39T/achHwS6qnf8D2jf8fIvI1kFvPdm8B14nIBuAK7EcmxpiVwAJgvYjsAD6k4SSmVKNpNVCllPJTegSglFJ+ShOAUkr5KU0ASinlpzQBKKWUn9IEoJRSfkoTgFJK+SlNAEop5af+P3xUnkOwusGGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# which fits data best?\n",
    "X,y = generate_data(30, w_nonneg)\n",
    "Xtest,ytest = generate_data(20,w_nonneg) # also try w_sparse and w_nonneg\n",
    "\n",
    "yhat_ridge = Ridge(alpha = alpha_ridge).fit(X,y).predict(Xtest)\n",
    "yhat_lasso = Lasso(alpha = alpha_lasso).fit(X,y).predict(Xtest)\n",
    "yhat_nonneg = LinearRegression(positive=True).fit(X,y).predict(Xtest)\n",
    "\n",
    "plt.plot(ytest,ytest)\n",
    "plt.scatter(ytest, yhat_ridge, label=\"ridge\")\n",
    "plt.scatter(ytest, yhat_lasso, label=\"lasso\")\n",
    "plt.scatter(ytest, yhat_nonneg, label=\"NNLS\")\n",
    "plt.xlabel(\"true value\")\n",
    "plt.ylabel(\"predicted value\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcceefc7b80>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0cElEQVR4nO3dd5wUVbbA8d/pYXDIGZYMKogJEQdEVERBQFRQnwEjImFVXN19qytsMOyuK66uruHByppQCWJAgoIkUQyEISigJCUNIFGQDDN93h+3eqaBYaZhuru6p8/38xm66vatrlMM1Olbt+4tUVWMMcaYwgT8DsAYY0zis2RhjDGmSJYsjDHGFMmShTHGmCJZsjDGGFOkUn4HECvVq1fXRo0a+R2GMcYklXnz5m1V1RpHlpfYZNGoUSOysrL8DsMYY5KKiKwpqNwuQxljjCmSJQtjjDFFsmRhjDGmSCW2z6Ighw4dIjs7m/379/sdStxkZGRQr1490tPT/Q7FGJPEUipZZGdnU6FCBRo1aoSI+B1OzKkq27ZtIzs7m8aNG/sdjjEmiaXUZaj9+/dTrVq1lEgUACJCtWrVUqolZYyJjZRKFkDKJIqQVDteY0xspFyyMMaYEmvrSpj6OMTg0ROWLHzWtWtXduzYcVT5Y489xjPPPBP/gIwxyWndXHj1cpj/JvyyPuofn1Id3IlGVZkwYQKBgOVsY0wxLJsE794JFWrBbR9ApXpR30XMzlIi8pqIbBaRxWFlVUVkiois8F6rhL03UERWisgyEekcVn6eiCzy3ntBkvwi/OrVqzn99NO59957admyJWlpaWzduhWAJ554gtNOO42OHTuybNmyvG3mzp1L8+bNueCCC3jooYc466yzAMjNzeWhhx6iVatWNG/enJdfftmXYzLG+GjBcBh1C9RsBr2nQLVTYrKbWLYs3gBeAt4MKxsATFPVQSIywFt/WETOAHoAZwJ1gKki0lRVc4EhQD9gFvAx0AWYWNzgHh+/hO82/FLcjznMGXUq8ujVZxZZb9myZbz++usMHjyY0GSH8+bNY9SoUSxYsICcnBxatmzJeeedB0CvXr0YOnQobdu2ZcCAAXmf8+qrr1KpUiXmzp3LgQMHuPDCC+nUqZPdJmtMqpj3Box/AE6+FG56G04qH7NdxaxloaqfA9uPKO4ODPOWhwHXhJWPUtUDqroKWAm0FpHaQEVV/Vrdw8LfDNsmaTVs2JA2bdocVjZz5kyuvfZaypYtS8WKFenWrRsAO3bsYNeuXbRt2xaAW265JW+byZMn8+abb9KiRQvOP/98tm3bxooVK+J3IMYY/8x91SWKUy+Hm0fFNFFA/PssaqnqRgBV3SgiNb3yuriWQ0i2V3bIWz6yvEAi0g/XCqFBgwaFBhJJCyBWypUrV2B5QVfYtJC7GlSVF198kc6dOx+zjjGmBJo9FCY+BE2vgBuHQamTYr7LROlZLagfQgspL5CqDlXVTFXNrFHjqOnYE1q7du0YM2YM+/btY9euXYwfPx6AKlWqUKFCBWbNcrl01KhRedt07tyZIUOGcOjQIQCWL1/Onj174h+8MSZ+vnrRJYpmV8GNb8YlUUD8WxabRKS216qoDWz2yrOB+mH16gEbvPJ6BZSXOC1btuSmm26iRYsWNGzYkIsvvjjvvVdffZW+fftSrlw52rdvT6VKlQDo06cPq1evpmXLlqgqNWrU4MMPP/TpCIwxMff5MzD9b3DmtXDdfyEtfnO+SWGXOYr94SKNgAmqepa3/jSwLayDu6qq/kFEzgRGAK1xHdzTgCaqmisic4HfALNxHdwvqurHRe07MzNTj3z40ffff8/pp58evQOMk927d1O+vLseOWjQIDZu3Mjzzz8f8fbJetzGGI8qzBgEnw2Cs2+Ea4ZAWmy+64vIPFXNPLI8Zi0LERkJtAeqi0g28CgwCBgtIr2BtcANAKq6RERGA98BOUB/704ogHtwd1aVwd0FVew7oZLNRx99xJNPPklOTg4NGzbkjTfe8DskY0y8qMLUR+HL56HFrdDtRQikxT2MmLYs/FSSWhbFlarHbUzSU4VJA2D2fyCzN3R9BmI8iDfuLQtjjDHFEMyFCb+D+cOgTX/o/AT4OCbZkoUxxiSanAPwQV/4bixc/CBc9mdfEwVYsjDGmMRyYBeMuhVWfQadnoC29/kdEWDJwhhjEseuTTDyJtj4LVzzH2hxs98R5bFkEWfly5dn9+7dfodhjEk0m76DETfC3m3QYwSc1sXviA5jycIYY/z2w3QY3RPSy0Cvj6HOuX5HdJREme4j5ezevZsOHTrQsmVLzj77bMaOHQvAnj17uPLKKznnnHM466yzeOeddwAYMGAAZ5xxBs2bN+fBBx8EYM2aNXTo0IHmzZvToUMH1q5d69vxGGNO0Lw34O3roVJ96DMtIRMFpHLLYuIA+GlRdD/zV2fDFYMiqpqRkcGYMWOoWLEiW7dupU2bNnTr1o1JkyZRp04dPvroIwB27tzJ9u3bGTNmDEuXLkVE8p6sd99993HHHXfQs2dPXnvtNe6//36b7sOYZBEMwtRH3FxPp3aE61+HjIp+R3VM1rLwiaryxz/+kebNm9OxY0fWr1/Ppk2bOPvss5k6dSoPP/wwM2fOpFKlSlSsWJGMjAz69OnDBx98QNmyZQH4+uuv86Ysv/322/niiy/8PCRjTKQO7oHRt7tE0aov3PxOQicKSOWWRYQtgFgZPnw4W7ZsYd68eaSnp9OoUSP2799P06ZNmTdvHh9//DEDBw6kU6dOPPLII8yZM4dp06YxatQoXnrpJaZPn37UZyb5QwSNSQ2/bIARN8GmxdDlKWhzt98RRSR1k4XPdu7cSc2aNUlPT+fTTz9lzZo1AGzYsIGqVaty2223Ub58ed544w12797N3r176dq1K23atOHUU08FoG3btowaNYrbb7+d4cOHc9FFF/l5SMaYoqyfDyNvdi2Lm9+Bpp38jihilix8cuutt3L11VeTmZlJixYtaNasGQCLFi3ioYceIhAIkJ6ezpAhQ9i1axfdu3dn//79qCrPPfccAC+88AJ33XUXTz/9NDVq1OD111/385CMMYVZMgbG3APlakDvMVDrDL8jOi42kWAKSNXjNiYhqMJnT8GMJ6F+G/es7PKJ+3A2m0jQGGPi7eBeGHuva1Wccwtc/e+4Pdku2ixZGGNMLOxcD6NudlN3XP5XaHu/75MBFkfKJQtVTam7hkrqZUZjEtq6ufDOra5lcfOohJu640Sk1DiLjIwMtm3bljInUFVl27ZtZGRk+B2KMalj4Uh440o3dUefKSUiUUCKtSzq1atHdnY2W7Zs8TuUuMnIyKBevXp+h2FMyZeb4x5/+vVL0OhiuGEYlKvmd1RRk1LJIj09ncaNG/sdhjGmpNn3M7zXG36YBq37Qed/QFq631FFlS+XoUTkARFZLCJLROS3XllVEZkiIiu81yph9QeKyEoRWSYinf2I2RhjCrT5e/jvZbDqc7j6eej6dIlLFOBDshCRs4C+QGvgHOAqEWkCDACmqWoTYJq3joicAfQAzgS6AINFJC3ecRtjzFG+nwCvdIQDu+HOCXDenX5HFDN+tCxOB2ap6l5VzQE+A64FugPDvDrDgGu85e7AKFU9oKqrgJW4RGOMMf4IBuHTf7g7nqo3hV9/Bg3a+B1VTPmRLBYD7USkmoiUBboC9YFaqroRwHut6dWvC6wL2z7bKzuKiPQTkSwRyUqlTmxjTBzt2+HGT3z2FLS4FXpNhIp1/I4q5uLewa2q34vIU8AUYDfwDZBTyCYFDYoo8N5XVR0KDAU33UcxQzXGmMNtWQ4je8CONdD1GWjVJ6kH2h0PXzq4VfVVVW2pqu2A7cAKYJOI1AbwXjd71bNxLY+QesCGeMZrjDGs+hxe7QgHfoGeE6B135RJFODf3VA1vdcGwHXASGAc0NOr0hMY6y2PA3qIyEki0hhoAsyJb8TGmJS2cAS8dR1UqO0efdrwAr8jiju/xlm8LyLVgENAf1X9WUQGAaNFpDewFrgBQFWXiMho4Dvc5ar+qprrU9zGmFSiCp8/A5/+HRpfAje+CWUq+x2VL3xJFqp6cQFl24AOx6j/BPBErOMyxpg8wSBM/hPMGgzNe0C3F6FUab+j8k1KjeA2xpiI5ObA+Pth4XBo/WvoMggCKTWV3lEsWRhjTLhD++H93rB0ArQfCJc8nFId2cdiycIYY0L2/wKjboHVM6HLU9Dmbr8jShiWLIwxBmD3Fhj+P7BpCVz3X2h+o98RJRRLFsYYs2EhvHM77NniHlbU5HK/I0o4qd1jY4wxC0fCa51Bg9DrI0sUx2AtC2NMajq4Fz75I8x7HRq3g+tfh3LV/Y4qYVmyMMaknvXz4IN+sG0ltL0fOjwKaXY6LIz97RhjUkduDnzxLMwYBBV+BXeMg5Mv8TuqpGDJwhiTGrYshzG/hg3z4azr4cpnoEyVorczgCULY0xJFwzCnJdh6mOQXsb1TZx1nd9RJR1LFsaYkuvnNTC2vxtk16QzdHvBXX4yx82ShTGm5FGF+cPgkz8BAt1egnNvs2k7isGShTGmZNmZDePuhx+muVtiu/8fVG7gd1RJz5KFMaZkUIUFb7nWRDDHPfY0s3fKzxYbLZYsjDHJb8c6N6X4D9Oh4UXQ/SWo2tjvqEoUSxbGmOQVDLoR2FMecS0La03EjCULY0xy2rrStSbWfOkeedrtBajSyO+oSixf0q+I/E5ElojIYhEZKSIZIlJVRKaIyArvtUpY/YEislJElolIZz9iNsYkiNxDMPNfMKQt/LTYPe70jrGWKGIs7slCROoC9wOZqnoWkAb0AAYA01S1CTDNW0dEzvDePxPoAgwWkbR4x22MSQDr58HQS2HaX6FpZ7hvDrS8w26JjQO/LuyVAsqISCmgLLAB6A4M894fBlzjLXcHRqnqAVVdBawEWsc3XGOMrw7ucXc5vdLRPXPiprfhprdsgF0cxb3PQlXXi8gzwFpgHzBZVSeLSC1V3ejV2SgiNb1N6gKzwj4i2ys7ioj0A/oBNGhg91UbUyIs/wQ++j3sXAeZd0HHxyCjkt9RpRw/LkNVwbUWGgN1gHIiclthmxRQpgVVVNWhqpqpqpk1atQofrDGGP/s+glG94QRN0LpctBrElz1nCUKn/hxN1RHYJWqbgEQkQ+AtsAmEanttSpqA5u9+tlA/bDt6+EuWxljSqJgLsx9Fab/DXIOwKV/hgsfgFKl/Y4spfmRLNYCbUSkLO4yVAcgC9gD9AQGea9jvfrjgBEi8iyuJdIEmBPvoI0xcbBhIUz4LWxYACe3hyufhWqn+ByUAX/6LGaLyHvAfCAHWAAMBcoDo0WkNy6h3ODVXyIio4HvvPr9VTU33nEbY2Jo3w749AmY+wqUrQbXvQJnX293OSUQUS3w8n/Sy8zM1KysLL/DMMYURhW+HQ2T/wx7t0KrPnDpn6BMZb8jS1kiMk9VM48stxHcxhh//LQYPn4I1n4Fdc+DW9+FOi38jsocgyULY0x87dsBM56EOf91dzZ1exFa3GbzOSU4SxbGmPgIBt0U4tMeh30/w3m94LI/Q9mqfkdmImDJwhgTe2tnw8Q/wMaF0OACuOKfULu531GZ42DJwhgTOzvXw9RHYdG7UKGO3eWUxCxZGGOi7+Ae+PIF+PJ50CC0ewgu+p0biW2SkiULY0z0BIOwaDRMfRx2bYAzr4WOj0OVhn5HZorJkoUxJjpWfe7GS2z8Bmq3gOtfg4YX+B2ViRJLFsaY4tm81PVLLJ8ElerDdf+Fs663W2FLGEsWxpgTs+sn+PQf7nbY0uWhw6PQ5h5IL+N3ZCYGLFkYY47P/l/gqxfg6/+D3IPQuh+0+wOUq+Z3ZCaGCk0WIhIArlfV0XGKxxiTqHIOQtar8PnTsHcbnHmdG1Rns8KmhEKThaoGReQ+wJKFMalKFb4f7/oltv8Ijdu5O5zqtvQ7MhNHkVyGmiIiDwLv4J45AYCqbo9ZVMaYxHBgF4y8GVbPhBrN4Nb34NSONqguBUWSLO7yXvuHlSlwcvTDMcYkjGAQxtwNa76Crs+4uZzSrJszVRX5m1fVxvEIxBiTYD5/GpZOgC6DoHVfv6MxPisyWYhIOnAP0M4rmgG8rKqHYhiXMcZPSz+GGf+Ac26G8+/2OxqTACJpUw4B0oHB3vrtXlmfWAVljPHRlmXwQT83Cvuq56x/wgCRJYtWqnpO2Pp0EfkmVgEZY3y072cY2QPSM6DHcBtgZ/JEMh4/V0TybqQWkZOB3BPdoYicJiILw35+EZHfikhVEZkiIiu81yph2wwUkZUiskxEOp/ovo0xhcjNgffugh3r4Ka3oVI9vyMyCSSSlsWDwKci8iMgQEOg14nuUFWXAS0ARCQNWA+MAQYA01R1kIgM8NYfFpEzgB7AmUAdYKqINFXVE05YxpgCTHsMfpgOVz8PDdr4HY1JMEWN4E4DzgGaAKfhksVSVT0Qpf13AH5Q1TUi0h1o75UPw3WkPwx0B0Z5+1wlIiuB1sDXUYrBGLNwJHz1IrTqC+fd6Xc0JgEVehnK+/beTVUPqOq3qvpNFBMFuBbDSG+5lqpu9Pa7EajpldcF1oVtk+2VHUVE+olIlohkbdmyJYphGlOCrZ0F4++HxpdAlyf9jsYkqEj6LL4SkZdE5GIRaRn6Ke6ORaQ00A14t6iqBZRpQRVVdaiqZqpqZo0aNYobojEl389rYNStULkB3DgM0tL9jsgkqEj6LNp6r38NK1PgsmLu+wpgvqpu8tY3iUhtVd0oIrWBzV55NlA/bLt6wIZi7tsYc2CXu/MpeAhufgfKVCl6G5OyIumzGKeqz8Vg3zeTfwkKYBzQExjkvY4NKx8hIs/iOribAHNiEI8xqSM3B97t5cZU3PY+VD/V74hMgouozyLaOxWRssDlwAdhxYOAy0VkhffeIC+GJbhZb78DJgH97U4oY4rpkz/Cyilw5b/glEv9jsYkgUguQ30lIi9x9Kyz8090p6q6F6h2RNk23N1RBdV/AnjiRPdnjAkz+2WY8zK0/Q1knvBd8CbF+NlnYYyJt+WfwKQB0Owq90wKYyIUyayz1kY1piTYsMD1U/yqOVw3FAJpfkdkkkiRt86KSC0ReVVEJnrrZ4hI79iHZoyJmh1rYcRNULYa3DIaSpfzOyKTZCIZZ/EG8AnuTiSA5cBvYxSPMSba9u2A4TfCof1w62ioUMvviEwSiiRZVFfV0UAQQFVzKMZEgsaYOMo5AO/cBttWwk1vQc3T/Y7IJKlIOrj3iEg1vFHTItIG2BnTqIwxxRcMwof3uOdnXzsUTr7E74hMEoskWfwvbmDcKSLyJVADuD6mURljim/qI7D4fej4GJxzk9/RmCQXyd1Q80XkEvJnnV1mj1Q1JsHNGuJmkW3dDy78rd/RmBIgkpZFqJ9iSYxjMcZEw6L33FiK06+GLoPssagmKiLp4DbGJIsfpsOYu6HhRXDdKzaWwkSNJQtjSor182HUbVCjGdw8wj1H25goiWRQ3rRIyowxPtqyHIZfD+WqwW3vQUYlvyMyJcwx+yxEJAMoC1QXkSrkP4SoIvkD9IwxftuxDt66BiQNbv8QKvzK74hMCVRYB/evcSO16wDzyE8WvwD/F9uwjDER2b3FJYoDu6HXR1DtFL8jMiXUMZOFqj4PPC8iv1HVF+MYkzEmEvt3wtvXwc71cPsY+NXZfkdkSrBIOrh/EpEKACLyZxH5IBrP4DbGFMPBPW6+p83fu2k8Gl7gd0SmhIskWfxFVXeJyEVAZ2AYMCS2YRljjik031P2HPif/0KTy/2OyKSASJJFaNLAK4EhqjoWKB27kIwxx5SbA+/d5cZTdHsRzrzW74hMiogkWawXkZeBG4GPReSkCLc7JhGpLCLvichSEfleRC4QkaoiMkVEVnivVcLqDxSRlSKyTEQ6F2ffxiStYC58eDcsnQBdnoJzb/M7IpNCIjnp34h7nkUXVd0BVAUeKuZ+nwcmqWoz4Bzge2AAME1VmwDTvHVE5AygB3Am0AUYLCI2LNWklmAQxj8Ai96FDo9Cm7v9jsikmCKTharuBTYDF3lFOcCKE92hiFQE2gGvep9/0EtC3XH9IXiv13jL3YFRqnpAVVcBK4HWJ7p/Y5KOKkx6GBa8Be3+ABf/r98RmRQUyQjuR4GHgYFeUTrwdjH2eTKwBXhdRBaIyCsiUg6opaobAbzXml79usC6sO2zvbKCYu0nIlkikrVly5ZihGhMglCFT/4Ic4bCBffBpX/0OyKToiK5DHUt0A3YA6CqG4AKxdhnKaAlrrP8XO9zBxRSv6ApM7Wgiqo6VFUzVTWzRo0axQjRmASgClP+ArMGw/n3QKe/2wyyxjeRJIuDqqrkPymvuE96zwayVXW2t/4eLnlsEpHa3j5q4y59herXD9u+HrChmDEYk9hUYdrj7pkUrfpClyctURhfRZIsRnt3Q1UWkb7AVOCVE92hqv4ErBOR07yiDsB3uKfx9fTKegJjveVxQA8ROUlEGgNNgDknun9jEl4oUXzxHJzXC674pyUK47tInpT3jIhcjpsT6jTgEVWdUsz9/gYYLiKlgR+BXrjENVpEegNrgRu8/S8RkdG4hJID9FfV3II/1pgkpwpTH4Mv/w2Zd0HXf0HAniRg/CfuClMhFUSeUtWHiypLNJmZmZqVleV3GMZEThWmPAJfvQCZvaHrM5YoTNyJyDxVzTyyPJJ/iQXNJXBF8UMyxuRRhUkDXaJo1QeutBaFSSyFPc/iHuBe4GQR+TbsrQrAl7EOzJiUEQzCxw9C1qvurifrzDYJqLA+ixHAROBJDr+1dZeqbo9pVMakimCuG5m94C248LfQ8TFLFCYhFfY8i53ATuDm+IVjCnRwD+zeBPt+dj+1W0C56n5HZYor9xCM+TUsft+NzL70j5YoTMIq8m4oE0eqsHcbbF8F21bC+nmwbhZsWgIazK9XNxN6T7Fr2sns0H54905YPhE6Pg4X/dbviIwplCULPx3cCxN+Bz+vgt2bYc8WOLg7//3S5aHueXDxg1D1ZChTBTZ/5+7B/2aEzTqarA7shlG3wKrPXUd2qz5+R2RMkSxZ+Gn5JPh2FNRv45JCuRpQuYFLDKGftCN+RU06ue2mPgbNroIylf2I3Jyovdth+PWwYSFc+x84p4ffERkTEUsWflr6EZStBr0+hkCEs64HAm5E79D2MGMQXDEopiGaKNq5Ht66Fn5eDTe9Dc26+h2RMRGzi95+yT0EK6ZA0ysiTxQhdVpAZi83E+mm72ISnomyrSvgtS6wayPc/oElCpN0LFn4ZfUXcGDniZ80LvsLZFR09+cXMQrf+Cw7C17tBDn7oOd4aHRR0dsYk2AsWfhl6UdQqgycfOmJbV+2Klz+V1jzJSwcHt3YTPQs/wSGXQ0ZlaD3ZNcqNCYJWbLwgyosmwinXAaly57457S4DRq0hcl/hj1boxefiY75b8HIm6F6E5coqp7sd0TGnDBLFn7Y+A38kl3869aBAFz1nLsV85M/RSc2U3yqMP0JGHcfnHwJ3PkRlK9Z9HbGJDBLFn5Y+hFIAJp2Kf5n1WzmBnR9Owp+nFH8zzPFk3MQPrwXPv+na/ndMhpOKs6DJY1JDJYs/LDsYze2IlpTdlz8e3eJY/wDbmoQ4499P8Pb17kBk+0HQveXIC3d76iMiQpLFvH28xrYtBiaXRm9z0wvA91edPfvT/tb9D7XRG77j/DK5bB2FlwzBNoPsHmeTIliySLeVs90r6d2jO7nNroIWveD2f+BNV9F97NN4dbOglc6uula7vgQWtzid0TGRJ0li3hb+7Wb46l60+h/dodH3XQhY/u7eadM7C0YDm9c5W6N7TPNxlCYEsuSRbytnQ31z4/NjLEnlXfXybf/CNPtclRMBXPdLctj74WGbV2iqH6q31EZEzO+JAsRWS0ii0RkoYhkeWVVRWSKiKzwXquE1R8oIitFZJmIdPYj5qjYsxW2rYAGbWK3j8btoFVfmDUYfvg0dvtJZft2wMge8NWL7lnZt73vBkkaU4L52bK4VFVbhD0YfAAwTVWbANO8dUTkDKAHcCbQBRgsIsc5mVKCWDfbvdaPYbIAN7K7elP48B43y6mJni3L4ZUO8MN0N734Vc/aHU8mJSTSZajuwDBveRhwTVj5KFU9oKqrgJVA6/iHFwVrv4a00lDn3Njup3RZ+J9XXEtm3G9s7qhoWTbRJYp9O+COcfYcCpNS/EoWCkwWkXki0s8rq6WqGwG819CQ17rAurBts72yo4hIPxHJEpGsLVu2xCj0Ylg72yWK9IzY76v2OdDxUVg6Aea/Gfv9lWTBXJj+d3fpqWpj6DcDGl3od1TGxJVfz7O4UFU3iEhNYIqILC2kbkE3qxf4VVlVhwJDATIzMxPr6/ShfbBhAbS5J377bNMfVk6FSQOgfmuoeXr89l1S7N0OH/R1f48tboMrn3HjWkyJFAwqOUElN6jkBIPkBpVDuYev5wSVnNzC13O99bzPClt3y0pObv56bui9o+qG11EO5XoxHfZefpyh7SY90I7SpaLbFvAlWajqBu91s4iMwV1W2iQitVV1o4jUBjZ71bOB+mGb1wM2xCq2xet3cjA3SOm0AKXShPS0AOkBt1wqTfKW09MCpAWEUgFBIhl8tWEBBA9BgwtiFfrRAgG49mX4z8Uw+g7o+6m7Y8pEZv08GH2newbFVc/Beb2SfqBdMKjkqjvxBEOvQdyJTt2ye3XvhZZDJ73QNrmHve+2d+9x2Hu5waBX5r2qkpsbJFc57L2cYNh+NHSy1byTb34M+duGTqCh48k54qSeq/ll+Sf14GGfG34s7hj8+92kp4l3TvHON0cupwW8VyEt4C0HhLKlS5EWENLT8utrwd+niyXuyUJEygEBVd3lLXcC/gqMA3oCg7zXsd4m44ARIvIsUAdoAsyJVXy/fWchKzfvLrpimFDSKBXwftlp+b/INO8XeMvB9+gL9JgYZM+UL0jz6qYFhDRxr4GAkCbu8wJhZQFx5fnLQiBAfh0RRNx6wHsVb1kEGjZ6jOuX3MfyV+7ik6Z/JxBw9cXbTghtn78M3vvknx+F8G047MQZWgp9XkHCz7PH6kZRNO+9w6po/j9/VVCvkobWCSvT/M/RI9bBnTBD5cHQ56oS9OoFg8q5P73HZWufZ096NT44bQjZP51JcMJ3efsOetsG89ZDy4ev56qi6vbtTrb57+eG1w3m1w16dcOXg5p/olbvc/OWvROjhp3IC9rOzxNhpEp5/+bTw/9/eCfG0HqozpFlpQIBAgEoW6rU4eVp7v9I+Be8Ummh/3cB0gLk/Z919dz/4TSRvBN1WthJO/yEnl8nkHfCTgurk54WOp7DT/rhcYS2TXR+tCxqAWO8b+OlgBGqOklE5gKjRaQ3sBa4AUBVl4jIaOA7IAfor6q5sQruqf85m1/253AoJ5jXrDuU676R5DX9cpVDwSC5ufnNwfBvQ+HN19A3ojbrVvCT1qd81V9R5ohvRznBIAdyjv5Pnv+fPew/fjD/JBQMez/0nkLeyUcJvVeZNYHreWjLaN7aUJfhuVEePV5CBARqyE4eL/U6nQJz+DR4LgNy+vPL4gqkBdZ5CTI/mQuhBO4SZCAsAReUxAXvC0EgbNnbJi1vO3fiCG2XFsj/AhD6zNA+08KXvRZuWni9UBwBObyuhOLILy8VVi8t7ItI6LMOez/t8P2ETsb5n0feCTDtiG2OfC90gg6tm8QloW9jJU1mZqZmZWX5HYYTDMI/G8PpV7tBcz7FoCNugFWfk3v7WIL12uR/+w775h36pq0KeOVw+Lf3vPc9ed/59RidSUfUDznWFR3J++PwVkpei4b8Vk9+uXszVCe8VXR4a4m8E31+60tcgIvehYl/cJMxXvZnuOA3sRk8aUwCE5F5YUMa8vjVwZ1ati6H/TtiOxivKIEAct1/4ZUOlBp9O/SdDlUa+hdPItm7HcbeB8s+gnqtXUKvcZrfURmTUOxrUzysm+VeYz0Yryhlq7rnK+QecreBHtjlbzyJYPWXMORCWDkFOj0Bd02yRGFMASxZxMO6uVCmKlQ7xe9I3CM+b3wDtiyD93pDbo7fEfkjmAsznoJhV7lbYXtPgbb3QSA5JwcwJtYsWcRD9hw3ziFRbrs85TLo+k9Y8QmMv9/1qaSSXzbCm91hxj/g7Bvg159BnRZ+R2VMQrM+i1jbu931WTS/ye9IDteqj5sOZMaTcFJF6PJk4iSzWFo+GT682w2SvGaIPXvCmAhZsoi1bO+OrPoJOJ3VJQ+7eY5mD3HPY7h0oN8Rxc6h/TD1UfdwqFpnwQ1vuEtyxpiIWLKItew5IAGo09LvSI4mAp3/AQd+gc8GudlT2z3od1TRt+k7eL8PbF4C598NHR+Pz/xcxpQglixibd0c9002UafZCATg6hcg96D3wCSFdg/5HVV0BIOu1TT1ccioCLe+B00u9zsqY5KSJYtYCua6+YXO6eF3JIVLK+XmkJKAm11VFS75g99RFc+Ode55HqtnQtMroNsLUL5m0dsZYwpkySKWNn8HB3e7gV6JLpDmOnwR+PQJF3fHx5Ov01sVFg6HSQNBg9DtRTj39uQ7DmMSjCWLWFrnzXeYiJ3bBQmkwTWD3cOTvnwe9m6Dq553LY9k8MtGdyvwisnQoK07lqqN/Y7KmBIhSc4CSSp7LpSrAVUa+R1J5AJpcOWzUK6m6/Teu909da90Ob8jOzZVWDgCPhkIOQehy1PQup/N62RMFFmyiKV1s90lqGS7BCLibqMtW81NrPfyJS5hJOLAtZ/XwPgH4MdPXWui+0uJMVLemBLGvnrFyp6tsP1HqN/K70hO3Pn94I6xbhbWVzq6S1OJMj2IKswbBoMvcC24rs/AnR9ZojAmRixZxEr2XPda/3x/4yiuky+Be76E07rAlEfg32fB9Cdgx1r/Ytq7Hd65zfVP1DsP7p0FrfvaZSdjYsguQ8XKmi8hrTTUOdfvSIqvbFW48S1YPgmyXoPPn3Y/jS92cyud3g3KVI5PLCunwof9Yd926PR395xxSxLGxJwli1hZ/QXUa+VmNC0JROC0K9zPjrWwYDgsGg3jfgMf/R5O6eAe7nTaFS65RNvBPTD5L5D1KtRoBre+C7WbR38/xpgCWbKIhf07YeM3JWck9JEqN3Ad4O0HwIb5sOg9+H48LJ8IkgYNLoAmHeHUy6HWmcXv4F/9hUtK21fBBffBZX+x6TqMiTPfkoWIpAFZwHpVvUpEqgLvAI2A1cCNqvqzV3cg0BvIBe5X1U98CTpSa2e5AWGNLvI7ktgSgbrnuZ/O/4CNC13SWDEZpj7mfsrVcMmj4YXuSYG1zop83Maera418c0Il6DunFDy/06NSVB+tiweAL4HKnrrA4BpqjpIRAZ46w+LyBlAD+BMoA4wVUSaqmquH0FHZPVM119RL4nvhDpeIq5/ps650OER+GWD619Y/QWs+Rq+H+fqlS7vkku9TJc4ftXcDZwLPXQoGHRTpCwdD/PfhAO74eLfw8UPusGCxhhf+JIsRKQecCXwBPC/XnF3oL23PAyYATzslY9S1QPAKhFZCbQGvo5jyMdn1cyS1V9xIirWgZZ3uB9wczWtm+1aXetmwRf/hlC+lzTXz1G2mpsyffdPECjl+kEu/yvUbObXURhjPH61LP4N/AGoEFZWS1U3AqjqRhEJzfpWF5gVVi/bKzuKiPQD+gE0aNAgyiFHaN8O+OlbaJfkE/FFW+X67ufs6936of2wZSn8tAh+XuVuh9233SWJJp2haScoU8XfmI0xeeKeLETkKmCzqs4TkfaRbFJAmRZUUVWHAkMBMjMzC6wTc6nSX1Fc6RluRHgijgo3xhzFj5bFhUA3EekKZAAVReRtYJOI1PZaFbWBzV79bKB+2Pb1gA1xjfh4rJ4JaSelVn+FMabEi/toJlUdqKr1VLURruN6uqreBowDenrVegJjveVxQA8ROUlEGgNNgDlxDjtyeeMr7NZOY0zJkUhDXwcBl4vICuBybx1VXQKMBr4DJgH9E/ZOqFB/hV2CMsaUML4OylPVGbi7nlDVbUCHY9R7AnfnVGKz/gpjTAmVSC2L5PfDNChVxvorjDEljiWLaFGFZRPhlMusv8IYU+JYsoiWnxbBznVuIj1jjClhLFlEy7KJgEDTzn5HYowxUWfJIlqWfez6KsrXLLquMcYkGUsW0bBzvZtxtVlXvyMxxpiYsGQRDcsnutfTLFkYY0omSxbRsGwiVD0Zqjf1OxJjjIkJSxbFdWAXrPrctSqK+0Q4Y4xJUJYsimvlNMg9aJegjDElmiWL4pr/pnt0aP3z/Y7EGGNixpJFcayf76b4uKB/5M+VNsaYJGTJojhm/gsyKkNmb78jMcaYmLJkcaI2LYGlE+D8uyGjot/RGGNMTFmyOFEz/wWly8P5v/Y7EmOMiTlLFidi60pYMgZa9YGyVf2OxhhjYs6SRbjcHPjkTzD75WPXCebCJwMhrbTr2DbGmBRgySJcIA22Loepj8OOtQXXmfY4rJgMnZ+wSQONMSnDkkU4EbjyWbc84X/dA43CLRwBXz7vLj+16hP/+IwxxidxTxYikiEic0TkGxFZIiKPe+VVRWSKiKzwXquEbTNQRFaKyDIRie0DIyrXhw5/gZVTYPH7+eU/fgbjH4DG7aDLoJiGYIwxicaPlsUB4DJVPQdoAXQRkTbAAGCaqjYBpnnriMgZQA/gTKALMFhE0mIaYet+UKclTBrg5n0aeQu82Q0q1YcbhkFaekx3b4wxiSbuyUKd3d5quvejQHdgmFc+DLjGW+4OjFLVA6q6ClgJtI5pkIE06PYC7N0Ow66GNV/ApX+GfjPs7idjTEryZY4Kr2UwDzgV+D9VnS0itVR1I4CqbhSRUO9xXWBW2ObZXllBn9sP6AfQoEGD4gX5q7Oh24uwawO06gtlKhfv84wxJon5kixUNRdoISKVgTEiclYh1Qua91sLKENVhwJDATIzMwusc1zOvbXYH2GMMSWBr3dDqeoOYAauL2KTiNQG8F43e9Wygfphm9UDNsQvSmOMMX7cDVXDa1EgImWAjsBSYBzQ06vWExjrLY8DeojISSLSGGgCzIlr0MYYk+L8uAxVGxjm9VsEgNGqOkFEvgZGi0hvYC1wA4CqLhGR0cB3QA7Q37uMZYwxJk5Ejxx4VkJkZmZqVlaW32EYY0xSEZF5qpp5ZLmN4DbGGFMkSxbGGGOKZMnCGGNMkSxZGGOMKVKJ7eAWkS3AmuPYpDqwNUbhJKpUPGZIzeNOxWOG1Dzu4h5zQ1WtcWRhiU0Wx0tEsgq6A6AkS8VjhtQ87lQ8ZkjN447VMdtlKGOMMUWyZGGMMaZIlizyDfU7AB+k4jFDah53Kh4zpOZxx+SYrc/CGGNMkaxlYYwxpkiWLIwxxhQppZKFiHQRkWUislJEBhTwvojIC97734pISz/ijLYIjvtW73i/FZGvROQcP+KMpqKOOaxeKxHJFZHr4xlfrERy3CLSXkQWisgSEfks3jFGWwT/viuJyHgR+cY75l5+xBlNIvKaiGwWkcXHeD/65zJVTYkfIA34ATgZKA18A5xxRJ2uwETc0/naALP9jjtOx90WqOItX5Hsxx3JMYfVmw58DFzvd9xx+l1Xxk3338Bbr+l33HE45j8CT3nLNYDtQGm/Yy/mcbcDWgKLj/F+1M9lqdSyaA2sVNUfVfUgMArofkSd7sCb6swCKoee3pfEijxuVf1KVX/2VmfhnkaYzCL5XQP8Bnif/KcyJrtIjvsW4ANVXQugqsl+7JEcswIVRESA8rhkkRPfMKNLVT/HHcexRP1clkrJoi6wLmw92ys73jrJ5niPqTfuG0kyK/KYRaQucC3wnzjGFWuR/K6bAlVEZIaIzBORO+IWXWxEcswvAafjHse8CHhAVYPxCc83UT+X+fGkPL9IAWVH3jccSZ1kE/ExiciluGRxUUwjir1IjvnfwMOqmuu+cJYIkRx3KeA8oANQBvhaRGap6vJYBxcjkRxzZ2AhcBlwCjBFRGaq6i8xjs1PUT+XpVKyyAbqh63Xw33TON46ySaiYxKR5sArwBWqui1OscVKJMecCYzyEkV1oKuI5Kjqh3GJMDYi/Te+VVX3AHtE5HPgHCBZk0Ukx9wLGKTuYv5KEVkFNAPmxCdEX0T9XJZKl6HmAk1EpLGIlAZ6AOOOqDMOuMO7k6ANsFNVN8Y70Cgr8rhFpAHwAXB7En/DDFfkMatqY1VtpKqNgPeAe5M8UUBk/8bHAheLSCkRKQucD3wf5zijKZJjXotrSSEitYDTgB/jGmX8Rf1cljItC1XNEZH7gE9wd1C8pqpLRORu7/3/4O6K6QqsBPbivpEktQiP+xGgGjDY+6ado0k8U2eEx1ziRHLcqvq9iEwCvgWCwCuqWuDtl8kgwt/134A3RGQR7vLMw6qa1NOWi8hIoD1QXUSygUeBdIjducym+zDGGFOkVLoMZYwx5gRZsjDGGFMkSxbGGGOKZMnCGGNMkSxZGGOMKZIlC1MieTOrTjjObeqIyHsnsK/KInJvcT8nmXh/v239jsPEjyULkxC8wUO+/XsUkVKqukFVT2Sq8spAXrIoxudElYjEchxVe9xsxRGLcTwmxixZGN+ISCMR+V5EBgPzgfoi8pCIzPXm4H88rO5fRGSpiEwRkZEi8qBXPkNEMr3l6iKyuoD9tBb3nI4F3utpXvmdIvKuiIwHJnvxLPbee0XcMx8WisgWEXlURMqLyDQRmS8ii0QkNLvpIOAUr+7TR3xOhoi87tVf4M2/Fdr3ByIySURWiMg/j/F3tFpEnhKROd7PqV751SIy2/vMqd7IZETkMREZKiKTgTe9WGZ6Mc8PtQa8lsFnIjJaRJaLyCBxzzWZ48V6ilevhoi87/1O5orIhSLSCLgb+J13zBcXVK+geE7wn4pJBH7Py24/qfsDNMKNIm7jrXfCPWxecF9kJuDm7c/ETQRXBqgArAAe9LaZAWR6y9WB1d5ye2CCt1wRKOUtdwTe95bvxM2hUzUsnsVHxNgQWOq9lgIqhu1rpRfrYduFrwO/B173lpvhpp7I8Pb9I1DJW18D1C/g72g18Cdv+Y6wY6pC/qDaPsC/vOXHgHlAGW+9LJDhLTcBssL+fnYAtYGTgPXA4957DwD/9pZHABd5yw2A78P282BYnIXVy4vHfpL3x5qFxm9r1M23Dy5ZdAIWeOvlcSe4CsBYVd0H4LUEjkclYJiINMHNvJke9t4UVS3wuQAikgG8C9ynqmtEJB34h4i0wyW5ukCtIvZ9EfAigKouFZE1uGnCAaap6k5vX9/hEtK6Aj5jZNjrc95yPeAdcc8oKA2sCqs/LvR35R3rSyLSAsgN2zfAXPXmCxKRH4DJXvki4FJvuSNwhuTPzFtRRCoUEGNh9cLjMUnKkoXx256wZQGeVNWXwyuIyO8K2T6H/MupGceo8zfgU1W91ruEMuMY+z/Sf3APCprqrd+Ke9Laeap6yLvkdax9hhQ2//mBsOVcjv3/UQtYfhF4VlXHiUh73Df4kPBj+h2wCTezbADYf4z9B8PWg2GxBIALjjzZy9HTuhdWr7C/Y5MkrM/CJJJPgLtEpDy4BxSJSE3gC+Bq7/p/eeDKsG1W457PAHCsTuVKuMss4C7/FElE+gMVVHXQEZ+z2UsUl+JaAgC7cK2fgnyOSzKISFPcJZplkcQQ5qaw16/DYgkdU89Ctq0EbFT3sJ/bcZPtHY/JwH2hFa+FAkcf87HqmRLCkoVJGKo6GXft+2txM4S+hzthz8VNufwNbir1LGCnt9kzwD0i8hWuH6Eg/wSeFJEvifxk+SBwdlgn993AcCBTRLJwCWCpF/c24EsRWSwiTx/xOYOBNO943gHuVNUDHJ+TRGQ2ri8h1Mp6DHhXRGYChc2gOhjoKSKzcJegjvdb/v24Y/7Wu1R2t1c+Hrg21MFdSD1TQtissyYpiEh5Vd0t7hkMnwP9VHW+33HFmnepK1OTfEptk/ysz8Iki6Eicgauj2BYKiQKYxKJtSyMMcYUyfosjDHGFMmShTHGmCJZsjDGGFMkSxbGGGOKZMnCGGNMkf4fx2nOkkCIrT4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cross validate over regularization parameter alpha\n",
    "np.random.seed(6)\n",
    "\n",
    "w = randn(40)\n",
    "X,y = generate_noisy_data(30, w)\n",
    "Xtest,ytest = generate_noisy_data(30, w)\n",
    "\n",
    "ridge_error = []\n",
    "lasso_error = []\n",
    "alphas = np.linspace(.01,1,100)\n",
    "for alpha in alphas:\n",
    "    yhat_ridge = Ridge(alpha = alpha).fit(X,y).predict(Xtest)\n",
    "    ridge_error.append(sum((ytest - yhat_ridge)**2))\n",
    "    yhat_lasso = Lasso(alpha = alpha).fit(X,y).predict(Xtest)\n",
    "    lasso_error.append(sum((ytest - yhat_lasso)**2))\n",
    "plt.plot(alphas, ridge_error, label=\"ridge\")\n",
    "plt.plot(alphas, lasso_error, label=\"lasso\")\n",
    "plt.ylabel(\"test error\")\n",
    "plt.xlabel(\"regularization parameter\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[486.8744087820671,\n",
       " 486.6381095030756,\n",
       " 486.4137391169036,\n",
       " 486.20096011523566,\n",
       " 485.99944655883013,\n",
       " 485.8088836064389,\n",
       " 485.6289670659078,\n",
       " 485.4594029662677,\n",
       " 485.2999071497014,\n",
       " 485.1502048823566,\n",
       " 485.01003048299606,\n",
       " 484.8791269685702,\n",
       " 484.7572457158444,\n",
       " 484.6441461382387,\n",
       " 484.5395953771271,\n",
       " 484.4433680068402,\n",
       " 484.3552457527017,\n",
       " 484.27501722143006,\n",
       " 484.2024776433035,\n",
       " 484.13742862550004,\n",
       " 484.07967791605574,\n",
       " 484.0290391779454,\n",
       " 483.9853317727668,\n",
       " 483.9483805535866,\n",
       " 483.91801566649457,\n",
       " 483.8940723604596,\n",
       " 483.8763908050878,\n",
       " 483.86481591591235,\n",
       " 483.85919718686824,\n",
       " 483.8593885296015,\n",
       " 483.8652481193156,\n",
       " 483.87663824683636,\n",
       " 483.89342517662715,\n",
       " 483.9154790104646,\n",
       " 483.942673556536,\n",
       " 483.974886203705,\n",
       " 484.0119978007142,\n",
       " 484.05389254010595,\n",
       " 484.1004578466554,\n",
       " 484.1515842701106,\n",
       " 484.20716538205664,\n",
       " 484.26709767672173,\n",
       " 484.3312804755535,\n",
       " 484.3996158354025,\n",
       " 484.47200846016295,\n",
       " 484.54836561571625,\n",
       " 484.62859704804316,\n",
       " 484.7126149043644,\n",
       " 484.80033365718964,\n",
       " 484.8916700311481,\n",
       " 484.98654293248774,\n",
       " 485.0848733811303,\n",
       " 485.18658444517786,\n",
       " 485.29160117777525,\n",
       " 485.3998505562207,\n",
       " 485.5112614232481,\n",
       " 485.625764430383,\n",
       " 485.74329198329036,\n",
       " 485.86377818903924,\n",
       " 485.98715880520393,\n",
       " 486.11337119073085,\n",
       " 486.24235425849895,\n",
       " 486.37404842951435,\n",
       " 486.5083955886623,\n",
       " 486.6453390419739,\n",
       " 486.7848234753312,\n",
       " 486.9267949145698,\n",
       " 487.07120068691495,\n",
       " 487.2179893837042,\n",
       " 487.36711082435613,\n",
       " 487.51851602152,\n",
       " 487.67215714737836,\n",
       " 487.82798750105434,\n",
       " 487.98596147707656,\n",
       " 488.1460345348729,\n",
       " 488.3081631692475,\n",
       " 488.47230488180827,\n",
       " 488.63841815330585,\n",
       " 488.80646241686213,\n",
       " 488.9763980320385,\n",
       " 489.14818625973265,\n",
       " 489.3217892378549,\n",
       " 489.49716995777806,\n",
       " 489.67429224150976,\n",
       " 489.85312071958117,\n",
       " 490.0336208096155,\n",
       " 490.2157586955583,\n",
       " 490.399501307541,\n",
       " 490.5848163023606,\n",
       " 490.77167204455117,\n",
       " 490.96003758803164,\n",
       " 491.14988265829885,\n",
       " 491.3411776351606,\n",
       " 491.5338935359841,\n",
       " 491.72800199943777,\n",
       " 491.92347526972,\n",
       " 492.1202861812488,\n",
       " 492.31840814380445,\n",
       " 492.51781512810277,\n",
       " 492.71848165179154]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ridge_error"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}