lasso.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CS 531 Spring 2020\n",
    "# Author: Alina Ene (aene@bu.edu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applications in statistics and machine learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Least squares, least absolute shrinkage and selection operator (LASSO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be using cvxpy to solve the LASSO optimization. The default solver for cvxpy is numerically unstable, so we will use the CVXOPT solver instead. We can specify the solver to cvxpy as follows:\n",
    "\n",
    "$\\texttt{prob.solve(solver=CVXOPT)}$\n",
    "\n",
    "If you are running this notebook on your machine, you will need to install cvxpy and cvxopt. You can install them using conda as follows:\n",
    "\n",
    "$\\texttt{conda install -c cvxgrp cvxpy libgcc}$\n",
    "\n",
    "$\\texttt{conda install -c anaconda cvxopt}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import cvxpy as cvx\n",
    "import numpy as np\n",
    "from numpy import linalg\n",
    "from cvxpy import *\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Least squares\n",
    "Recall that, in the least squares problem, we are given $A \\in \\mathbb{R}^{m \\times n}$, $b \\in \\mathbb{R}^m$, and the goal is to solve the problem\n",
    "$$\n",
    "\\min_{x \\in \\mathbb{R}^n} \\|Ax - b\\|^2 = \\min_{x \\in \\mathbb{R}^n} \\sum_{i = 1}^m (a_i^T x - b_i)^2 \n",
    "$$\n",
    "where $a_i^T$ is the $i$-th row of $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# generate a random instance\n",
    "def random_instance(m, n):\n",
    "    A = np.random.randn(m,n)\n",
    "    b = np.random.randn(m,1)\n",
    "    return (A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# solve a least squares instance using CVX\n",
    "def solve_least_squares(A, b):\n",
    "    x = cvx.Variable(len(A[0]))\n",
    "    obj = cvx.Minimize(cvx.norm(A*x - b))\n",
    "    prob = cvx.Problem(obj)\n",
    "    prob.solve(solver=CVXOPT)\n",
    "    return (prob.value, x, prob.status)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "status: optimal\n",
      "optimal value 3.4295091476370145\n",
      "optimal var [[ 0.12187844]\n",
      " [-0.1258244 ]\n",
      " [ 0.0444059 ]\n",
      " [ 0.20489149]\n",
      " [ 0.30398088]\n",
      " [-0.19723836]\n",
      " [ 0.20083373]\n",
      " [ 0.52730448]]\n"
     ]
    }
   ],
   "source": [
    "# run least squares on a random instance\n",
    "(A, b) = random_instance(16, 8)\n",
    "(value, x, status) = solve_least_squares(A, b)\n",
    "\n",
    "print(\"status:\", status)\n",
    "print(\"optimal value\", value)\n",
    "print(\"optimal var\", x.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LASSO\n",
    "In the LASSO problem, we are solving a variation of the least squares problem where the objective has an additional term that encourages the solution to be sparse (have a small number of non-zeros).\n",
    "$$\n",
    "    \\min_{x \\in \\mathbb{R}^n} \\|Ax - b\\|^2 + \\gamma \\cdot \\|x\\|^2_1\n",
    "$$\n",
    "$\\gamma > 0$ is a positive value that we choose (more on how to set $\\gamma$ later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# solve a LASSO instance using CVX\n",
    "def solve_lasso(A, b, gamma_value):\n",
    "    gamma = cvx.Parameter(sign=\"positive\", value=gamma_value)\n",
    "    x = cvx.Variable(len(A[0]))\n",
    "    obj = cvx.Minimize((cvx.norm(A*x - b))**2+gamma*(cvx.norm(x,1))**2)\n",
    "    prob = cvx.Problem(obj)\n",
    "    prob.solve(solver=CVXOPT)  \n",
    "    return (prob.value, x, prob.status)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "status: optimal\n",
      "optimal value 4.854015704236237\n",
      "optimal var [[-0.47171799]\n",
      " [ 0.00927019]\n",
      " [-0.02328763]\n",
      " [ 0.63281547]\n",
      " [ 0.01590656]\n",
      " [-0.6108903 ]\n",
      " [-0.39388531]\n",
      " [-0.05401183]]\n"
     ]
    }
   ],
   "source": [
    "# run lasso on a random instance\n",
    "(A, b) = random_instance(16, 8)\n",
    "(value, x, status) = solve_lasso(A, b, 0.1)\n",
    "\n",
    "print(\"status:\", status)\n",
    "print(\"optimal value\", value)\n",
    "print(\"optimal var\", x.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now consider varying $\\gamma$ and see what happens to the least squares error term $\\|Ax - b\\|^2$ and the number of non-zeros. Given the LASSO solution $x$, we will threshold down to $0$ all of the entries of $x$ that are at most $0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gamma_range(start, end, step):\n",
    "    gamma = start\n",
    "    while gamma < end:\n",
    "        yield gamma\n",
    "        gamma += step\n",
    "\n",
    "# set each entry of x that is smaller than 0.1 to 0\n",
    "def round_down(x):\n",
    "    x_rounded = [0 for i in range(len(x))]\n",
    "    x_nonzeros = 0\n",
    "    for i in range(len(x)):\n",
    "        if (float(x[i]) < 0.1):\n",
    "            x_rounded[i] = 0.0\n",
    "        else:\n",
    "            x_rounded[i] = float(x[i])\n",
    "            x_nonzeros += 1\n",
    "\n",
    "    return (x_rounded, x_nonzeros)\n",
    "\n",
    "# returns the ell_2 norm squared of Ax - b\n",
    "def least_squares_evaluate(A, x, b):\n",
    "    m = len(A)\n",
    "    n = len(A[0])\n",
    "    Ax = [0.0 for i in range(m)]\n",
    "    for i in range(m):\n",
    "        for j in range(n):\n",
    "            Ax[i] += (float(A[i][j]) * float(x[j]))\n",
    "    z = [0.0 for i in range(m)]\n",
    "    for i in range(m):\n",
    "        z[i] = Ax[i] - b[i]\n",
    "    value = 0.0\n",
    "    for i in range(m):\n",
    "        value += (float(z[i]) * float(z[i]))\n",
    "    return value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# solve the LASSO instances for varying gamma\n",
    "\n",
    "(A, b) = random_instance(16, 8) \n",
    "\n",
    "gamma_values = [x for x in gamma_range(0.1, 2.0, 0.1)]\n",
    "error_values = [0.0 for x in range(len(gamma_values))]\n",
    "num_nonzeros = [0 for x in range(len(gamma_values))]\n",
    "\n",
    "for i in range(len(gamma_values)):\n",
    "    (value, x, status) = solve_lasso(A, b, gamma_values[i])\n",
    "    (xnew, nnz) = round_down(x.value)\n",
    "    error_values[i] = least_squares_evaluate(A, xnew, b)\n",
    "    num_nonzeros[i] = nnz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xt4VNd57/HvKyEkga6AAAlxMTfb\ngG2wFZv4bsdJ7PiW1ElqN27iNC1P2iTNpWl6fJo6jtvktGnT5OQkeRLXTWPn5iYkaV03bpPaweAL\ntkHYYC5G5iYQAkYgIQkkoct7/thb40GWYEDaM5rR7/M887Bn77Vn3tmM9jtrrb3WNndHREQEICfd\nAYiIyOihpCAiInFKCiIiEqekICIicUoKIiISp6QgIiJxSgoiIhKnpCAZzcx2m9kNp9h+jpn1mdm3\nB9l2u5m9bGatZtZkZk+a2ZxwW5mZfc/MDphZm5ltN7O/SNjXzOzPzazOzDrMrN7M/tbM8qP4nCKp\noqQg2e6DQDNwZ+IJ28zmA48AfwaUAucA3wb6wiJfA4qA88PttwE7El73G8CK8PWLgZuA64GfjkTQ\nZjYumXVn+hoip6OkINnug8DngW7g1oT1S4Fd7v6kB9rc/efuXh9ufwvwY3dvdvc+d9/m7isBzGwB\n8CfAB9z9eXfvcffNwB3AjWZ2/WCBmFmpmf2zmTWaWYOZ/Y2Z5Ybb7jGzZ83sa2Z2BLh/iHU5ZvZ5\nM9tjZofM7BEzKw1fY46ZuZl9xMzqgadG+mBK9lNSkKxlZlcB1cCjBL/gP5iwuRY4LzzhXmdmRQN2\nXwt8ycw+HCaBRG8D9rn7i4kr3X1vuN/bhwjpYaAHmA8sA94B/GHC9suAncBU4EtDrLsnfFwHzCWo\nzXxzwPtcQ1DDeecQcYgMSUlBstmHgCfcvRn4MXCTmU0FcPedwLXADIKE0WRm309IDp8AfgR8HNhi\nZq+b2U3htilA4xDv2RhuP4mZTSNoYvqUux9z90METVR3JhTb7+7/L6x5dAyx7gPAP7r7TndvB+4l\naBpLbCq6P3yPDkTOkJKCZCUzKwTeR3Bix92fB+qB3+sv4+5r3f397l4BXAVcDfxluK3D3b/s7pcA\nkwkSx8/MbBLQBFQO8daV4faBZgN5QKOZtZhZC/BdghpAv72D7DdwXRWwJ+H5HmAcMO00ryOSFCUF\nyVbvAUqAb4dXEB0gqBV8cLDC7v4S8AtgySDbWoEvAxMJOqSfAmaa2aWJ5cxsJrAceHKQt9gLdAFT\n3L0sfJS4++LEtxostAHP9xMkmH6zCJqkDp7mdUSSoqQg2SDPzAoSHuMImo6+B1xA0Km8FLgCWGpm\nF5jZlWb2R/3NSWZ2HsEVRmvD539lZm8xs/FmVgB8EmgBXnP37cB3gB+Z2XIzyzWzxcDPgf9x9/8Z\nGKC7NwK/Br5qZiVhh/E8M7vmDD/rT4BPh5faFhEkq391954zfB2RQSkpSDb4FdCR8Pgngs7gr7v7\ngYTHeuC/CBJGC0ES2GRm7eH6XwJfCV/TgX8haAraT9B5fHPYjg9BX8NDwA+B/v1XEVyBNJQPAuOB\nLQSXya5k6GaooXwP+AGwGtgFdBL0f4iMCNNNdkREpJ9qCiIiEqekICIicUoKIiISp6QgIiJxGTdh\n1pQpU3zOnDnpDkNEJKOsX7++KRyoeUoZlxTmzJnDunXr0h2GiEhGMbM9py+l5iMREUmgpCAiInFK\nCiIiEqekICIicUoKIiISp6QgIiJxSgoiIhKXceMURETGgs7uXhqPdrK/pYOGlg72t3TwtvOmcUF1\naaTvq6QgIpJi7k5T+wn2hyf74KQfJID9R4N1Te0n3rTflKJ8JQURkUzT1dNLY0snDS0dNDR3sC88\n+ccfRzs50dN30j6FebnMKC+kqqyQxVUlVJUGy1VlhcwoK2RaaT7543Ijj11JQUTkDLV2dtPQ/Mav\n/P4Tf0Nz8DzW1nVSeTOYWpwfnPBnlPLOxdPjJ/yqsgJmlBVSWpiHmaXpE71BSUFEJIG7c+TYCfaF\nJ/h9zcfjJ/v+dW2dJ98Se3xuDlVlBVSVFXLtwgpmlAe/7meUF1JdNoHppQWMH5cZ1/UoKYjImNLX\n5zS1d7Gv/yTfHJ74E553dPeetE9x/rj4if4tcyYNOOkXMqUon5yc9P/KHwlKCiKSVXp6+zjQ2hnv\nuB14wt/X0vGm9vyyCXlUlxcyr2Ii1yysYEZZIdXlb/zSL52Ql6ZPk3pKCiKSUdq7eoK2/LAppyGh\nA7ehuYMDrZ30+cn7TCkaz4zyCZxfWcLbF00LTvblhcwom8CM8kKK8nUq7KcjISKjRmd3L4dauzjQ\n2smB1k4aEy7ZbGjppKH5OK0D2vPH5RiVZQVUlRayfN7koFmn7OQrdwrHR3/VTrZQUhCRyLk7zce7\nOXC0k4PhCX+w5ebj3W/at7hgXPxEXzO7PH7Z5oyyAmaUTaCiOJ/cLGnPHw2UFERkxD2/4zA/emFP\n/KR/sLXrTe34ZjB5Yj7TS/OpLi/kktnlTC8pYFppAdNLCpheGjxKCsZOe/5ooKQgIiPuq79+jW0H\n2lhcVcLFs8KTfXii7/93anE+ebmZcZnmWBJ5UjCzXGAd0ODutwzY9lHgY0Av0A6scPctUcckItFp\n7exmw94W/viaeXz2neemOxw5Q6lI058Etg6x7cfufoG7LwW+AvxjCuIRkQg9v+MwvX3OVQumpDsU\nOQuRJgUzqwZuBh4abLu7tyY8nQj4YOVEJHOs3h5j4vhcls0qT3cochaibj76OvA5oHioAmb2MeAz\nwHjg+iHKrABWAMyaNWvkoxSREbOmrom3zpucMdM6yMki+18zs1uAQ+6+/lTl3P1b7j4P+Avg80OU\nedDda9y9pqKiIoJoRWQk7Dl8jPojx7lqgf5OM1WUqfwK4DYz2w08ClxvZj88RflHgXdHGI+IRGx1\nXRMAVy9UUshUkSUFd7/X3avdfQ5wJ/CUu9+dWMbMFiQ8vRmoiyoeEYnemu0xqssLmTN5QrpDkbOU\n8nEKZvYAsM7dHwM+bmY3AN1AM/ChVMcjIiOju7eP53cc5paLqkbFfQHk7KQkKbj7KmBVuHxfwvpP\npuL9RSR6r+xtoa2rh6t1KWpG0+UBIjIiVtc1kWNw+TwlhUympCAiI2JNXYyLZpaNqXsPZCMlBREZ\ntqPHu3llb4suRc0CSgoiMmzP7miiz1F/QhZQUhCRYVtTF6M4fxxLZ5alOxQZJiUFERkWd2f19iYu\nnz+ZcZoKO+Ppf1BEhmVX0zEaWjrUn5AllBREZFjW9E9toaSQFZQURGRY1tTFmD15ArM0tUVWUFIQ\nkbN2oieY2kI31MkeSgoictZq65s5dqJX/QlZRElBRM7amroYuTnGW+dNTncoMkKUFETkrK2pa+Li\nWWWUFGhqi2yhpCAiZ+XIsRNsajiqpqMso6QgImfl2debcEedzFlGSUFEzsqauhglBeO4sFpTW2QT\nJQUROWPuzpq6Jq5cMIXcHN1lLZsoKYjIGdsRa6fxaKf6E7JQ5EnBzHLNbIOZPT7Its+Y2RYz22hm\nT5rZ7KjjEZHhe3p7MLXFlfPVn5BtUlFT+CSwdYhtG4Aad78QWAl8JQXxiMgwramLMXfKRGZO0tQW\n2SbSpGBm1cDNwEODbXf337r78fDpWqA6ynhEZPi6enpZu1NTW2SrqGsKXwc+B/QlUfYjwBODbTCz\nFWa2zszWxWKxkYxPRM7Q+t3NdHb3cfVC9Sdko8iSgpndAhxy9/VJlL0bqAH+frDt7v6gu9e4e01F\nhb6IIum0uq6JvFxj+VxNbZGNxkX42lcAt5nZu4ACoMTMfujudycWMrMbgL8ErnH3rgjjEZERsKYu\nxsWzypmYH+XpQ9IlspqCu9/r7tXuPge4E3hqkISwDPgucJu7H4oqFhEZGU3tXWze36qmoyyW8nEK\nZvaAmd0WPv17oAj4mZm9bGaPpToeEUnes68Hl6Kqkzl7paT+5+6rgFXh8n0J629IxfuLyMhYvb2J\n8gl5LK4qTXcoEhGNaBaRpARTW8S4Yr6mtshmSgoikpTXDrZxqK2LqzW1RVZTUhCRpKwJp7a4aqH6\nE7KZkoKIJGV1XYwFU4uoLC1MdygSISUFETmtzu5eXtx1RLOijgFKCiJyWi/tPkJXT5+ajsYAJQUR\nOa01dU2Mz83hsnMmpTsUiZiSgoic1urtMWrmlDNhvKa2yHZKCiJySodaO9l2oE39CWOEkoKInNKa\nOk1tMZYoKYjIKa2pizGlaDyLKkvSHYqkgJKCiAypr8955vUmrpw/hRxNbTEmKCmIyJC2Hmilqf2E\n+hPGECUFERmS+hPGHiUFERnSmroY500vZmpJQbpDkRRRUhCRQXWc6OWlXc2qJYwxSgoiMqi1uw5z\nordP/QljTORJwcxyzWyDmT0+yLarzazWzHrM7L1RxyIiyVuzvYnx43K4VFNbjCmpqCl8Etg6xLZ6\n4B7gxymIQ0TOwJq6GJedM4mCvNx0hyIpFGlSMLNq4GbgocG2u/tud98I9EUZh4icmcajHdQdatdd\n1sagqGsKXwc+xzBP+ma2wszWmdm6WCw2MpGJyJDil6JqquwxJ7KkYGa3AIfcff1wX8vdH3T3Gnev\nqajQLxeRqK2pa6KiOJ9zpxWnOxRJsShrClcAt5nZbuBR4Hoz+2GE7yciI6Cvz3mmLsZVC6Zgpqkt\nxprIkoK73+vu1e4+B7gTeMrd747q/URkZGze30rz8W71J4xRKb9jhpk9AKxz98fM7C3AL4Fy4FYz\n+6K7L051TCJjWU9vH6/H2tnc0Mrm/a08tyPoT7hivvoTxqKUJAV3XwWsCpfvS1j/ElCdihhEBDq7\ne9l2oI3N+4/yakMrW/YfZduBNrp6gmtBCvJyWFRZwuduPJeK4vw0RyvpoHvriWSp1s5utuxv5dWG\no2zZH9QCXo+109vnAJQUjGPJjFI++NbZLJlRyuKqEs6ZUkSupsge05QURLLAobZONu9vDU/+R9m8\nv5U9h4/Ht08tzmfJjFLesXgai6uCBFBdXqiOZHkTJQWRDNLX59QfOc7mhJP/5v2tNLV3xcvMmjSB\nJTNKeH/NTBZXlbC4qlRNQZI0JQWRUepETx91h9pOqgFsbWyjvasHgHE5xvypRVyzsCI8+ZdwflUJ\nJQV5aY5cMpmSgsgo0NbZzbYDbSc1/9QdbOdEb9ABXJiXy/mVxbxn2Yz4r/8F04o0L5GMOCUFkRRy\ndw60drIl/PW/pTF4JLb/T5o4nsVVJXz4yjksriplUWUJ50yZqA5gSQklBZGIdPf2sSPW/qYE0HK8\nO15mzuQJLK4q4X2XVHN+ZQmLqkqYXlKgDmBJGyUFkRFwtKObbY2tbG184+S//cAbzT/543I4d3ox\nNy6ezqKqEhZVlnBeZQlF+foTlNFF30iRM9Db5+xqOsa2A61sa2xj24FWtja20dDSES8Tb/65Yk48\nAZwzZSLjcnWjQxn9lBREhtB87ARbE07+2w608VrC6N/cHGNexUQumV3OB5bP4vzpQfPP1OJ8Nf9I\nxlJSkDHvRE9f/Nf/1v4E0NjGgdbOeJnJE8dzfmUJv798NudVlnB+ZTHzpxaRP05X/0h2UVKQMaOn\nt489R45Td7CN1w60s/1QG9sPtLGr6Rg94dQPebnG/KnFXD5vMudVFnN+ZQnnTS/R4C8ZM5QUJOv0\n9Tl7m4+z/WA72w+2hY92dsTaORE2/ZjBzPIJLJxWzNsXTWPhtCABzK2YSJ7a/mUMU1KQjNXb5zQ0\nd7Aj1h4/8W8/2Mbrh9rp6O6Nl5tRVsiCaUVcvWAKC6YVc+60oOmncLyafkQGUlKQUa+9q4edsXZ2\nxo6xIxb84t8ZO8bOpmPxX/4QTPp27vRi7rp0FgunFbFwejELphZRrGkfRJKmpCCjQl9fMNJ3R6yd\nHYfa2RE7xs6mdnYcOnZSh2+OwezJE5k7ZSJXL6xgXsVE5lYUsWBqEWUTxqfxE4hkh9MmBTPLBf7U\n3b+WgnhkDHF3vvSfW3l+52F2xo6d1ORTXDCOeRVFXD5/MvMqisLHRGZNnqArfkQidNqk4O69ZnY7\noKQgI2r/0U4eemYXF1aXctels5g3dSLzKoqYWzGRiiJd6y+SDsk2Hz1rZt8E/hU41r/S3WtPt2NY\n01gHNLj7LQO25QOPAJcAh4HfdffdScYkGa52TzMAX37PBSyZUZrmaEQEkk8Kl4f/PpCwzoHrk9j3\nk8BWoGSQbR8Bmt19vpndCfwd8LtJxiQZrra+mYK8YE4gERkdkkoK7n7d2by4mVUDNwNfAj4zSJHb\ngfvD5ZXAN83M3N3P5v0ks9TWt3BhdZnGBYiMIkn9NZpZqZn9o5mtCx9fNbNk6vtfBz4H9A2xfQaw\nF8Dde4CjwORB3n9F/3vHYrFkQpZRrrO7ly37j3LxrPJ0hyIiCZL9ifY9oA14f/hoBf7lVDuY2S3A\nIXdff6pig6x7Uy3B3R909xp3r6moqEgyZBnNNu8/Snevc/GssnSHIiIJku1TmOfudyQ8/6KZvXya\nfa4AbjOzdwEFQImZ/dDd704osw+YCewzs3FAKXAkyZgkg9XuaQFgmWoKIqNKsjWFDjO7sv+JmV0B\ndJyiPO5+r7tXu/sc4E7gqQEJAeAx4EPh8nvDMupPGANq65uZOalQE82JjDLJ1hQ+CjyS0I/QzBsn\n8zNiZg8A69z9MeCfgR+Y2esENYQ7z+Y1JfNsqG/hsrmT0h2GiAyQzIjmHOBcd7/IzEoA3L31TN7E\n3VcBq8Ll+xLWdwLvO5PXksy3v6WDA62d6mQWGYVO23zk7n3Ax8Pl1jNNCCID1dYHg9aWqZNZZNRJ\ntk/hN2b2WTObaWaT+h+RRiZZq3ZPCwV5OZxfOdh4RhFJp2T7FP4g/PdjCescmDuy4chYsGFvMxfO\n0KA1kdEo2T6Fu9392RTEI1muq6eXzQ2tfPjKOekORUQGkWyfwj+kIBYZA15taOVEbx/LZqqTWWQ0\nSrb+/mszu8M0l7EM04awk/ni2epkFhmNku1T+AwwAeg1s06C6Snc3dVTKGektr6Z6vJCphYXpDsU\nERlEskmhFPgAcI67P2Bms4DK6MKSbLWhvoW3zNGFayKjVbLNR98ClgN3hc/bgG9GEpFkrcajHTQe\n7dT4BJFRLNmawmXufrGZbQBw92Yz013S5Yz0T4Knkcwio1eyNYXu8LaaDmBmFQx9jwSRQdXWN5M/\nToPWREazZJPCN4BfAlPN7EvAM8CXI4tKstKG+mYurC5l/DgNWhMZrZK9HeePzGw98DaCK4/e7e5b\nI41MskpXTy+vNrRyzxVz0h2KiJxCsn0KuPs2YFuEsUgW27w/GLSmO62JjG6qx0tK1O4JB62pk1lk\nVFNSkJTYsLeFGWWFTC3RoDWR0UxJQVJiw55mjU8QyQBKChK5A0c72X9Ud1oTyQSRJQUzKzCzF83s\nFTPbbGZfHKTMbDN70sw2mtkqM6uOKh5Jn9r4JHhKCiKjXZQ1hS7gene/CFgK3GhmyweU+QfgEXe/\nEHgA+D8RxiNpsqG+mfHjclikQWsio15kScED7eHTvPDhA4otAp4Ml38L3B5VPJI+tfUtXDBDg9ZE\nMkGkf6VmlmtmLwOHgN+4+wsDirwC3BEuvwcoNrPJg7zOCjNbZ2brYrFYlCHLCDvR08emhqManyCS\nISJNCu7e6+5LgWrgUjNbMqDIZ4Frwon2rgEagJ5BXudBd69x95qKioooQ5YRtnn/UU709KmTWSRD\nJD2ieTjcvcXMVgE3Aq8mrN8P/A6AmRUBd7j70VTEJKmxoT6cGVWdzCIZIcqrjyrMrCxcLgRuYMA0\nGWY2xcz6Y7gX+F5U8Uh61NY3U1VawDQNWhPJCFE2H1UCvzWzjcBLBH0Kj5vZA2Z2W1jmWuA1M9sO\nTAO+FGE8kgYb6ltYplqCSMaIrPnI3TcCywZZf1/C8kpgZVQxSHodbO2koaWDP7jynHSHIiJJ0jWC\nEpk3JsHTlUcimUJJQSKzYW8L43NzWFSlQWsimUJJQSJTu6eZJTNKyB+Xm+5QRCRJSgoSiRM9fWxs\nOKrxCSIZRklBIrGlsTUYtKYrj0QyipKCRGJDODOq7qEgklmUFCQStfUtVJYWUFlamO5QROQMKClI\nJGr3NKs/QSQDKSnIiDsUDlpT05FI5lFSkBFXG06Ct0w1BZGMo6QgI25DfTPjc3NYMkOD1kQyjZKC\njLja+mYWa9CaSEZSUpARdaKnj437NGhNJFMpKciI2nagla6ePnUyi2QoJQUZUW/MjKqagkgmUlKQ\nEVVb38L0kgKqyjRoTSQTKSnIiKqtb+bi2Wo6EslUUd6jucDMXjSzV8xss5l9cZAys8zst2a2wcw2\nmtm7oopHoneorZN9zR0sm6mmI5FMFWVNoQu43t0vApYCN5rZ8gFlPg/81N2XAXcC344wHonYhnDQ\nmmoKIpkryns0O9AePs0LHz6wGNA/wqkU2B9VPBK92vpm8nKNxVWl6Q5FRM5SpH0KZpZrZi8Dh4Df\nuPsLA4rcD9xtZvuAXwGfGOJ1VpjZOjNbF4vFogxZhmHDnhYWV5VSkKdBayKZKtKk4O697r4UqAYu\nNbMlA4rcBXzf3auBdwE/MLM3xeTuD7p7jbvXVFRURBmynKXu3j42NrRofIJIhkvJ1Ufu3gKsAm4c\nsOkjwE/DMs8DBcCUVMQkI2tbYxud3X0anyCS4aK8+qjCzMrC5ULgBmDbgGL1wNvCMucTJAW1D2Wg\n2vBOa7r9pkhmi6yjGagEHjazXILk81N3f9zMHgDWuftjwJ8B/2RmnybodL4n7KCWDFNb38y0knyq\nSgvSHYqIDEOUVx9tBJYNsv6+hOUtwBVRxSCpU1vfzLKZ5ZhZukMRkWHQiGYZtlhbF3uPdGh8gkgW\nUFKQYdtQr0nwRLKFkoIMW219C3m5xpIZGrQmkumUFGTYauubWVRZokFrIllASUGGpae3j437Wlim\npiORrKCkIMOy7UA4aE3jE0SygpKCDEt80JqmtxDJCkoKMiy1e5qpKM5nhu60JpIVlBRkWDbsbeHi\nWWUatCaSJZQU5Kw1tXex5/BxjU8QySJKCnLW3rjTmpKCSLZQUpCzVlvfzLgc4wINWhPJGkoKclZ6\n+5wXdx1hUZUGrYlkEyUFOSPdvX38fP0+3v61p1m/p5lrz52a7pBEZARFeT8FySJdPb2sXL+P7zy9\ng71HOjhvejHf/L1l3LSkMt2hicgIUlKQUzp+ooefvLiXB1fv4GBrF0tnlnH/rYu5/rypugxVJAsp\nKcig2jq7eeT5PfzzM7s4cuwEy+dO4qvvW8oV8ycrGYhkMSUFOUnzsRP8y7O7+P5zu2nt7OHacyv4\n+HXzqZkzKd2hiUgKRJYUzKwAWA3kh++z0t2/MKDM14DrwqcTgKnurkl00uBQWycPrdnFD9fu4fiJ\nXm5cPJ2PXTefC6p1uanIWBJlTaELuN7d280sD3jGzJ5w97X9Bdz90/3LZvYJBrmns0SroaWD7z69\ng0df2ktPbx+3XVTFn1w3n4XTitMdmoikQWRJwd0daA+f5oUPP8UudwFfOMV2GUE7Y+185+kd/KK2\nATO44+JqPnrNPOZMmZju0EQkjSLtUzCzXGA9MB/4lru/MES52cA5wFNDbF8BrACYNWtWNMGOAX19\nzuq6GA8/t5vfvhYjf1wOdy+fzYqr51KlWU5FhIiTgrv3AkvNrAz4pZktcfdXByl6J0GfQ+8Qr/Mg\n8CBATU3NqWobMoi2zm5+vn4fjzy/h51Nx6gozudTNyzgA5fNpqI4P93hicgokpKrj9y9xcxWATcC\nQyWFj6UilrFkV9MxHn5uNyvX76O9q4elM8v4v3cu5aYllYwfp8HsIvJmUV59VAF0hwmhELgB+LtB\nyp0LlAPPRxXLWNLX5zwdNhGtei1GXq5xy4VVfOjyOSydqQu7ROTUoqwpVAIPh/0KOcBP3f1xM3sA\nWOfuj4Xl7gIeDTum5Sz1NxE9/PwedoVNRJ++YSF3XTaTqcUF6Q5PRDJElFcfbWSQS0zd/b4Bz++P\nKoaxYGesnUee38PP1u3l2Ilels1SE5GInD2NaM5A/U1E3392N09vD5qIbg2biC5SE5GIDIOSQgbZ\n13ycn69vYGXtXvYe6WCqmohEZIQpKYxyHSd6+e/NB/jZ+r08t+Mw7nD5vMl89h3nqolIREacksIo\n5O7U1rewcv1eHn+lkbauHmZOKuRTb1vIHZfMoLp8QrpDFJEspaQwihxs7eQXtQ2sXL+XHbFjFObl\nctMF03nfJTO57JxJ5ORoymoRiZaSQpp19fTy5NZD/GzdXp7eHqPPoWZ2OX93x1xuvrCKonz9F4lI\n6uiMkwbuzub9raxcv49/e7mBluPdTC8p4KPXzOO9l1Qzt6Io3SGKyBilpJBCe48c54lXG/lFbQPb\nDrQxflwO71g0jfdeUs1VCyrIVfOQiKSZkkLEdjcd41evNvLEpgNsajgKwIXVpfz17Yu59aIqyiaM\nT3OEIiJvUFKIwI5YO09sauQ/Nx1ga2MrABfNLOPem87jpiWVzJqsq4dEZHRSUhgB7k7doXZ+tSmo\nEbx2sA2AS2aX8/mbz+fGJdN1GamIZAQlhbPk7mxtbOOJVxv51aZGdsSOYQZvmTOJ+29dxI1LKple\nqlHGIpJZlBTOQP9VQ/+5qZEnNjWy+/BxcgwuO2cy91xxDu9cPE3TTYhIRlNSOAV3Z0fsGGt3Hg4f\nR2hq7yI3x7h83mRWXD2PdyyexpQi3b1MRLKDkkKCoZIAwPSSAq6cP5nL50/h7edPo3yirhoSkewz\nppNCkATaeX7nEdbuPMwLA5LAVQumsHzuJJbPncysSRMw0zgCEcluYyopvDkJHKap/QSgJCAiAtHe\no7kAWA3kh++z0t2/MEi59wP3Aw684u6/F0U8j75Yzz/8+rV4EqgsLeCqBRVKAiIiCaKsKXQB17t7\nu5nlAc+Y2RPuvra/gJktAO4FrnD3ZjObGlUw00oKuHpBBcvnTmb53MnMnFSoJCAiMkCU92h2oD18\nmhc+fECxPwK+5e7N4T6HoornuvOmct15keUcEZGsEOltu8ws18xeBg4Bv3H3FwYUWQgsNLNnzWyt\nmd0YZTwiInJqkSYFd+9196WMdv4TAAAHYElEQVRANXCpmS0ZUGQcsAC4FrgLeMjM3nTneTNbYWbr\nzGxdLBaLMmQRkTEtJTf4dfcWYBUwsCawD/h3d+92913AawRJYuD+D7p7jbvXVFRURB6viMhYFVlS\nMLOK/l/9ZlYI3ABsG1Ds34DrwjJTCJqTdkYVk4iInFqUVx9VAg+bWS5B8vmpuz9uZg8A69z9MeC/\ngXeY2RagF/hzdz8cYUwiInIKFlwklDlqamp83bp16Q5DRCSjmNl6d685XbmU9CmIiEhmUFIQEZG4\njGs+MrMYsCfdcZzGFKAp3UEkQXGOrEyJEzInVsU5cma7+2kv38y4pJAJzGxdMm136aY4R1amxAmZ\nE6viTD01H4mISJySgoiIxCkpROPBdAeQJMU5sjIlTsicWBVniqlPQURE4lRTEBGROCUFERGJU1I4\nA2Z2o5m9Zmavm9n/GmT7Z8xsi5ltNLMnzWx2wrZeM3s5fDyW5jjvMbNYQjx/mLDtQ2ZWFz4+FGWc\nScb6tYQ4t5tZS8K2VB7T75nZITN7dYjtZmbfCD/HRjO7OGFbyo5pEnF+IIxvo5k9Z2YXJWzbbWab\nwuMZ6VwyScR5rZkdTfj/vS9h2ym/MymO888TYnw1/E5OCrel7HiOKHfXI4kHkAvsAOYC44FXgEUD\nylwHTAiX/xj414Rt7aMoznuAbw6y7ySCWWonAeXhcnk6Yx1Q/hPA91J9TMP3uhq4GHh1iO3vAp4A\nDFgOvJCmY3q6OC/vf3/gpv44w+e7gSmj5HheCzw+3O9M1HEOKHsr8FQ6judIPlRTSN6lwOvuvtPd\nTwCPArcnFnD337r78fDpWoKbC6XaaeM8hXcS3CHviAe3SP0Nb74Hxkg601jvAn4SYTxDcvfVwJFT\nFLkdeMQDa4EyM6skxcf0dHG6+3NhHJC+72gyx3Mow/l+n7EzjDNt38+RpKSQvBnA3oTn+8J1Q/kI\nwS/HfgXh3ePWmtm7owgwlGycd4RNCCvNbOYZ7jtSkn6/sCnuHOCphNWpOqbJGOqzpPqYnomB31EH\nfm1m681sRZpiSvRWM3vFzJ4ws8XhulF5PM1sAkGy/3nC6tF2PJMS5f0Uso0Nsm7Q63nN7G6gBrgm\nYfUsd99vZnOBp8xsk7vvSFOc/wH8xN27zOyjwMPA9UnuO5LO5P3uBFa6e2/CulQd02QM9VlSfUyT\nYmbXESSFKxNWXxEez6nAb8xsW/hLOR1qCebqaTezdxHckGsBo/R4EjQdPevuibWK0XQ8k6aaQvL2\nATMTnlcD+wcWMrMbgL8EbnP3rv717r4//Hcnwa1Jl6UrTnc/nBDbPwGXJLvvCDuT97uTAVXzFB7T\nZAz1WVJ9TE/LzC4EHgJu94SbWiUcz0PALwmaatLC3VvdvT1c/hWQZ8HdGUfd8Qyd6vuZ9uN5RtLd\nqZEpD4Ja1U6CJoz+Dq7FA8osI+gEWzBgfTmQHy5PAeqIqHMsyTgrE5bfA6wNlycBu8J4y8PlSek8\npmG5cwk67SwdxzThPecwdMfozZzc0fxiOo5pEnHOAl4HLh+wfiJQnLD8HHBjGuOc3v//TXAyrQ+P\nbVLfmVTFGW4vJeh3mJjO4zlSDzUfJcnde8zs4wS3EM0luApms518e9G/B4qAn5kZQL273wacD3zX\nzPoIamd/6+5b0hjnn5rZbUAPwZf5nnDfI2b218BL4cs94CdXh9MRKwQdeI96+BcWStkxBTCznxBc\nETPFzPYBXwDyws/xHeBXBFcgvQ4cBz4cbkvpMU0izvuAycC3w+9ojweze04DfhmuGwf82N3/K41x\nvhf4YzPrATqAO8P//0G/M2mME4IfVr9292MJu6b0eI4kTXMhIiJx6lMQEZE4JQUREYlTUhARkTgl\nBRERiVNSEBGROCUFERGJU1IQEZE4DV4TAczsr4APEEy21gSsB44CKwhGzr4O/L67Hzez7xMMqDoP\nmE0wUO1DwFsJpqK+J3zNduBbwA1AM/C/ga8QjCr+lLs/ZmZzgB8QjHoF+Li7PxftpxUZmmoKMuaZ\nWQ1wB8E0Jb9DMJkhwC/c/S3ufhGwlWACuX7lBJMIfppggsGvAYuBC8xsaVhmIrDK3S8B2oC/Ad5O\nMAL2gbDMIeDt7n4x8LvANyL5kCJJUk1BJJgp9N/dvQPAzP4jXL/EzP4GKCOYvuS/E/b5D3d3M9sE\nHHT3TeG+mwnmynkZOAH0T22wCehy9+5wnznh+jzgm2Ei6QUWRvMRRZKjpCAy+HTMAN8H3u3ur5jZ\nPQRz4PTrn2W2L2G5/3n/31V3wnxN8XLu3mdm/WU+DRwELiKouXee9acQGQFqPhKBZ4BbzazAzIoI\nZjwFKAYazSyPoL8hCqVAo7v3Ab9PMMmbSNqopiBjnru/ZGaPEUzDvAdYR9DJ/FfAC+G6TQRJYqR9\nG/i5mb0P+C1w7DTlRSKlWVJFADMr8uAuXxOA1cAKd69Nd1wiqaaagkjgQTNbBBQADyshyFilmoKI\niMSpo1lEROKUFEREJE5JQURE4pQUREQkTklBRETi/j9HVsGCvy8wzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f4c5f98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X+UXWV97/H3Z37k5xwgIckZCQkB\nyZkKahCnIEIVrNJAVerVtnCRgpVmwdKqver1R1ehRW176622FimNgqhV1BaoUUFhVZECFyRQIEDM\nD0MwIYEJCYT8TibzvX/sfejpOD92krPPnnPm81prr+zz7Gfv8z17TuY7ez/Pfh5FBGZmZqNpKzoA\nMzNrDk4YZmaWiROGmZll4oRhZmaZOGGYmVkmThhmZpaJE4a1NEk3SPp0Qe8tSV+R9LyknxURg1k9\nOWFYQ0laK+lZSVNryi6VdGeBYeXlDOAtwNERcUrRwZgdKicMK0IH8MGigzhQktoPcJdjgLURsSOP\neBpNUkfRMVixnDCsCJ8FPiLpiMEbJM2TFLW/nCTdKenSdP0SSfdI+rykFyStkfT6tHydpD5JFw86\n7AxJd0jaJumnko6pOfavpdu2SFoh6fdqtt0g6R8l3SppB3DWEPEeJWlJuv9qSX+Ulr8X+DJwmqTt\nkv5iiH0vkXS3pP+b3rZ6UtI5ox073fbnkr4j6Wvp53pcUu9wJzw9V9vTZUd6juel294q6eG0zr2S\nXl2z31pJH5P0KLBDUoekV6Q/kxfS9317Tf1zJT2RxvS0pI8MF5M1oYjw4qVhC7AWeDNwM/DptOxS\n4M50fR4QQEfNPncCl6brlwD9wHuAduDTwC+BLwITgbOBbUBXWv+G9PUb0u1/D9ydbpsKrEuP1QGc\nDDwHnFiz71bgdJI/riYN8Xl+ClwDTAJOAjYBv1kT690jnItLgH3AH6Wf5XJgA6AMx/5zYDdwbrrv\nXwH3ZfwZ/CVwF9CZfuY+4NT0OBenP6OJNT+vh4E5wOR0n9XAJ4EJwJvS89uT1t8I/Ea6Pg04uejv\nnJf6Lb7CsKJcAfyxpJkHse+TEfGViNgPfJvkl9lVEbEnIm4H9gLH19T/QUTcFRF7gD8l+at/DvBW\nkltGX4mI/oh4CLgJeFfNvt+NiHsiYiAidtcGkR7jDOBjEbE7Ih4muaq46AA+y1MR8aX0s3wVeBlQ\nznjsuyPi1nTfrwMLRnszSb8P/E/gnRFRTVb/FBH3R8T+iPgqsAd4Xc1uX4iIdRGxKy3vAv46IvZG\nxI+B7wMXpHX3ASdIOiwink/PqbUIJwwrREQ8RvKL5uMHsfuzNeu70uMNLuuqeb2u5n23A1uAo0ja\nGE5Nb628IOkF4EKge6h9h3AUsCUittWUPQXMPoDP8kxNbDvT1a6Mx36mZn0nMCm9ZXRhze2n26oV\nJL0GuBp4R0RsSouPAT486BzMSd+/qvYcHAWsi4iBYeJ6J8lVz1Pp7b/TspwEaw5uxLIiXQk8BPxt\nTVm1gXgK8GK6XvsL/GDMqa5I6gKmk9z6WQf8NCLeMsK+Iw3nvAGYLqlU84t9LvD0IcZ7SMeOiG8A\n36gtS6/kbgHeHxH/WbNpHfCZiPjMSIccFNccSW01SWMusDJ97weA8yR1Au8HvkPN+bfm5isMK0xE\nrCa5pfSBmrJNJL8U3y2pXdIfAi8/xLc6V9IZkiYAnwLuj4h1JFc4FUkXSepMl1+X9IqM8a8D7gX+\nStKktLH4vQz6ZX0w6nnstAPBTcA3IuLbgzZ/CbhM0qlKTJX025JKwxzufpKk/r/T83Um8DbgW5Im\npFc3h6e3u14E9h9ovDZ2OWFY0a4iaXyu9UfAR4HNwIkkvzgPxTdJrma2AK8lue1E+pf72cD5JH85\nPwP8H5LG8awuIGmo30DyF/yVEXHHIcZb72MfDfwG8KGaW1XbJc2NiKUk5/tq4HmSBu1LhjtQROwF\n3g6cQ9JB4BrgDyLi52mVi4C1kl4ELgPefRDx2hhV7Y1hZmY2Il9hmJlZJk4YZmaWiROGmZll4oRh\nZmaZtNRzGDNmzIh58+YVHYaZWdN48MEHn4uITCMutFTCmDdvHkuXLi06DDOzpiHpqax1fUvKzMwy\nccIwM7NMnDDMzCwTJwwzM8vECcPMzDLJLWGkI2z+TNIj6TSOQ01ROVHSt9PpJ++vThmZbvtEWr5C\n0m/lFaeZmWWT5xXGHuBNEbGAZHrJhZJeN6jOe4HnI+J44PMkI4Ui6QSSEURPBBYC10hqzzFWMzMb\nRW7PYUQyDO729GVnugweGvc8krmJAf4VuFqS0vJvpVNqPilpNXAK8P/yiPUL/76K/v0Do1e0zGYd\nNokLT51L8uM0s1aQ64N76VXBgyTzK38xIu4fVGU26fSPEdEvaStwZFp+X0299Qwz7aWkRcAigLlz\n5x5UnNf+9Bfs2ud5XuqlOmL+GcfPYN6MwVNdmFmzyjVhpJPTnyTpCOAWSa9M53KuGurPzxihfKj3\nWAwsBujt7T2oyT2euGrhwexmw3hk3Quc98V7WPHsNicMsxbSkF5SEfECcCdJe0St9aTz/abTSB5O\nMivaS+Wpo0lmHbMmcPysLgBWPrNtlJpm1kzy7CU1M72yQNJk4M3AzwdVWwJcnK6/C/hx2vaxBDg/\n7UV1LDAf+FlesVp9TZ3YwZzpk1nxrBOGWSvJ85bUy4Cvpu0YbcB3IuL7kq4ClkbEEuA64Otpo/YW\nkp5RRMTjkr4DPAH0A+9Lb29Zk+gpl1jphGHWUvLsJfUo8Johyq+oWd8N/O4w+38G+Exe8Vm+KuUS\nd67YxN7+ASZ0+PlQs1bg/8mWi57uEv0DwdrNO4oOxczqxAnDcjF/VgmAFW74NmsZThiWi+NmTqW9\nTW7HMGshThiWi0md7cw7coqvMMxaiBOG5aanu8Sqvu2jVzSzpuCEYbmplEus3byD3R52xawlOGFY\nbirlEhGw2lcZZi3BCcNyUym7p5RZK3HCsNzMO3IKE9rbWNnnhGHWCpwwLDcd7W28fFaXByE0axFO\nGJarnnIXK591G4ZZK3DCsFzNL5d4+oVdbNu9r+hQzOwQOWFYrnrShm9fZZg1PycMy1VPd5IwVnmI\nELOm54RhuZp9xGSmTGj3ZEpmLcAJw3LV1ibmz+ryIIRmLSDPKVrnSPqJpOWSHpf0wSHqfFTSw+ny\nmKT9kqan29ZKWpZuW5pXnJa/SrnEimfchmHW7PK8wugHPhwRrwBeB7xP0gm1FSLisxFxUkScBHwC\n+GlEbKmpcla6vTfHOC1nPd0lntu+hy079hYdipkdgtwSRkRsjIiH0vVtwHJg9gi7XADcmFc8VpzK\nSz2lfFvKrJk1pA1D0jyS+b3vH2b7FGAhcFNNcQC3S3pQ0qIRjr1I0lJJSzdt2lS/oK1uqj2lnDDM\nmlvuCUNSF0ki+FBEvDhMtbcB9wy6HXV6RJwMnENyO+sNQ+0YEYsjojciemfOnFnX2K0+ZpUmctik\nDg9CaNbkck0YkjpJksU3IuLmEaqez6DbURGxIf23D7gFOCWvOC1fkujpLvkKw6zJ5dlLSsB1wPKI\n+NwI9Q4H3gh8t6ZsqqRSdR04G3gsr1gtf0lPqW1ERNGhmNlByvMK43TgIuBNNV1nz5V0maTLauq9\nA7g9InbUlJWBuyU9AvwM+EFE/DDHWC1nPd0lXtzdT9+2PUWHYmYHqSOvA0fE3YAy1LsBuGFQ2Rpg\nQS6BWSFqJ1MqHzap4GjM7GD4SW9rCHetNWt+ThjWENOnTmBG10T3lDJrYk4Y1jA93R5TyqyZOWFY\nw1TKJVb1bWdgwD2lzJqRE4Y1TE+5xM69+3n6hV1Fh2JmB8EJwxpmfk1PKTNrPk4Y1jCVcheAJ1My\na1JOGNYwpUmdzD5ishu+zZqUE4Y1VKXcxcpnPZmSWTNywrCGqnSX+EXfdvr3DxQdipkdICcMa6jK\nrBJ79w+wdvPOokMxswPkhGEN5cmUzJqXE4Y11PGzupCcMMyakROGNdSkznbmHTnVCcOsCTlhWMPN\nn9Xlh/fMmpAThjVcT3eJtZt3snvf/qJDMbMDkOcUrXMk/UTSckmPS/rgEHXOlLS1Zka+K2q2LZS0\nQtJqSR/PK05rvEq5xP6BYM2mHaNXNrMxI7cZ94B+4MMR8VA6P/eDku6IiCcG1fuPiHhrbYGkduCL\nwFuA9cADkpYMsa81oWpPqVV92zjhqMMKjsbMssrtCiMiNkbEQ+n6NmA5MDvj7qcAqyNiTUTsBb4F\nnJdPpNZo846cSme73I5h1mQa0oYhaR7wGuD+ITafJukRSbdJOjEtmw2sq6mznmGSjaRFkpZKWrpp\n06Y6Rm15mdDRxrEz3FPKrNnknjAkdQE3AR+KiBcHbX4IOCYiFgD/APxbdbchDjXkrDsRsTgieiOi\nd+bMmfUK23JWKZc8aq1Zk8k1YUjqJEkW34iImwdvj4gXI2J7un4r0ClpBskVxZyaqkcDG/KM1Rqr\np1xi3ZZd7NjTX3QoZpZRnr2kBFwHLI+Izw1Tpzuth6RT0ng2Aw8A8yUdK2kCcD6wJK9YrfEqacP3\n6j6PXGvWLPLsJXU6cBGwTNLDadkngbkAEXEt8C7gckn9wC7g/IgIoF/S+4EfAe3A9RHxeI6xWoP1\nVGffe3YbC+YcUXA0ZpZFbgkjIu5m6LaI2jpXA1cPs+1W4NYcQrMxYM70KUzsaGOle0qZNQ0/6W2F\naG8T88tdbvg2ayJOGFaYSrnkrrVmTcQJwwrTUy7x7It72LpzX9GhmFkGThhWmGpPqZV9vsowawZO\nGFaYSrWnlBu+zZqCE4YV5qjDJ9E1scPtGGZNwgnDCiOJStmTKZk1CycMK1RPd9JTKnle08zGMicM\nK1SlXOL5nft4bvveokMxs1E4YVihqg3fbscwG/ucMKxQ7ill1jycMKxQM7omMH3qBF9hmDUBJwwr\nVLWnlBOG2djnhGGFS8aU2u6eUmZjnBOGFa5SLrF9Tz8btu4uOhQzG0GeM+7NkfQTScslPS7pg0PU\nuVDSo+lyr6QFNdvWSlom6WFJS/OK04rXUx1Tyg3fZmNanlcY/cCHI+IVwOuA90k6YVCdJ4E3RsSr\ngU8BiwdtPysiToqI3hzjtIJVZrlrrVkzyHPGvY3AxnR9m6TlwGzgiZo699bsch9wdF7x2Nh1+JRO\nug+b5MmUzMa4hrRhSJoHvAa4f4Rq7wVuq3kdwO2SHpS0aIRjL5K0VNLSTZs21SNcK8B895QyG/Ny\nTxiSuoCbgA9FxIvD1DmLJGF8rKb49Ig4GTiH5HbWG4baNyIWR0RvRPTOnDmzztFbo/SUS6x6djv7\nB9xTymysyjVhSOokSRbfiIibh6nzauDLwHkRsblaHhEb0n/7gFuAU/KM1YpV6S6xp3+AX27ZWXQo\nZjaMTAlD0sslTUzXz5T0AUlHjLKPgOuA5RHxuWHqzAVuBi6KiJU15VMllarrwNnAY1litebU4zGl\nzMa8rFcYNwH7JR1PkgSOBb45yj6nAxcBb0q7xj4s6VxJl0m6LK1zBXAkcM2g7rNl4G5JjwA/A34Q\nET88gM9lTWZ+uQtw11qzsSxrL6mBiOiX9A7g7yLiHyT950g7RMTdgEapcylw6RDla4AFv7qHtaop\nEzqYM32ye0qZjWFZrzD2SboAuBj4flrWmU9INl71lEu+JWU2hmVNGO8BTgM+ExFPSjoW+Of8wrLx\nqFIusWbTDvb2DxQdipkNIVPCiIgngI8AyyS9ElgfEX+da2Q27vR0l+gfCNZu3lF0KGY2hKy9pM4E\nVgFfBK4BVg73XITZwfJkSmZjW9ZG778Fzo6IFQCSKsCNwGvzCszGn+NmTqW9TW7HMBujsrZhdFaT\nBUD6zIQbva2uJna0M+/IKb7CMBujsl5hLJV0HfD19PWFwIP5hGTjWU93iSc2DDmCjJkVLOsVxuXA\n48AHgA+SjDh72Yh7mB2ESrnEU1t2snvf/qJDMbNBRr3CkNQOXBcR7waGHOLDrF56yiUiYHXfdl45\n+/CiwzGzGqNeYUTEfmCmpAkNiMfGufnuKWU2ZmVtw1gL3CNpCfBSJ/nhBhU0O1jzjpzChPY295Qy\nG4OyJowN6dIGlPILx8a7jvY2Xj6ry2NKmY1BmRJGRPwFJEONR4Qfw7Vc9ZS7eGDt80WHYWaDZH3S\n+zRJTwDL09cLJF2Ta2Q2bs0vl3j6hV1s272v6FDMrEbWbrV/B/wWsBkgIh4BPDSI5eK/JlPaXnAk\nZlYr8xStEbFuUJE7ylsuero9+57ZWJQ1YayT9HogJE2Q9BHS21PDkTRH0k8kLZf0uKQPDlFHkr4g\nabWkRyWdXLPtYkmr0uXiA/pU1tRmHzGZKRPa3bXWbIzJ2kvqMuDvgdnAeuB24H2j7NMPfDgiHkrn\n535Q0h3pUOlV5wDz0+VU4B+BUyVNB64EeoFI910SEW4JHQfa2sT8colVfU4YZmNJ1oQxOSIurC2Q\n1D3SDhGxEdiYrm+TtJwk4dQmjPOAr0VEAPdJOkLSy4AzgTsiYkv6XncAC0lGyLVxoDKrix8s28jl\n/+why8aS9jbxgd+c/9JQ9Da+ZE0YT0r6F+API2JXWnYrcPII+7xE0jzgNcD9gzbNBmrbRtanZcOV\nD3XsRcAigLlz52YJx5rAWxccxbKnt/KLTW74HktW921n9rTJfOKcVxQdihUga8JYBvwHcLek34uI\nXwDKsqOkLuAm4EMRMXgY0qGOESOU/2phxGJgMUBvb++Qdaz5vLEykzdWZhYdhg2y8O/uYpV7r41b\nWRu9IyKuIRmt9nuS3sYwv8BrSeokSRbfiIibh6iyHphT8/pokifKhys3swL1dJfcGWEcy5owBBAR\n9wC/CXwU+LURd5AEXAcsH2HMqSXAH6S9pV4HbE3bPn4EnC1pmqRpwNlpmZkVqOKHKse1rLekzq2u\nRMRGSW8CXj/KPqcDFwHLJD2cln0SmJse51qSdpBzgdXATuA96bYtkj4FPJDud1W1AdzMilNt7F7V\nt52T504rOBprtKxjSW2srkv6fkS8FbhrlH3uZpR2jrR31JDdcyPieuD6LPGZWWO89BT+M9ucMMah\nzE961xiyt5KZtb6jp01mcme7h20Zpw4mYfxn3aMws6bQ1iYq5S4P2zJOHXDCiIg/zCMQM2sO88sl\nz1cyTmUd3vx0SXdIWilpjaQnJa3JOzgzG3t6yiU2bdvD8zv2Fh2KNVjWXlLXAX8CPIhHqTUb1yo1\nowmfetyRBUdjjZT1ltTWiLgtIvoiYnN1yTUyMxuT/mu+Et+WGm+yXmH8RNJngZuBPdXCiHgol6jM\nbMwqHzaRwyZ1uB1jHMqaME5N/+2tKQvgTfUNx8zGOklUyiVWPuOuteNN1gf3zso7EDNrHpXuEj94\ndCMRQTIKkI0HWXtJHS7pc5KWpsvfSjo87+DMbGzqKZfYumsfm7btGb2ytYysjd7XA9uA30uXF4Gv\n5BWUmY1t1TGl3I4xvmRNGC+PiCsjYk26/AVwXJ6BmdnYVSl3AXio83Ema8LYJemM6gtJpwO7Rqhv\nZi3syK6JzOia4K6140zWXlKXA1+tabd4Hrg4n5DMrBlUyiVWeBDCcSVrwlgO/A3wcuAIYCvwO8Cj\nOcVlZmNcpVziX5auY2AgaGtzT6nxIOstqe8CbwN2A08D24EdeQVlZmNfT3eJHXv38/QLvjs9XmS9\nwjg6IhYeyIElXQ+8FeiLiFcOsf2jwIU1cbwCmJnOtreWpFfWfqA/InoH729mxao2fK98dhtzpk8p\nOBprhKxXGPdKetUBHvsGYNgkExGfjYiTIuIk4BPATwdNw3pWut3JwmwMmu+uteNO1iuMM4BLJD1J\nMpaUSGZYffVwO0TEXZLmZTz+BcCNGeua2Rhw2KROjjp8EivdtXbcyJowzskrAElTSK5E3l9THMDt\nkgL4p4hYPML+i4BFAHPnzs0rTDMbQqW75Olax5GsY0k9lWMMbwPuGXQ76vSI2CBpFnCHpJ9HxF3D\nxLYYWAzQ29sbOcZpZoP0lEvc+4vN9O8foKP9YGZ8tmYyFn7C5zPodlREbEj/7QNuAU4pIC4zG8X8\ncom9/QM8tWVn0aFYAxSaMNIHAd9I0m23WjZVUqm6DpwNPFZMhGY2kpcmU3I7xriQtQ3jgEm6ETgT\nmCFpPXAl0AkQEdem1d4B3B4Rtc90lIFb0iGTO4BvRsQP84rTzA7e8bO6kJKeUue86mVFh2M5yy1h\nRMQFGercQNL9trZsDbAgn6jMrJ4mT2jnmOlTWOWG73FhLLRhmFkTS8aU8i2p8cAJw8wOSaVc4snn\ndrCnf3/RoVjOnDDM7JBUukvsHwjWbPLwcq3OCcPMDslLPaV8W6rlOWGY2SE5dsZUOtrkhDEOOGGY\n2SGZ0NHGcTOnsuIZ95RqdU4YZnbI5pdLvsIYB5wwzOyQ9ZRL/HLLTnbu7S86FMuRE4aZHbJK2vC9\nus+3pVqZE4aZHbKe7nQyJY8p1dKcMMzskM2dPoWJHW1ux2hxThhmdsja28Txs7pY4TGlWpoThpnV\nRU+55GHOW5wThpnVRaW7xDMv7mbrrn1Fh2I5ccIws7qoDhGyyu0YLSu3hCHpekl9koacLU/SmZK2\nSno4Xa6o2bZQ0gpJqyV9PK8Yzax+5pe7ADzUeQvL8wrjBmDhKHX+IyJOSperACS1A18EzgFOAC6Q\ndEKOcZpZHcw+YjJTJ7S7HaOF5ZYwIuIuYMtB7HoKsDoi1kTEXuBbwHl1Dc7M6k4SlW5PptTKim7D\nOE3SI5Juk3RiWjYbWFdTZ31aNiRJiyQtlbR006ZNecZqZqPoKZc8XWsLKzJhPAQcExELgH8A/i0t\n1xB1Y7iDRMTiiOiNiN6ZM2fmEKaZZVUpl9i8Yy/Pbd9TdCiWg8ISRkS8GBHb0/VbgU5JM0iuKObU\nVD0a2FBAiGZ2gKpjSrkdozUVljAkdUtSun5KGstm4AFgvqRjJU0AzgeWFBWnmWVX6XZPqVbWkdeB\nJd0InAnMkLQeuBLoBIiIa4F3AZdL6gd2AedHRAD9kt4P/AhoB66PiMfzitPM6mdm10SmTen0mFIt\nKreEEREXjLL9auDqYbbdCtyaR1xmlh9JVMolVrrhuyUV3UvKzFpMT3cyplRyw8BaiROGmdXV/HKJ\nbXv62bh1d9GhWJ05YZhZXVXHlHLDd+txwjCzuqqkY0p5EMLW44RhZnV1xJQJlA+byIpn3PDdapww\nzKzukp5SvsJoNU4YZlZ3lXKJVX3b2D/gnlKtxAnDzOqup1xi974B1m3ZWXQoVkdOGGZWd5XudEwp\n35ZqKU4YZlZ382clPaWcMFqLE4aZ1d3UiR0cPW0yKzxESEtxwjCzXPSUSx7mvMU4YZhZLirdJdY8\nt519+weKDsXqxAnDzHLRUy6xb3+w9rkdRYdideKEYWa5qHhMqZbjhGFmuThu5lTa5OlaW0luCUPS\n9ZL6JD02zPYLJT2aLvdKWlCzba2kZZIelrQ0rxjNLD+TOtuZN2OqrzBaSJ5XGDcAC0fY/iTwxoh4\nNfApYPGg7WdFxEkR0ZtTfGaWsx7PvtdScksYEXEXsGWE7fdGxPPpy/uAo/OKxcyKUSmXeGrzDnbv\n2190KFYHY6UN473AbTWvA7hd0oOSFo20o6RFkpZKWrpp06ZcgzSzA9PTXWIgYHWfrzJaQeEJQ9JZ\nJAnjYzXFp0fEycA5wPskvWG4/SNicUT0RkTvzJkzc47WzA5EdTIlDxHSGgpNGJJeDXwZOC8iNlfL\nI2JD+m8fcAtwSjERmtmhOObIqUxob3PDd4soLGFImgvcDFwUEStryqdKKlXXgbOBIXtamdnY1tne\nxnEzp7prbYvoyOvAkm4EzgRmSFoPXAl0AkTEtcAVwJHANZIA+tMeUWXglrSsA/hmRPwwrzjNLF89\n3SWWrn1+9Io25uWWMCLiglG2XwpcOkT5GmDBr+5hZs2oUi7x3Yc3sG33PkqTOosOxw5B4Y3eZtba\nqkOErHJPqabnhGFmuepJE4bbMZqfE4aZ5eroaZOZ3NnuJ75bgBOGmeWqrU1Uyl1+FqMFOGGYWe4q\n5ZKfxWgBThhmlrtKucSmbXvYsmNv0aHYIXDCMLPcVbrThm9fZTQ1Jwwzy121p9QqJ4ym5oRhZrkr\nHzaRwyZ1uB2jyTlhmFnuJFEpl1j5jLvWNjMnDDNriEp30lMqIooOxQ6SE4aZNURPucTWXfvo27an\n6FDsIDlhmFlDVMeUck+p5uWEYWYNUZ19b4XHlGpaThhm1hBHdk1kRtcEX2E0MScMM2uYZIgQ95Rq\nVrkmDEnXS+qTNOQUq0p8QdJqSY9KOrlm28WSVqXLxXnGaWaNUSmXWPXsNgYG3FOqGeV9hXEDsHCE\n7ecA89NlEfCPAJKmk0zpeipwCnClpGm5RmpmuevpLrFz736efmFX0aHYQchtilaAiLhL0rwRqpwH\nfC2Sjtn3STpC0stI5gK/IyK2AEi6gyTx3JhnvGaWr2pPqQu+dB+TO9sLjqZ1TJsyge9cdlru75Nr\nwshgNrCu5vX6tGy48l8haRHJ1Qlz587NJ0ozq4tXzT6cC0+dy/M7PWptPR3WoLnSi04YGqIsRij/\n1cKIxcBigN7eXt8YNRvDJnS08Zl3vKroMOwgFd1Laj0wp+b10cCGEcrNzKwgRSeMJcAfpL2lXgds\njYiNwI+AsyVNSxu7z07LzMysILnekpJ0I0kD9gxJ60l6PnUCRMS1wK3AucBqYCfwnnTbFkmfAh5I\nD3VVtQHczMyKkXcvqQtG2R7A+4bZdj1wfR5xmZnZgSv6lpSZmTUJJwwzM8vECcPMzDJxwjAzs0zU\nStMlStoEPFV0HCOYATxXdBAZNEuc0DyxOs76a5ZYx3qcx0TEzCwVWyphjHWSlkZEb9FxjKZZ4oTm\nidVx1l+zxNoscWbhW1JmZpaJE4aZmWXihNFYi4sOIKNmiROaJ1bHWX/NEmuzxDkqt2GYmVkmvsIw\nM7NMnDDMzCwTJ4w6kbRQ0gpJqyV9fIjt/0vSE5IelfTvko6p2bZf0sPpsqTgOC+RtKkmnktrtl0s\naVW6XFxwnJ+viXGlpBdqtjXyfF4vqU/SY8Nsl6QvpJ/jUUkn12xr5PkcLc4L0/gelXSvpAU129ZK\nWpaez6V5xpkx1jMlba35GV/MwgagAAAFqElEQVRRs23E702D4/xoTYyPpd/L6em2hp7TuokIL4e4\nAO3AL4DjgAnAI8AJg+qcBUxJ1y8Hvl2zbfsYivMS4Ooh9p0OrEn/nZauTysqzkH1/xi4vtHnM32v\nNwAnA48Ns/1c4DaSWSRfB9zf6POZMc7XV98fOKcaZ/p6LTBjDJ3TM4HvH+r3Ju84B9V9G/Djos5p\nvRZfYdTHKcDqiFgTEXuBbwHn1VaIiJ9ExM705X0kswg22qhxjuC3gDsiYktEPA/cASwcI3FeANyY\nUywjioi7gJHmajkP+Fok7gOOkPQyGns+R40zIu5N44Divp/VWEY7p8M5lO/3ATvAOAv7jtaTE0Z9\nzAbW1bxen5YN570kf3VWTZK0VNJ9kn4njwBTWeN8Z3pr4l8lVafKPdDPeCgyv1d6a+9Y4Mc1xY06\nn1kM91kaeT4P1ODvZwC3S3pQ0qKCYhrsNEmPSLpN0olp2Zg8p5KmkPwxcFNN8Vg8p6PKdQKlcURD\nlA3ZX1nSu4Fe4I01xXMjYoOk44AfS1oWEb8oKM7vATdGxB5JlwFfBd6Ucd96OZD3Oh/414jYX1PW\nqPOZxXCfpZHnMzNJZ5EkjDNqik9Pz+cs4A5JP0//ui7KQyTjH22XdC7wb8B8xug5JbkddU/891lD\nx9o5zcRXGPWxHphT8/poYMPgSpLeDPwp8PaI2FMtj4gN6b9rgDuB1xQVZ0RsrontS8Brs+7byDhr\nnM+gS/0Gns8shvssjTyfmUh6NfBl4LyI2FwtrzmffcAtJLd+ChMRL0bE9nT9VqBT0gzG4DlNjfQd\nHRPnNLOiG1FaYSG5UltDcmuk2th24qA6ryFpkJs/qHwaMDFdnwGsIqeGuoxxvqxm/R3Afen6dODJ\nNN5p6fr0ouJM6/WQNB6qiPNZ857zGL6B9rf5743eP2v0+cwY51xgNfD6QeVTgVLN+r3AwjzjzBBr\nd/VnTvKL9pfp+c30vWlUnOn2w0naOaYWfU7rsfiWVB1ERL+k9wM/IumpcX1EPC7pKmBpRCwBPgt0\nAf8iCeCXEfF24BXAP0kaILni++uIeKLAOD8g6e1AP8kX/ZJ03y2SPgU8kB7uqvjvl9iNjhOShsRv\nRfo/L9Ww8wkg6UaSXjszJK0HrgQ6089xLXArSU+p1cBO4D3ptoadz4xxXgEcCVyTfj/7IxlhtQzc\nkpZ1AN+MiB/mFWfGWN8FXC6pH9gFnJ9+B4b83hQYJyR/dN0eETtqdm34Oa0XDw1iZmaZuA3DzMwy\nccIwM7NMnDDMzCwTJwwzM8vECcPMzDJxwjAzs0ycMMzMLBM/uGc2Ckl/BlxIMrDdc8CDwFZgEckT\nxauBiyJip6QbSB4m+zXgGJIH9S4GTiMZMvyS9JjbgS8CbwaeBz4J/A3JE9cfioglkuYBXyd5Ghjg\n/RFxb76f1mx4vsIwG4GkXuCdJEO7/A+SgSMBbo6IX4+IBcBykgH7qqaRDNj4JySDOX4eOBF4laST\n0jpTgTsj4rXANuDTwFtIngy+Kq3TB7wlIk4Gfh/4Qi4f0iwjX2GYjewM4LsRsQtA0vfS8ldK+jRw\nBMmQLz+q2ed7ERGSlgHPRsSydN/HScYeehjYC1SHg1gG7ImIfek+89LyTuDqNMnsByr5fESzbJww\nzEY21JDZADcAvxMRj0i6hGRMoarqaL8DNevV19X/c/tqxsB6qV5EDEiq1vkT4FlgAcndgN0H/SnM\n6sC3pMxGdjfwNkmTJHWRjD4LUAI2Suokad/Iw+HAxogYAC4iGVDPrDC+wjAbQUQ8IGkJyVDZTwFL\nSRq8/wy4Py1bRpJA6u0a4CZJvwv8BNgxSn2zXHm0WrNRSOqKZHa3KcBdwKKIeKjouMwazVcYZqNb\nLOkEYBLwVScLG698hWFmZpm40dvMzDJxwjAzs0ycMMzMLBMnDDMzy8QJw8zMMvn/HQPpVOAhe2EA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x102229198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the results\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(gamma_values, error_values)\n",
    "plt.title(\"LASSO error\")\n",
    "plt.xlabel(\"gamma\")\n",
    "plt.ylabel(\"error\")\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(gamma_values, num_nonzeros)\n",
    "plt.title(\"Number of non-zeros\")\n",
    "plt.xlabel(\"gamma\")\n",
    "plt.ylabel(\"non-zeros\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}