05-ADMM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Distributed Optimization with the Alternating Direction Method of Multipliers\n",
    "\n",
    "In this lecture we will see the how the ADMM algorithm can be used for distributed optimization. Using the notatation from the course slides, we aim to solve the following optimization problem\n",
    "\n",
    "![](ADMM_problem.png)\n",
    "\n",
    "\n",
    "The algorithm is described by the following iteration:\n",
    "\n",
    "![](ADMM_algorithm.png)\n",
    "\n",
    "where the first N updates can be done in parallel.\n",
    "\n",
    "\n",
    "## Structure of this session:\n",
    "\n",
    "   * ADMM with 1 split (sequential algorithm).\n",
    "   * ADMM with 2 splits or more.\n",
    "   * Parallelize the splitting algorithm in Python using multiprocessing.\n",
    "   * We will see one application of the lasso, that of tomography reconstruction. \n",
    "   \n",
    "\n",
    "![](tomography.png)\n",
    "\n",
    "\n",
    "## References\n",
    "\n",
    "* [DSOML course notes](http://marcocuturi.net/Teaching/ENSAE/2017/DSOML/DSOML5.pdf)\n",
    "\n",
    "* [Distributed Optimization and Statistical\n",
    "Learning via the Alternating Direction\n",
    "Method of Multipliers](http://stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf), S. Boyd et al. 2010"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "# .. needed imports ..\n",
    "import numpy as np\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# Practical session: Lasso problem\n",
    "\n",
    "We will now use ADMM to solve a nonsmooth optimization problem.\n",
    "\n",
    "\n",
    "The Lasso is staistical model that has found many applications in machine learning, signal processing, biomedical imaging and genetics, to name a few. Given a regression dataset $A \\in \\mathbb{R}^{n \\times p}, b \\in \\mathbb{R}^n$, the Lasso is determined as the solution to the optimization problem\n",
    "\n",
    "$$\\text{argmin}_{\\theta} \\frac{1}{2n}\\|A \\theta - b\\|^2 + \\lambda \\|\\theta\\|_1 \\quad,$$\n",
    "where $\\lambda$ is a regularization parameter. The Lasso model is so useful because, given enough regularization, it gives solutions which are *sparse*, i.e., with many zero coefficients.\n",
    "\n",
    "\n",
    "\n",
    "# Part 1: a simple implementation\n",
    "\n",
    "We can reformulate the Lasso problem in the following equivalent form\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\text{argmin}_{\\theta} \\frac{1}{2n}\\|A \\theta - b\\|^2 + \\lambda \\|\\rho\\|_1 \\quad,\\\\\n",
    "&\\text{ subject to }~ \\theta=\\rho\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Now that we have written it in ADMM-friendly fashion, we can use the the above algorithm with $N=1$. This gives the following iteration:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\theta^{t+1} = \\text{argmin}_{\\theta} \\frac{1}{2n}\\|A \\theta - b\\|^2 + \\frac{\\tau}{2}\\|\\theta - \\rho^t + u^t\\|^2 \\\\\n",
    "&\\rho^{t+1} = \\text{argmin}_{\\theta} \\lambda \\|\\theta\\|_1 + \\frac{\\tau}{2} \\|\\theta - \\theta^{t+1} - u^t\\|^2\\\\\n",
    "& u^{t+1} = u^t + \\theta^{t+1} - \\rho^{t+1}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "It turns that for this particular problem many of the steps have a closed form solution. Using this we have the following closed form updates:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\theta^{t+1} = \\frac{1}{n}(A^T A / n + \\tau I)^{-1}(A^T b / n + \\tau (\\rho^t - u^t)) \\\\\n",
    "&\\rho^{t+1} = (\\theta^{t+1} + u^t - \\frac{\\tau}{\\lambda})_+ - (\\frac{\\tau}{\\lambda} - \\theta^{t+1} - u^t)_+\\\\\n",
    "& u^{t+1} = u^t + \\theta^{t+1} - \\rho^{t+1}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $(\\cdot)_+$ denote the positive part, i.e., $(x)_+ := \\text{max}(x, 0)$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "# .. create some fake data ..\n",
    "np.random.seed(0)\n",
    "\n",
    "n_samples, n_features = 100, 100\n",
    "A = np.random.randn(n_samples, n_features)\n",
    "w = np.random.randn(n_features)\n",
    "b = A.dot(w) + np.random.randn(n_samples)\n",
    "\n",
    "def obj_fun(theta, A, b, lamb):\n",
    "    # .. the lasso objective function ..\n",
    "    loss = (0.5 / n_samples) * np.linalg.norm(A.dot(theta) - b)**2\n",
    "    return loss + lamb * np.sum(np.abs(theta))\n",
    "\n",
    "def lasso_ADMM(A, b, max_iter=100, lam=1.):\n",
    "    # .. initialize variables ..\n",
    "    tau = 1.\n",
    "    n_samples, n_features = A.shape\n",
    "    rho = np.zeros(n_features)\n",
    "    u = np.zeros(n_features)\n",
    "\n",
    "    # .. to keep track of progress ..\n",
    "    obj_fun_history = []\n",
    "\n",
    "    # .. cache inverse matrix ..\n",
    "    AtA_inv = np.linalg.pinv(A.T.dot(A) / n_samples + tau * np.eye(n_features))\n",
    "    for i in range(max_iter):\n",
    "        theta = AtA_inv.dot(A.T.dot(b) / n_samples + tau * (rho - u))\n",
    "        rho = np.fmax(theta + u - lam /tau, 0) - np.fmax(-lam/tau - theta - u, 0)\n",
    "        u = u + theta - rho\n",
    "        obj_fun_history.append(obj_fun(theta, A, b, lam))\n",
    "\n",
    "    return theta, obj_fun_history\n",
    "    \n",
    "theta, func_vals = lasso_ADMM(A, b, lam=1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUXXV99/H355w5c8uQDElgjOZGhILWcjHhIqKdgPpQ\n6yP1jpU+oHVl2QulrT60LnvDrlVbpVZdai1YC7UtKXhD44OCmAEtGMhoCLcot3AnAWJCJiG3me/z\nx94TDpNzZs5MZs/J7P15rXVWztlnX77fEOY7v/377d9PEYGZmRVXqdkBmJlZc7kQmJkVnAuBmVnB\nuRCYmRWcC4GZWcG5EJiZFZwLgZlZwbkQmJkVnAuBmVnBtTQ7gEbMnTs3Fi9ePKFjd+zYwYwZMyY3\noGmgiHkXMWcoZt5FzBnGn3d/f/8zEXHEWPtNi0KwePFi1q5dO6Fj+/r66O3tndyApoEi5l3EnKGY\neRcxZxh/3pIebmQ/3xoyMys4FwIzs4JzITAzKzgXAjOzgnMhMDMrOBcCM7OCcyEwMyu4XBeCG+/d\nxKoH9zQ7DDOzQ1quC8HNv3ia6x7a2+wwzMwOabkuBG2VMnsHmx2FmdmhLdeFoL2lxN4hiIhmh2Jm\ndsjKdSFoq5QJYM/gULNDMTM7ZOW6ELRXygDs2utCYGZWT64LQVtLkt5udxSYmdWV60Iw3CLYvc8t\nAjOzenJeCJL0drlFYGZWV64LQVuL+wjMzMaS60Iw3CLYvc8tAjOzenJeCNwiMDMbS74Lwf5bQ24R\nmJnVk+tC0DbcWexbQ2ZmdWVeCCSVJf1M0qr082xJN0i6L/3z8KyuPdwi2O1bQ2ZmdU1Fi+Ai4N6q\nz38O3BgRxwA3pp8z0e4WgZnZmDItBJLmA78JfLlq8znAlen7K4Hfyur6Hj5qZja2rFsEnwEuBqp/\nEvdExJPp+6eAnqwu3uYHyszMxtSS1YklvQXYHBH9knpr7RMRIanmHNGSVgArAHp6eujr6xt3DBGB\nCH5x/4P06bFxHz+dDQwMTOjvbDorYs5QzLyLmDNkl3dmhQB4LfBWSW8G2oGZkv4D2CRpXkQ8KWke\nsLnWwRFxGXAZwLJly6K3t3dCQVRu+C4vedkCentfMaHjp6u+vj4m+nc2XRUxZyhm3kXMGbLLO7Nb\nQxHx0YiYHxGLgXOBH0bEecC3gfPT3c4Hrs0qBoBK2beGzMxG04znCP4eeKOk+4A3pJ8zUynJs4+a\nmY0iy1tD+0VEH9CXvn8WOGsqrgvQ6haBmdmocv1kMUCl5OGjZmajyX0haC3JD5SZmY0i94WgUvYU\nE2Zmo8l9IXCLwMxsdLkvBMnwUbcIzMzqyX8hKMFujxoyM6sr94WgteznCMzMRpP7QpAMH3WLwMys\nntwXglYXAjOzUeW+EFTKYpdvDZmZ1ZX7QtBahsGhYN+gi4GZWS25LwSVkgDcKjAzq6MAhSD50/0E\nZma15b4QtCbLFnsIqZlZHfkvBMO3htwiMDOrKfeFoJK2CFwIzMxqy38h2N9H4FtDZma15L4QtJaT\nW0O7PQOpmVlNuS8Ewy0Cr0lgZlZb7gvBcIvAfQRmZrXlvhDs7yPwrSEzs5pyXwj2P0fgW0NmZjXl\nvhBU/ByBmdmocl8IWvffGnKLwMysltwXgopvDZmZjSr3haAkpWsS+NaQmVktuS8EAO0tZfcRmJnV\nUYhC0FYpe4oJM7M6MisEktol3SbpDkl3S7ok3X6ipJ9IWidpraRTsophWHul5CkmzMzqaMnw3LuB\nMyNiQFIF+LGk64CPA5dExHWS3gx8EujNMA7aWkruLDYzqyOzQhARAQykHyvpK9LXzHT7LOCJrGIY\n1l5xH4GZWT1ZtgiQVAb6gaOBL0TEGkl/DHxf0qUkt6ZOzzIGSAuBbw2ZmdWk5Bf3jC8idQPfBC4E\nVgA3RcTXJb0bWBERb6hxzIp0X3p6epauXLlyQtceGBjgC/eU2TcEHzutY8I5TDcDAwN0dXU1O4wp\nVcScoZh5FzFnGH/ey5cv74+IZWPtNyWFAEDSXwE7gb8EuiMiJAnYFhEzRzt22bJlsXbt2gldt6+v\nj3/fOIPN23ex6sLXTegc01FfXx+9vb3NDmNKFTFnKGbeRcwZxp+3pIYKQZajho5IWwJI6gDeCGwg\n6RP49XS3M4H7sophWHul5OGjZmZ1ZNlHMA+4Mu0nKAFXR8QqSVuBz0pqAXaR3v7Jkh8oMzOrL8tR\nQ+uBk2ps/zGwNKvr1tJWKbHbk86ZmdVUjCeL3SIwM6urEIWgvVL2A2VmZnUUohC0tZTYMzjE0NDU\njJAyM5tOxlUIJB0u6fisgslKe7oogfsJzMwONGYhkNQnaaak2cBPgcslfTr70CZPe7qCvfsJzMwO\n1EiLYFZEPAe8Hfj3iDgVOOBJ4EPZcIvA00yYmR2okULQImke8G5gVcbxZKKtJUnTHcZmZgdqpBB8\nHPg+8EBE3C5pCVPwNPBkcovAzKy+MR8oi4hrgGuqPj8IvCPLoCbbC30EbhGYmY3USGfxr0i6UdJd\n6efjJf1F9qFNnraWtEXgzmIzswM0cmvocuCjwF7YP3XEuVkGNdmGWwQePmpmdqBGCkFnRNw2Ytu+\nLILJilsEZmb1NVIInpH0cpIlJpH0TuDJTKOaZPs7i10IzMwO0Mjso38AXAYcJ+lx4CHgvEyjmmT7\nh4/61pCZ2QEaGTX0IPAGSTOAUkRszz6sybV/igm3CMzMDjBmIUiXmKz+DEBEfDyjmCadh4+amdXX\nyK2hHVXv24G3APdmE0423FlsZlZfI7eG/rH6s6RLSZ40njYqZVGS+wjMzGqZyHoEncD8yQ4kS5Jo\nr3iVMjOzWhrpI7iTdOgoUAaOIJl/aFppr5Q915CZWQ2N9BG8per9PmBTREyrB8ogGULqzmIzswPV\nLQTpQjQAI4eLzpRERGzJLqzJ114pu4/AzKyG0VoE/SS3hFTjuwCWZBJRRpIWgW8NmZmNVLcQRMRR\nUxlI1txZbGZWWyN9BEg6HDiG5DkCACLi5qyCykJbS8krlJmZ1dDIqKEPAheRDBldB5wG3AqcmW1o\nk6u9Umbrzj3NDsPM7JDTyHMEFwEnAw9HxHLgJGBrplFloL3iUUNmZrU0Ugh2RcQuAEltEbEBODbb\nsCZfW4ufIzAzq6WRQvCYpG7gW8ANkq4FHh7rIEntkm6TdIekuyVdUvXdhZI2pNs/OfHwG9decR+B\nmVktjcw19Lb07d9IWg3MAr7XwLl3A2dGxICkCvBjSdcBHcA5wAkRsVvSkROMfVz8ZLGZWW2NdBZ/\nDlgZEbdExE2NnjgiAhhIP1bSVwC/B/x9ROxO99s87qgnwMNHzcxqa+TWUD/wF5IekHSppGWNnlxS\nWdI6YDNwQ0SsAX4FeJ2kNZJuknTyxEIfn6QQDDE0FGPvbGZWIEp+cW9gx2TKiXcA5wILI+KYhi+S\n9DF8E7gQWAmsBv6IZDTSfwNLYkQgklYAKwB6enqWrly5stHLvcjAwABdXV18f+Nertqwhy+c1cmM\nSq2HpfNlOO8iKWLOUMy8i5gzjD/v5cuX90fEmL+8N/RAWepo4DhgEeNcmCYitqb9C2cDjwHfSH/w\n3yZpCJgLPD3imMtI1kpm2bJl0dvbO55L7tfX10dvby/P9D/GVRvu4FWvPoVFc2ZM6FzTyXDeRVLE\nnKGYeRcxZ8gu7zFvDUn6pKT7SKaevhNYFhH/u4HjjkhbAkjqAN4IbCAZfbQ83f4rQCvwzIQzaFB3\nRwWArTv3Zn0pM7NppZEWwQPAayJivD+s5wFXSiqTFJyrI2KVpFbgK5LuAvYA54+8LZSF7s60EDzv\nQmBmVq2R4aP/MpETR8R6kqeQR27fA5w3kXMejOFCsM2FwMzsRSayVOW0NKujFYBtnm/IzOxFClQI\n3EdgZlZLo9NQl4Ge6v0j4pGsgspCa0uJztay+wjMzEZo5MniC4G/BjYBw5P1BHB8hnFloruj4haB\nmdkIjbQILgKOjYhnsw4ma7M6W9n2vPsIzMyqNdJH8CiwLetApkJ3R8WjhszMRmikRfAg0CfpuyQz\nigIQEZ/OLKqMdHdWuH/zwNg7mpkVSCOF4JH01Zq+pq3uzoo7i83MRmjkgbJLACR1pZ+n7a/Uszpa\n2bZzLxGBlP+J58zMGtHIXEOvkvQz4G7gbkn9kn41+9Am36yOCnsGh3je6xKYme3XSGfxZcCfRsSi\niFgEfBi4PNuwsrF/viEPITUz26+RQjAjIlYPf4iIPmBazuM8PAOpRw6Zmb2goVFDkv4S+Gr6+TyS\nkUTTziy3CMzMDtBIi+ADwBHAN9LXEem2aad7eOI5P1RmZrZfI6OGfkmyrOS05z4CM7MD1S0Ekj4T\nEX8s6Tskcwu9SES8NdPIMuDFaczMDjRai2C4T+DSqQhkKnRUylTKcovAzKxK3UIQEf3p2xMj4rPV\n30m6CLgpy8CyICl5qMx9BGZm+zXSWXx+jW0XTHIcU6a70xPPmZlVG62P4L3AbwNHSfp21VeHAVuy\nDiwrXpPAzOzFRusjuAV4EpgL/GPV9u3A+iyDylJ3Z4Untu5qdhhmZoeM0foIHgYelvQ+4ImI2AUg\nqQOYD2yckggn2ayOVu59cnuzwzAzO2Q00kdwNS8sUQkwCFyTTTjZ6+6ssHWnO4vNzIY1UghaImL/\nT870/bRdl2BWR4UdewbZOzg09s5mZgXQSCF4WtL+h8cknQM8k11I2Rp+qMwjh8zMEo1MOvch4D8l\nfR4QyRrG/yfTqDI0q+OFaSbmdrU1ORozs+ZrZK6hB4DT8rBCGUB3pyeeMzOrNmYhkNQGvANYDLQM\nL/EYER/PNLKMdHd44jkzs2qN9BFcC5wD7AN2VL1GJald0m2S7pB0t6RLRnz/YUkhae5EAp8oz0Bq\nZvZijfQRzI+Isydw7t3AmRExIKkC/FjSdRHxE0kLgDcBj0zgvAdleE0Cz0BqZpZopEVwi6RfG++J\nIzHcn1BJX8PTWf8TcDE1prfO2mHtLUgeNWRmNqyRQnAG0C/p55LWS7pTUkNTTEgqS1oHbAZuiIg1\n6fDTxyPijoOIe8JKJTGzvcI2P1RmZgaAIkb/pVzSolrb0ykoGruI1A18E7gIuBx4U0Rsk7QRWBYR\nBzyXIGkFsAKgp6dn6cqVKxu93IsMDAzQ1dX1om0X37yTJbNKfOiE9gmdczqolXfeFTFnKGbeRcwZ\nxp/38uXL+yNi2Vj7NdJHcNC3byJiq6TVJJ3ORwF3pKOP5gM/lXRKRDw14pjLgMsAli1bFr29vRO6\ndl9fHyOPnXfXj2nrbKW395QJnXM6qJV33hUxZyhm3kXMGbLLu5FC8F2SYiCgneQH+c+BXx3tIElH\nAHvTItABvBH4h4g4smqfjdRpEWRpVmerO4vNzFKNPFD2oo5iSa8Gfr+Bc88DrpRUJumLuDoiVk0o\nyknW3VHhkWfHHAFrZlYIjbQIXiQifirp1Ab2Ww+cNMY+i8d7/cngVcrMzF7QyJPFf1r1sQS8Gngi\ns4imwKyOpBAMDQWlkpodjplZUzUyfPSwqlcbSZ/BOVkGlbXuzlaGws8SmJnB6GsWfzUifgfYGhGf\nncKYMjf/8A4AHvvl8xw+Y9ourWBmNilGaxEslfRS4AOSDpc0u/o1VQFmYeHsTgAe3uIOYzOz0foI\nvgTcCCwB+kmGjw6LdPu0tCAtBI9s2dnkSMzMmq9uiyAiPhcRrwC+EhFLIuKoqte0LQIAXW0tzO1q\n5VEXAjOzsTuLI+L3piKQqbZgdicPP+tCYGbWyKihXFo0u9O3hszMKHAhWDi7kye2Ps+efUPNDsXM\nrKmKWwjmzGAo4Imtzzc7FDOzpipuIdg/hNS3h8ys2ApfCNxPYGZFV9hCcORhbbS1lDyE1MwKr7CF\noFRSOoTUTxebWbEVthDA8BBSdxabWbEVuhAsmN3JI8/uYKx1m83M8qzQhWDRnE527Blky449zQ7F\nzKxpCl0IPITUzKzghWDRnKQQeOSQmRVZoQvB/MPTZwk8+ZyZFVihC0F7pUzPzDbfGjKzQit0IQBY\nNHuGny42s0IrfCFIhpC6EJhZcRW+ECya08lTz+1i197BZodiZtYUhS8ERx/ZBcBtD21pciRmZs1R\n+EJw1iuOpGdmG1+66YFmh2Jm1hSFLwRtLWU+eMYSbnngWdY9urXZ4ZiZTbnCFwKA9566kFkdFb64\n+v5mh2JmNuUyKwSS2iXdJukOSXdLuiTd/ilJGyStl/RNSd1ZxdCorrYWzj99Mdffs4n7Nm1vdjhm\nZlMqyxbBbuDMiDgBOBE4W9JpwA3AqyLieOAXwEczjKFhF5y+mI5KmX92X4GZFUxmhSASA+nHSvqK\niLg+Ival238CzM8qhvGYPaOVc09ZwLXrnuBxL2hvZgWSaR+BpLKkdcBm4IaIWDNilw8A12UZw3h8\n8HVLiAi+euvDzQ7FzGzKaCoWZUn7Ab4JXBgRd6XbPgYsA94eNYKQtAJYAdDT07N05cqVE7r2wMAA\nXV1dDe//+Z/t4t4tg3y6t5O2siZ0zUPBePPOgyLmDMXMu4g5w/jzXr58eX9ELBtzx4iYkhfwV8BH\n0vcXALcCnY0cu3Tp0pio1atXj2v/Wx94Jhb92aq4as3DE77moWC8eedBEXOOKGbeRcw5Yvx5A2uj\ngZ+xWY4aOmJ4RJCkDuCNwAZJZwMXA2+NiENukp9Tj5rNcS85jCtu2eglLM2sELLsI5gHrJa0Hrid\npI9gFfB54DDgBknrJH0pwxjGTRLvf+1iNjy1nTWedsLMCqAlqxNHxHrgpBrbj87qmpPlnBNfxieu\n28AV/7OR05bMaXY4ZmaZ8pPFNbRXypx78kKuv+cpDyU1s9xzIajjvNMWAnDVmkeaHImZWbZcCOqY\nf3gnvcceyX+vfZS9g0PNDsfMLDMuBKN436kLeXr7bn5wz6Zmh2JmlhkXglH0HnskL53Vzn/d5ttD\nZpZfLgSjKJfEuacs5Ef3PcPGZ3Y0Oxwzs0y4EIzhPScvoFwSV93uVoGZ5ZMLwRh6ZrbzhlccydfW\nPsbufV7g3szyx4WgAe87dRHP7tjD9+56qtmhmJlNOheCBpxx9FyOmjuDf/ufjc0Oxcxs0rkQNKBU\nEue/ZhHrHt3Kzx75ZbPDMTObVC4EDXrnsgUc1tbCFbdsbHYoZmaTyoWgQV1tLbxr2QK+u/5JNj23\nq9nhmJlNGheCcbjg9MUMRvAfP/FSlmaWHy4E47BwTidnHdfDf615hF17PZTUzPLBhWCc3v/axTy7\nYw/X9D/W7FDMzCaFC8E4nf7yOZy2ZDaf+t4GNruvwMxywIVgnCTxibcfz659Q/zNd+5udjhmZgfN\nhWACjpo7g4vOOob/d+dTXH+3nzY2s+nNhWCCVrx+Cce95DD+8tq7eG7X3maHY2Y2YS4EE1Qpl/iH\ndxzP09t386Gv9rPteRcDM5ueXAgOwgkLurn0XSdw+8YtvPOfb+HRLTubHZKZ2bi5EBykt796Pld+\n4BQ2PbeLt33xf7jl/meaHZKZ2bi4EEyC018+l2/8/ul0trbw219ew/v/7TbuffK5ZodlZtYQF4JJ\ncvSRh3H9n7yeP/+N4+h/+Je8+XM/4oNXruXadY+zY/e+ZodnZlZXS7MDyJP2SpkP/frLee/JC/mX\nmx/g6z99jB/cu4n2SonXLJnDSQsP58QF3bzypTOZM6MVSc0O2czMhSALszorXHz2cXzkTcdy+8Yt\nfPfOJ7n1gWfp+8XTRCT7zGgts2B2Jy/r7mBuVxtzulqZPaOVme0VZna0cFh7hY7WMp2tZToqZdpa\nyrS2lGhtKVEpi0qpRKnkQmJmB8+FIEOlkjh1yRxOXTIHgOd27WX9o9v4xabtPLJlJ49u2ckT23Zx\n1xPbeHZgD/uGYnznF7SUSpRLoqUkSiVREpRLYu/evbTf8gNKEiJ5IloieZG+rzrXcOtk/7aqL+uV\nm6xbNOM9+46dO5nx05syieVQVsS8i5Tz37391zh58exMr+FCMIVmtlc445i5nHHM3AO+iwie27WP\n7bv28tzzyZ879w7y/J7ktWdwiN17B9m9b4h9Q8HewSH2DQaDEQymnyNgcCjZ9sTjT/CSeUcyOBQE\nEJFcI3jhzxeunf5ZFcv+7+olM76aNW4xgQts3vw8Rx7ZlUE0h7Yi5l2knDsq5cyvkVkhkNQO3Ay0\npdf5WkT8taTZwH8Di4GNwLsjovDrP0piVkeFWR0VOPzgz9fX9yy9vccf/Immkb6+Pnp7lzY7jClX\nxLyLmHOWshw1tBs4MyJOAE4EzpZ0GvDnwI0RcQxwY/rZzMyaJLNCEImB9GMlfQVwDnBluv1K4Ley\nisHMzMam6vvBk35yqQz0A0cDX4iIP5O0NSK60+8F/HL484hjVwArAHp6epauXLlyQjEMDAzQ1VWM\ne4nViph3EXOGYuZdxJxh/HkvX768PyKWjbljRGT+ArqB1cCrgK0jvvvlWMcvXbo0Jmr16tUTPnY6\nK2LeRcw5oph5FzHniPHnDayNBn5GT8mTxRGxNS0EZwObJM0DSP/cPBUxmJlZbZkVAklHSBq+BdQB\nvBHYAHwbOD/d7Xzg2qxiMDOzsWX5HME84Mq0n6AEXB0RqyTdClwt6XeBh4F3ZxiDmZmNIbNCEBHr\ngZNqbH8WOCur65qZ2fhkOmposkh6mqT1MBFzgSIuElDEvIuYMxQz7yLmDOPPe1FEHDHWTtOiEBwM\nSWujkeFTOVPEvIuYMxQz7yLmDNnl7fUIzMwKzoXAzKzgilAILmt2AE1SxLyLmDMUM+8i5gwZ5Z37\nPgIzMxtdEVoEZmY2ilwXAklnS/q5pPsl5XK6a0kLJK2WdI+kuyVdlG6fLekGSfelf07CKgeHFkll\nST+TtCr9XIScuyV9TdIGSfdKek3e85b0J+m/7bskXSWpPY85S/qKpM2S7qraVjdPSR9Nf7b9XNL/\nOphr57YQpE80fwH4DeCVwHslvbK5UWViH/DhiHglcBrwB2meRVj34SLg3qrPRcj5s8D3IuI44ASS\n/HObt6SXAX8ELIuIVwFl4FzymfMVJPOxVauZZ/r/+LnAr6bHfDH9mTchuS0EwCnA/RHxYETsAVaS\nrIWQKxHxZET8NH2/neQHw8vI+boPkuYDvwl8uWpz3nOeBbwe+FeAiNiTTuiY67xJZkDokNQCdAJP\nkMOcI+JmYMuIzfXyPAdYGRG7I+Ih4H6Sn3kTkudC8DLg0arPj6XbckvSYpJpPdYAPRHxZPrVU0BP\nk8LKymeAi4Ghqm15z/ko4Gng39JbYl+WNIMc5x0RjwOXAo8ATwLbIuJ6cpzzCPXynNSfb3kuBIUi\nqQv4OvDHEfFc9XfpvOS5GR4m6S3A5ojor7dP3nJOtQCvBv45Ik4CdjDilkje8k7viZ9DUgRfCsyQ\ndF71PnnLuZ4s88xzIXgcWFD1eX66LXckVUiKwH9GxDfSzXle9+G1wFslbSS55XempP8g3zlD8lvf\nYxGxJv38NZLCkOe83wA8FBFPR8Re4BvA6eQ752r18pzUn295LgS3A8dIOkpSK0nHyrebHNOkS5f7\n/Ffg3oj4dNVXuV33ISI+GhHzI2IxyX/XH0bEeeQ4Z4CIeAp4VNKx6aazgHvId96PAKdJ6kz/rZ9F\n0g+W55yr1cvz28C5ktokHQUcA9w24as0sozZdH0BbwZ+ATwAfKzZ8WSU4xkkzcX1wLr09WZgDsko\ng/uAHwCzmx1rRvn3AqvS97nPGTgRWJv+9/4WcHje8wYuIVnU6i7gq0BbHnMGriLpB9lL0vr73dHy\nBD6W/mz7OfAbB3NtP1lsZlZweb41ZGZmDXAhMDMrOBcCM7OCcyEwMys4FwIzs4JzIbBcktQnadLX\ndpX0qXQmzE+N2P7W4RluJf3WZE5wKOlESW+udS2zyeDho5ZLkvqAj0TE2kk+7zaSsdyDo+xzBcmz\nDV8bx3lbImJfne8uIJl98w/HGa5ZQ9wisKaRtDidU//y9Lfs6yV1pN/t/41e0tx0OgkkXSDpW+nc\n7Bsl/aGkP00nYfuJpNlVl/gdSevSeexPSY+fkc77flt6zAEz0irxqfS4OyW9J93+baAL6B/eVnXM\nBZI+L+l04K3Ap9Jrvzx9fU9Sv6QfSTouPeYKSV+StAb4pKRTJN2axnWLpGPTp+I/DrwnPd97hq9V\n9Xf4Q0nrJd0oaWHVuT+XnudBSe9Mt8+TdHPV38vrJuk/p01nzX6azq/ivoDFJOspnJh+vho4L33f\nR/JbMMBcYGP6/gKSKXcPA44AtgEfSr/7J5JJ94aPvzx9/3rgrvT931Vdo5vkyfMZI+J6B3ADydz3\nPSTTHMxLvxuok8sFwOfT91cA76z67kbgmPT9qSRTYgzvtwoop59nAi3p+zcAXx957hrX+g5wfvr+\nA8C3qs59Dckve68kmZId4MOkT9mn+R3W7H8HfjX/1TJqlTDL3kMRsS59309SHMayOpK1F7ant2q+\nk26/Ezi+ar+rIJnnXdJMSd3Am0gmrPtIuk87sJAXL3BzBnBVJLd/Nkm6CTiZCcxVlc4KezpwTTJV\nDpBMkTDsmnjhNtMs4EpJx5BMG1Jp4BKvAd6evv8q8Mmq774VEUPAPZKGpy++HfhKOlHht6r+7q3A\nXAis2XZXvR8EOtL3+3jh1mX7KMcMVX0e4sX/pkd2gAUg4B0R8fOJBjxOJWBrRJxY5/sdVe//lqTI\nvU3J2hJ9B3nt6r8nwf6i+HqSRX2ukPTpiPj3g7yOTXPuI7BD1UZgafr+nRM8x/C9/TNIFjTZBnwf\nuDCdyRJJJ9U47kck9+TLko4gubU0npkdt5PcuiKStSEekvSu9HqSdEKd42bxwlTCF9Q6Xw23kMzA\nCvC+NPa6JC0CNkXE5SSru7161EysEFwI7FB1KfB7kn5G0kcwEbvS479EMpMjJL91V4D1ku5OP4/0\nTZLZPe8AfghcHMkU0I1aCfzftNP35SQ/oH9X0h3A3dRfMvWTwCfSmKtbNquBVw53Fo845kLg/ZLW\nA79Dso7zaHqBO9JrvIdkDWQrOA8fNTMrOLcIzMwKzoXAzKzgXAjMzArOhcDMrOBcCMzMCs6FwMys\n4FwIzMyKjokUAAAADUlEQVQKzoXAzKzg/j86wIS8gG1rdQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11701c860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# lets plot the objective values of the function\n",
    "# to make sure it has converged\n",
    "plt.plot(func_vals)\n",
    "plt.ylabel('function values')\n",
    "plt.xlabel('numbe of iterations')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGn1JREFUeJzt3X2UXHWd5/H3h6aRgkh6WLAxDRqWZaIsGQlpFSczTDfq\nBB+Oxoy7wg6zsriTnTOr4K7ESTazM+KsJzmblZ2ZMz6AM466uvTZlRBZkAkPoYcRFUkMGp4yMjIK\nHQVc7GwCDXbCd/+o21DpdHXferh9q+79vM6p0123bt37/Vbdqu+9v/ur31VEYGZmdlTeAZiZWWdw\nQTAzM8AFwczMEi4IZmYGuCCYmVnCBcHMzAAXBDMzS7ggmJkZ4IJgZmaJo/MOoBEnnXRSLF68uKnn\nPvPMMxx//PHtDagLlDHvMuYM5cy7jDlD43nv3LnzZxFx8lzzdVVBWLx4MTt27GjquaOjowwNDbU3\noC5QxrzLmDOUM+8y5gyN5y3pR2nmc5ORmZkBLghmZpZwQTAzM8AFwczMEi4IZmYGuCCYmVnCBcHM\nus7WXWOs2LSd3WP7WLFpO1t3jeUdUiF01e8QzMy27hpj/ZbdTEwegtNgbHyC9Vt2A7Bq2UDO0XU3\nHyGYWVfZvG1PtRjUmJg8xOZte3KKqDhcEMysq+wdn2houqXngmBmXWVRX6Wh6ZaeC4KZdZW1K5dQ\n6e05bFqlt4e1K5fkFFFx+KSymXWVqRPH1XMG+xnoq7B25RKfUG4DFwQz6zqrlg2watkAo6OjfOi3\nh/IOpzDcZGRmZoALgpmZJVwQzMwMyLEgSDpN0p2SHpT0gKQr8orFzMzyPal8EPhIRHxX0suBnZJu\ni4gHc4zJzKy0cjtCiIifRMR3k//3Aw8B7jdmZpaTjjiHIGkxsAy4J99IzMzKSxGRbwDSAuBvgU9E\nxJYZHl8DrAHo7+9fPjIy0tR6Dhw4wIIFC1oJtSuVMe8y5gzlzLuMOUPjeQ8PD++MiMG55su1IEjq\nBW4CtkXE1XPNPzg4GDt27GhqXaOjowwNDTX13G5WxrzLmDOUM+8y5gyN5y0pVUHIs5eRgL8CHkpT\nDMzMLFt5nkNYAfwOcIGk+5Lb23OMx8ys1HLrdhoR3wCU1/rNzOxwHdHLyMzM8ueCYGZmgAuCmZkl\nXBDMzAxwQTAzs4QLgpmZAS4IZmaWcEEwMzPABcHMzBIuCGZmBrggmJlZwgXBzMwAFwQzM0u4IJiZ\nGeCCYGZmCRcEMzMDXBDMzCzhgmBmZoALgpmZJVwQzMwMcEEwM7OEC4KZmQEuCGZmlnBBMDMzwAXB\nzMwSLghmZga4IJiZWcIFwczMADg67wDMzGpt3TXG5m172Ds+waK+CmtXLmHVsoG8wyoFFwQz6xhb\nd42xfstuJiYPATA2PsH6LbsBXBTmgZuMzKxjbN6258ViMGVi8hCbt+3JKaJycUEws46xd3yioenW\nXi4IZtYxFvVVGppu7eWCYGYdY+3KJVR6ew6bVuntYe3KJTlFVC4+qWxmHWPqxLF7GeUj14Ig6fPA\nO4EnI+LsPGMxs86watmAC0BO8m4y+gJwYc4xmJkZOReEiLgLeDrPGMzMrCrvIwQzM+sQioh8A5AW\nAzfVO4cgaQ2wBqC/v3/5yMhIU+s5cOAACxYsaDLK7lXGvMuYM5Qz7zLmDI3nPTw8vDMiBuear+N7\nGUXEtcC1AIODgzE0NNTUckZHR2n2ud2sjHmXMWcoZ95lzBmyy9tNRmZmBuRcECRdB3wLWCLpcUkf\nyDMeM7Myy7XJKCIuznP9Zmb2EjcZmZkZ4IJgZmYJFwQzMwNcEMzMLNHxv0MwM2uGr83cOBcEMysc\nX5u5OW4yMrPC8bWZm+MjhCb5cNSsc2V5beYif/ZdEJrgw1Gzzraor8LYDF/+rV6bueiffTcZNaEb\nDke37hpjxabt7B7bx4pN29m6ayzvkMzmTVbXZu6Gz34rfITQhCwPR9vhsL2Y04q3F2M2l6yuzdzp\nn/1WuSA0IavD0XaZbS/GBaE8itzWnUYW12bu9M9+q9xk1ISsDkfbpeh7MTa3qaPEsfEJgpeOEt10\n2JpO/+y3KlVBUNUlkv4ouf8qSW/INrTOtWrZABtXL2Wgr4KAgb4KG1cv7Zi9r3p7K0XZi7G5Fb2t\nOy+d/tlvVdomo08DLwAXAB8H9gPXA6/PKK6Ol8XhaLusXbnksJ4QUKy9mE7Uac0zPkrMTprPfqdt\nD2mlLQhvjIhzJe0CiIifSzomw7isBbUn1GA/A120QXajTuyKWPS27k7WidtDWmnPIUxK6gECQNLJ\nVI8YrEOtWjbA3esuYOnAQu5ed0HHb4jdrBObZ4re1t3JOnF7SCvtEcKfAzcAr5D0CeC9wB9mFpVZ\nF+nE5pmsul3a3Dpxe0grVUGIiK9I2gm8GRCwKiIeyjQysy7Rqc0znXyeq8g6dXtIY9YmI0knTt2A\nJ4HrgP8JPJFMMys9N89YrW7eHuY6QthJ9byBgFcBP0/+7wN+DJyeaXRmXcDNM1arm7eHWQtCRJwO\nIOlzwA0R8fXk/tuAVdmHZ0XSrV3x0nDzTDnV26a7dXtIe1L5vIj43ak7EXGLpP+aUUxWQN3cFc+K\nrdkdlSJu02m7ne6V9IeSFie3DcDeLAOzYunmrnhWXK0M8VHEbTptQbgYOJlq19MbgFck08xS6eau\neFZcrXypF3GbTtvt9GngioxjsQLr5q54VlytfKkXcZtOO7jdnZK2T79lHZwVRzd3xbPiamUgyCJu\n02lPKl9Z8/+xwG8BB9sfjhVVN3fFs+JqZSDIIm7TaZuMdk6bdLek72QQjxVYt3bFs+Jq9Uu9aNt0\nqoIw7VfJRwHLgYWZRGRmNo+K9qXeirRNRrW/WD4IPAp8IKugzMxs/qUtCK+NiOdqJ0h6WQbxmJlZ\nTtIWhG8C506b9q0ZplmbFHmYBzPrTLMWBEmnAANARdIyqk1GACcAx2UcW2kV8SfxNj+8I2GtmOsI\nYSVwKXAqcHXN9P3Af8ooptKb7deT/nBbPeMTk6y/wzsS1rxZf5gWEV+MiGHg0ogYrrm9KyK2tLpy\nSRdK2iPpEUnrWl1eURTxJ/GWvSf2PVe4sXVsfs3VZHRJRHwZWCzpP05/PCKunuFpqSTXaP4U8Fbg\nceBeSTdGxIPNLnMmU4fQF522nw9fdSsSjD87ycJK75z/L+qrMPyak7nz4afYOz5xxGNTh+O1h+mz\nLSvt/FEnlwDOqckhTXx90Pb4Wp2flDE1+/9lZ0ywYdP2GeOrF/dsr2XabaWR16OV/+vFetkZLzDT\nPt7e8YnMY0r7uWnltazXBFYvtzTzp90u0zS/pdnOGt1mpr+WtZ/rLCii3tcPSPp3EXGNpD+e6fGI\nuKrpFUtvAj4WESuT++uTZW6s95zBwcHYsWNH6nXUtsV/ZOlBPrk77Tn0dCq9PfzW8gGu3zl2xJ5Z\nO+ZvR3z/eXnwJzuVSXzNzL9x9VKAI34d2k5T73W9+LJ6H+b7/Z2u3jbeV+nl+YMv5BJTs2Z77zau\nXvrSF/gtt7H+m4fq5nbE/NPOz9XOA7Nvl9OXVavecufaHprZZiq9PWz81R5Wve2tqeYHkLQzIgbn\nmm+uC+Rck/xt+ot/FgPAYzX3Hwfe2M4VzNQW304Tk4e47p7HODRLUW1l/lZNTB7i6WcOMTHZM/fM\nZJ9PbfPFfHw51Ysvq/dhvt/fNCq9PUjz83q302zvXe25tGozWf2W7+nzzzW66Wyv02zn8eotd67t\noZltZmLyEE/sm0w9fyNmPUJ4cSbpZOB3gcXUFJGIuKzpFUvvBS6MiH+b3P8d4I0R8cFp860B1gD0\n9/cvHxkZSb2O3WP7Xvy/vwJPlLAJvox5lzFnODLvY3qOon/hsTz29LP5BZWRpQPVgRKefHpfqvd6\nav7a74RW112rHcttRH8FXnFi+sEihoeHWz9CqPE14O+A24F27WqMAafV3D81mXaYiLgWuBaqTUZD\nQ0OpV7Bh0/YXh6fNoskIoEdqqLpnPf90Vy49xCd3pztCaGZ9jc4/kIwiOdOwwe1S+17Xi6/V17We\nrJabRm3eA30V7l53AQAraj4H3aTeaznQV+FDvz0EwDUj/4dP7p590Oba+TfUeS3Sbpe1y6pVb7lZ\nbQ/rz3mBf9nAd2FaaS+Qc1xE/EFE/K+IuH7q1uK67wXOlHS6pGOAi4AbW1zmYWYanradKr09XPzG\n01KvI+v5Z3r+icf3ZrY+AYciXvxxSprlr125JPP3pXZ9M+XT6uva6Prm2/TROufr9W6n2d672tz6\nFx47a25pXou02+Vso6DWW25W21n/wmPbuswpaXeZb5L09oj4ertWHBEHJX0Q2Ab0AJ+PiAfatXw4\nfCRD2E9fG3tL1PYQGHz1iQ31smnX/Kl6Ge37ARtXn9nW+H7+7CTipd5QU4NcBdR9jWfqpZFVrxc4\nyECdfOrl2c6eMWne3yx68UzPe6bPQbf1Mqr33k3pq/SycfVZqXsGpRndtJleRrMtd6bt4efP1j8H\n0DfLa1n7uc5C2nMI+4HjgeeBSZLPf0SckElUdTTay6jW6OgojTQ3FUUWeddrgqhtpsiT3+vy6Nac\nT19384zdywU8uukdcz6/0bzb0stoSkS8PPWarfD8wzmz1nTq5TfTXkLz3BluZ0hq/1la63itXHbQ\nzDr38ptpTyp/Gvg28Lnk9m3gfwN7JP1mRrFZh+rUjdmsW6xaNsDG1UsZ6Ksgqs2t9X70Np/S7uHv\nBT4wddJX0lnAx4GPAluAW7MJzzpRN19L1qOBWqfoxCu1pS0Iv1zbAygiHpT0moj4oZS206EVSSdu\nzHPxsOJms0tbEB6Q9Blg6mfC7wMeTK6als1vqM3azMOKl5ePDNNJWxAuBX4f+HBy/27gSqrFYLj9\nYZm1n3tHlVM7jgzLUlDSdjudAD6Z3KY70NaIzDLSqV39LFutHhmWqakxbbfTMyV9VdKDkn44dcs6\nOLN2cu+o7rZ11xgrNm3n9HU3s2LTdrbuOmLosxm1emQ41wipRZK22+lfA58BDlJtIvoS8OWsgjLL\nQqd29bO5Te2lj41PELy0lz4+MfcpzFZ/N1Ompsa05xAqEXGHJEXEj4CPSdoJ/FGGsZm1XTf2jrL6\ne+lprguwduWSGS9ek/bIsExNjWmPEJ6XdBTwA0kflPQeYEGGcZmZvaje3vgvDr0w53NbPTIsU1Nj\n2iOEK4DjgMuBPwEuAN6fVVBmZrXq7aUf05Nun7aVI8Nu/iFmo9L2Mro3+fcA8G+yC8fM7Ej1mn36\nFx4zL+svS1PjrAVB0qwXrImId7U3HDOzI9XbS8/qugBlNdcRwpuAx4DrgHsg9cWxzMzaaqa99NFR\nF4R2mqsgnAK8FbgY+FfAzcB17b6ymZmZ5W/WMzIRcSgi/iYi3g+cBzwCjCaXvjQzswKZ86RyMoDd\nO6geJSwG/hy4IduwzMxsvs11UvlLwNnA14GrIuL+eYnKzMzm3VxHCJcAz1D9HcLlNdc+EBARcUKG\nsZmZ2TyatSBERNpfMpuZWZfzF76ZmQEuCGZmlnBBMDMzwAXBzMwSLghmZga4IJiZWSLt9RDMusrW\nXWOlGL/erJ1cEKxwxicmWX/HS2PnT11/F3BRMJuFm4yscJ7Y99yM19/dvG1PThGZdQcXBCucetfZ\nrXddXjOrckGwwql3nd1FfZV5jsSsu7ggWOH0LzyWSm/PYdMqvT2sXbkkp4jMuoMLghVOX6WXjauX\nMtBXQcBAX4WNq5f6hLLZHHLpZSTpXwAfA14LvCEiduQRhxXXTNffNbPZ5dXt9H5gNXBNTuufd+4X\nb2adLpeCEBEPAdRccKfQtu4aY/0W94s3s87mcwjzYPO2Pe4Xb2YdTxGRzYKl24FTZnhoQ0R8LZln\nFLhytnMIktYAawD6+/uXj4yMNBXPgQMHWLBgQVPPbdXusX11H1s6sDDTdeeZd17KmDOUM+8y5gyN\n5z08PLwzIgbnmi+zgpBGmoJQa3BwMHbsaO788+joKENDQ009t1UrNm1nbIYfRQ30Vbh73QWZrjvP\nvPNSxpyhnHmXMWdoPG9JqQqCm4zmwdqVS9wv3sw6Xi4FQdJ7JD0OvAm4WdK2POKYL6uWDbhfvJl1\nvLx6Gd0A3JDHuvPifvFm1uncZGRmZoALgpmZJVwQzMwMcEEwM7OEC4KZmQEuCGZmlnBBMDMzwAXB\nzMwSLghmZga4IJiZWcIFwczMABcEMzNLuCCYmRnggmBmZgkXBDMzA1wQzMws4YJgZmaAC4KZmSVc\nEMzMDHBBMDOzhAuCmZkBLghmZpZwQTAzM8AFwczMEi4IZmYGuCCYmVnCBcHMzAAXBDMzS7ggmJkZ\n4IJgZmYJFwQzMwNcEMzMLOGCYGZmgAuCmZklXBDMzAzIqSBI2izpYUnfl3SDpL484jAzs5fkdYRw\nG3B2RPwK8PfA+pziMDOzRC4FISJujYiDyd1vA6fmEYeZmb2kE84hXAbckncQZmZlp4jIZsHS7cAp\nMzy0ISK+lsyzARgEVkedQCStAdYA9Pf3Lx8ZGWkqngMHDrBgwYKmntvNyph3GXOGcuZdxpyh8byH\nh4d3RsTgnDNGRC434FLgW8BxaZ+zfPnyaNadd97Z9HO7WRnzLmPOEeXMu4w5RzSeN7AjUnzHHt1c\nfWqNpAuBjwK/ERHP5hGDmZkdLq9zCH8BvBy4TdJ9kj6bUxxmZpbI5QghIv5ZHus1M7P6cikI1nm2\n7hpj87Y97B2fYFFfhbUrl7Bq2UDeYZnZPHJBMLbuGmP9lt1MTB4CYGx8gvVbdgO4KJiVSCf8DsFy\ntnnbnheLwZSJyUNs3rYnp4jMLA8uCMbe8YmGpptZMbkgGIv6Kg1NN7NickEw1q5cQqW357Bpld4e\n1q5cklNEZpYHn1S2F08cu5eRWbm5IBhQLQouAGbl5iYjMzMDXBDMzCzhgmBmZoALgpmZJVwQzMwM\ncEEwM7NEZpfQzIKkp4AfNfn0k4CftTGcblHGvMuYM5Qz7zLmDI3n/eqIOHmumbqqILRC0o5Ic03R\ngilj3mXMGcqZdxlzhuzydpORmZkBLghmZpYoU0G4Nu8AclLGvMuYM5Qz7zLmDBnlXZpzCGZmNrsy\nHSGYmdksSlEQJF0oaY+kRyStyzueLEg6TdKdkh6U9ICkK5LpJ0q6TdIPkr+/lHes7SapR9IuSTcl\n98uQc5+kr0p6WNJDkt5U9Lwl/Ydk275f0nWSji1izpI+L+lJSffXTKubp6T1yXfbHkkrW1l34QuC\npB7gU8DbgLOAiyWdlW9UmTgIfCQizgLOA/59kuc64I6IOBO4I7lfNFcAD9XcL0POfwb8TUS8Bngd\n1fwLm7ekAeByYDAizgZ6gIsoZs5fAC6cNm3GPJPP+EXAP0+e8+nkO68phS8IwBuARyLihxHxC2AE\neHfOMbVdRPwkIr6b/L+f6hfEANVcv5jM9kVgVT4RZkPSqcA7gL+smVz0nBcC5wN/BRARv4iIcQqe\nN9Xrt1QkHQ0cB+ylgDlHxF3A09Mm18vz3cBIRDwfEY8Cj1D9zmtKGQrCAPBYzf3Hk2mFJWkxsAy4\nB+iPiJ8kD/0U6M8prKz8KfBR4IWaaUXP+XTgKeCvk6ayv5R0PAXOOyLGgP8G/Bj4CbAvIm6lwDlP\nUy/Ptn6/laEglIqkBcD1wIcj4v/VPhbVLmWF6VYm6Z3AkxGxs948Rcs5cTRwLvCZiFgGPMO0ppKi\n5Z20mb+bajFcBBwv6ZLaeYqWcz1Z5lmGgjAGnFZz/9RkWuFI6qVaDL4SEVuSyU9IemXy+CuBJ/OK\nLwMrgHdJ+keqTYEXSPoyxc4ZqnuBj0fEPcn9r1ItEEXO+y3AoxHxVERMAluAX6XYOdeql2dbv9/K\nUBDuBc6UdLqkY6iegLkx55jaTpKotik/FBFX1zx0I/D+5P/3A1+b79iyEhHrI+LUiFhM9X3dHhGX\nUOCcASLip8BjkpYkk94MPEix8/4xcJ6k45Jt/c1Uz5MVOeda9fK8EbhI0ssknQ6cCXyn6bVEROFv\nwNuBvwf+AdiQdzwZ5fhrVA8jvw/cl9zeDvwTqr0SfgDcDpyYd6wZ5T8E3JT8X/icgXOAHcn7vRX4\npaLnDVwFPAzcD/wP4GVFzBm4jup5kkmqR4MfmC1PYEPy3bYHeFsr6/Yvlc3MDChHk5GZmaXggmBm\nZoALgpmZJVwQzMwMcEEwM7OEC4IVkqRTJI1I+gdJOyV9XdIvN7Gcy5PRRL+S9PW+XdJ9kt6XDBlR\nd6BESe9qdnTdZDTT32/muWbNcrdTK5zkh0vfBL4YEZ9Npr0OOCEi/q7BZT0MvCUiHpd0HvBfIuIt\nbQ/6yPUupvq7irOzXpfZFB8hWBENA5NTxQAgIr4HfEPS5mQ8/d2S3jf1uKS1ku6V9H1JVyXTPgv8\nU+AWSX8AfBl4fXKEcIakUUmDybwXSvqupO9JuiOZdqmkv0j+P1nS9ck67pW0Ipn+sWT8+1FJP5R0\neRLSJuCMZF2bJb1S0l3J/fsl/XrWL6KVz9F5B2CWgbOBmQa8W031F76vA04C7pV0F7CU6k/+3wAI\nuFHS+RHxe5IuBIYj4meS7gGujIh3AlQPRKpf9sDngPMj4lFJJ86w7j8D/ntEfEPSq4BtwGuTx15D\ntYi9HNgj6TNUB6s7OyLOSdbxEWBbRHwiGe/+uFZeILOZuCBYmfwacF1EHKI6WNjfAq+nem2B3wR2\nJfMtoFog7kq53POAu6I6Hj0RMX0se6gOznbWVBEBTkhGpgW4OSKeB56X9CQzD+F8L/D5ZADDrRFx\nX8rYzFJzQbAiegB4bwPzC9gYEddkFA9Um2fPi4jnDltxtUA8XzPpEDN8LiPiLknnU70Y0BckXR0R\nX8owXishn0OwItoOvEzSmqkJkn4FGAfep+o1mE+memTwHarNN5dN7bFLGpD0igbW923g/GS0Seo0\nGd0KfKgmnnPmWOZ+qk1IU/O/GngiIj5H9epw5zYQn1kqPkKwwomIkPQe4E+Tk8HPAf8IfJhqc9D3\nqI4M+9GoDiX9U0mvBb6V7LEfAC4h5dj6EfFUUny2SDoqed5bp812OfApSd+n+rm7C/i9WZb5fyXd\nreqF1m+hOsLnWkmTSXz/Ok1sZo1wt1MzMwPcZGRmZgkXBDMzA1wQzMws4YJgZmaAC4KZmSVcEMzM\nDHBBMDOzhAuCmZkB8P8B2/ZlhwcEWDsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11abdc7f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# let's also look at the coefficients of the solution\n",
    "plt.scatter(np.arange(len(theta)), theta)\n",
    "plt.ylabel('Magnitude')\n",
    "plt.xlabel('Coefficients')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They are **sparse**, i.e., with many coefficients set exactly to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above example, both steps are closed form but this is not necessarily the case. For example, in the case of L1-regularized logistic regression,\n",
    "$$\n",
    "    \\text{argmin}_{\\theta} \\frac{1}{n} \\sum_{i=1}^n -\\log(1 + e^{-b_i a_i^T x})\n",
    "$$\n",
    "\n",
    "\n",
    "<span style=\"color: red; font-weight: bold\">Exercise (10 min): Logistic Regression. Adapt the above code to optimize an L1 regularized logistic regression (i.e. when the loss is the logistic loss instead of the squared loss).</span>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sol.x is a solution to the logistic regression problem\n"
     ]
    }
   ],
   "source": [
    "from scipy import optimize\n",
    "\n",
    "# .. make b be {-1, 1} since its a classification problem ..\n",
    "b = np.sign(A.dot(w) + np.random.randn(n_samples))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# this solves the logistic regression problem\n",
    "# from scipy import optimize\n",
    "print('sol.x is a solution to the logistic regression problem')\n",
    "\n",
    "\n",
    "def ADMM_logloss(A, b, max_iter=100, lam=1.):\n",
    "    # .. initialize variables ..\n",
    "    tau = 1.\n",
    "    n_samples, n_features = A.shape\n",
    "    rho = np.zeros(n_features)\n",
    "    u = np.zeros(n_features)\n",
    "    theta = np.zeros(n_features)\n",
    "\n",
    "    # .. to keep track of progress ..\n",
    "    obj_fun_history = []\n",
    "\n",
    "    # .. cache inverse matrix ..\n",
    "    for i in range(max_iter):\n",
    "        def logistic_loss(theta):\n",
    "            p = b * A.dot(theta)\n",
    "            t = 1.0 + np.exp(-p)\n",
    "            return -np.sum(np.log(1.0 / t)) / n_samples\\\n",
    "            + (tau / 2.) * np.linalg.norm(theta - rho + u) ** 2\n",
    "\n",
    "        def logistic_grad(theta):\n",
    "            t = b * A.dot(theta)\n",
    "            z = 1. / (1 + np.exp(-t))\n",
    "            z0 = (z - 1) * b\n",
    "            grad_w = A.T.dot(z0) \n",
    "            return grad_w / n_samples\\\n",
    "               + tau * (theta - rho + u)\n",
    "\n",
    "        theta = optimize.minimize(\n",
    "            logistic_loss, theta, jac=logistic_grad).x\n",
    "        rho = np.fmax(theta + u - lam /tau, 0)\\\n",
    "                - np.fmax(-lam/tau - theta - u, 0)\n",
    "        u = u + theta - rho\n",
    "        obj_fun_history.append(logistic_loss(theta)\\\n",
    "                               + lam * np.sum(np.abs(theta)))\n",
    "\n",
    "    return theta, obj_fun_history\n",
    "    \n",
    "theta, func_vals = ADMM_logloss(A, b, lam=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG4tJREFUeJzt3X+Q3HWd5/Hnq7unAzOJkGRgjARIxCiGX4pjROXuBvFH\ncNdjz91y4c5zj9JKcQWnXt3WybpX5V35D1dleeLJGnOaBW+V3NZq7nKaFUG3jewhJEAwJCTrmLAh\nERIgSJwEMr/e90d/J3Rm+tc03ZmZ7/f1qArT/f3R3897El796Xd/+9uKCMzMLDtyMz0AMzM7vRz8\nZmYZ4+A3M8sYB7+ZWcY4+M3MMsbBb2aWMQ5+M7OMcfCbmWWMg9/MLGMKMz2Aanp7e2PZsmUt7Xvs\n2DF6enraO6BZLos1QzbrzmLNkM26p1vzI4888nxEnNPMtrMy+JctW8a2bdta2rdUKjEwMNDeAc1y\nWawZsll3FmuGbNY93Zol/WOz27rVY2aWMQ5+M7OMcfCbmWWMg9/MLGMc/GZmGePgNzPLmIbBL2m9\npMOSnqix/mJJD0o6IelPJ617StIOSdsltXZ+ppmZtVUzM/67gNV11h8BPg18qcb6ayLibRHRP82x\nTdtXf/Irdjw32unDmJnNaQ2DPyK2UA73WusPR8RWYKSdA2vF2p/9mideGJvpYZiZzWqd/uRuAPdL\nGgO+ERHram0oaQ2wBqCvr49SqTTtgynGePlEtLTvXDY0NJS5miGbdWexZshm3Z2sudPBf3VEHJR0\nLnCfpN3JK4gpkieFdQD9/f3Rysezz3zgfnKFMX+0OyOyWHcWa4Zs1t3Jmjt6Vk9EHEx+HgY2Aqs6\nebxiXoyNd/IIZmZzX8eCX1KPpAUTt4EPAlXPDGqXrkKO0fHo5CHMzOa8hq0eSfcAA0CvpAPAF4Au\ngIhYK+n1wDbgdcC4pM8CK4FeYKOkieN8NyJ+1IkiJnTlc4w6983M6moY/BFxY4P1zwJLq6w6ClzR\n4rha0pXPMeqTeszM6krVJ3fd4zczayxVwV9u9bjXY2ZWT/qC3zN+M7O6UhX8hbwc/GZmDaQq+Iv5\nHGPu9JiZ1ZWq4C+3epz8Zmb1pCv4Czmf1WNm1kC6gj8vf4DLzKyBVAV/0Wf1mJk1lKrg78rnGHOP\n38ysrtQFv1s9Zmb1pSv4Cz6P38yskVQFv3v8ZmaNpSr4u/I5AtznNzOrI1XBX8gLgBGfzG9mVlOq\ngr+YL5cz7OA3M6spVcHflQT/iBv9ZmY1pTP4faU2M7OaUhb87vGbmTXSMPglrZd0WNITNdZfLOlB\nSSck/emkdasl7ZE0KOm2dg26lmJhYsbv4Dczq6WZGf9dwOo6648Anwa+VLlQUh64E7gOWAncKGll\na8Nsjls9ZmaNNQz+iNhCOdxrrT8cEVuBkUmrVgGDEbE3IoaBDcD1r2Wwjbwa/J7xm5nV0ske/3nA\n0xX3DyTLOmaix+/TOc3MaivM9AAmSFoDrAHo6+ujVCpN+zF2vTAGwNZtj3J0b76dw5vVhoaGWvp9\nzXVZrDuLNUM26+5kzZ0M/oPA+RX3lybLqoqIdcA6gP7+/hgYGJj2AXueOgJbH+SSy67g6hW9095/\nriqVSrTy+5rrslh3FmuGbNbdyZo72erZCqyQtFxSEbgB2NTB41HI+XROM7NGGs74Jd0DDAC9kg4A\nXwC6ACJiraTXA9uA1wHjkj4LrIyIo5JuBe4F8sD6iNjZmTLKunzJBjOzhhoGf0Tc2GD9s5TbONXW\nbQY2tza06fN5/GZmjaXsk7sOfjOzRlIW/EmP39+/aGZWU6qCf+KyzCPjnvGbmdWSquD3ZZnNzBpL\nV/AXfK0eM7NG0hX8vmSDmVlD6Qr+nM/qMTNrJFXBn8uJvBz8Zmb1pCr4AfI59/jNzOpJX/ALhn1W\nj5lZTakL/kLOrR4zs3rSF/ySg9/MrI70Bb97/GZmdaUu+PNu9ZiZ1ZW64C/4dE4zs7rSF/w5udVj\nZlZHCoPfM34zs3pSF/w+j9/MrL7UBb9n/GZm9aUw+N3jNzOrp2HwS1ov6bCkJ2qsl6SvShqU9EtJ\nV1ase0rSDknbJW1r58Br8UXazMzqa2bGfxewus7664AVyZ81wNcnrb8mIt4WEf0tjXCaCjlfj9/M\nrJ6GwR8RW4AjdTa5Hvh2lP0COFvSknYNcLr8AS4zs/oKbXiM84CnK+4fSJY9AwRwv6Qx4BsRsa7W\ng0haQ/kVA319fZRKpdZGMzbK0LHx1vefg4aGhjJV74Qs1p3FmiGbdXey5nYEfz1XR8RBSecC90na\nnbyCmCJ5UlgH0N/fHwMDAy0d8O6d95Ir5Gh1/7moVCplqt4JWaw7izVDNuvuZM3tOKvnIHB+xf2l\nyTIiYuLnYWAjsKoNx6vLp3OamdXXjuDfBHwiObvnKuCliHhGUo+kBQCSeoAPAlXPDGqnvHw6p5lZ\nPQ1bPZLuAQaAXkkHgC8AXQARsRbYDHwYGASOAzclu/YBGyVNHOe7EfGjNo9/Cs/4zczqaxj8EXFj\ng/UB3FJl+V7gitaH1ppCDkbHg/HxIJfT6T68mdmsl75P7iZZPzLuWb+ZWTWpC/58Mst3n9/MrLrU\nBf/JGb+v0GlmVlXqgj+fVOQ3eM3Mqktd8BeSiny9HjOz6lIY/O7xm5nVk77gn+jxe8ZvZlZV6oLf\nPX4zs/pSF/yFk8HvVo+ZWTXpC35N9Pg94zczqyZ9wT8x4/d5/GZmVaUu+PM+ndPMrK7UBf+rZ/W4\nx29mVk36gj/nHr+ZWT2pC/68z+M3M6srdcF/8pINfnPXzKyq1Aa/e/xmZtWlLvjz7vGbmdWVuuD3\ntXrMzOprGPyS1ks6LOmJGusl6auSBiX9UtKVFetWS9qTrLutnQOvxa0eM7P6mpnx3wWsrrP+OmBF\n8mcN8HUASXngzmT9SuBGSStfy2CbUfBF2szM6moY/BGxBThSZ5PrgW9H2S+AsyUtAVYBgxGxNyKG\ngQ3Jth2Vk8jJwW9mVks7evznAU9X3D+QLKu1vOO68jlfssHMrIbCTA9ggqQ1lFtF9PX1USqVWnqc\noaEhcoh9T+2nVDrUxhHOXkNDQy3/vuayLNadxZohm3V3suZ2BP9B4PyK+0uTZV01llcVEeuAdQD9\n/f0xMDDQ0mBKpRJnzhuhb8kSBgYubekx5ppSqUSrv6+5LIt1Z7FmyGbdnay5Ha2eTcAnkrN7rgJe\niohngK3ACknLJRWBG5JtO66Qk3v8ZmY1NJzxS7oHGAB6JR0AvkB5Nk9ErAU2Ax8GBoHjwE3JulFJ\ntwL3AnlgfUTs7EANU7jHb2ZWW8Pgj4gbG6wP4JYa6zZTfmI4rYqFnM/jNzOrIXWf3AXoysvfwGVm\nVkNKgz/nHr+ZWQ2pDX73+M3Mqktl8BfzOUbd4zczqyqVwd9V8OmcZma1pDP43eM3M6sptcE/7FaP\nmVlVqQz+omf8ZmY1pTL4u/Lu8ZuZ1ZLK4C/kc/4Al5lZDakMfvf4zcxqS2XwF93qMTOrKZXB79M5\nzcxqS2fwFxz8Zma1pDP48+XLMpevGG1mZpVSGfzFvAAYHXfwm5lNlsrg78qXy3K7x8xsqnQH/6hn\n/GZmk6Uz+AvlsnxNfjOzqVIZ/BM9frd6zMymair4Ja2WtEfSoKTbqqxfKGmjpF9KeljSpRXrnpK0\nQ9J2SdvaOfhaCjn3+M3Maik02kBSHrgT+ABwANgqaVNE7KrY7PPA9oj4F5IuTra/tmL9NRHxfBvH\nXddEq8fBb2Y2VTMz/lXAYETsjYhhYANw/aRtVgI/BYiI3cAySX1tHek0TLR6hv3mrpnZFA1n/MB5\nwNMV9w8A75q0zePAR4GfS1oFXAgsBQ4BAdwvaQz4RkSsq3YQSWuANQB9fX2USqVplPGqoaEhBg/v\nBOAXW7dy+Kx8S48zlwwNDbX8+5rLslh3FmuGbNbdyZqbCf5m3A7cIWk7sAN4DBhL1l0dEQclnQvc\nJ2l3RGyZ/ADJE8I6gP7+/hgYGGhpIKVSiSvffAk8+jCXX/F2+pctaulx5pJSqUSrv6+5LIt1Z7Fm\nyGbdnay5meA/CJxfcX9psuykiDgK3AQgScA+YG+y7mDy87CkjZRbR1OCv50mzuP36ZxmZlM10+Pf\nCqyQtFxSEbgB2FS5gaSzk3UAnwK2RMRRST2SFiTb9AAfBJ5o3/CrKxaSSzb4mvxmZlM0nPFHxKik\nW4F7gTywPiJ2Sro5Wb8WeCtwt6QAdgKfTHbvAzaWXwRQAL4bET9qfxmn8iUbzMxqa6rHHxGbgc2T\nlq2tuP0g8OYq++0FrniNY5w2B7+ZWW2p/OTuqz1+t3rMzCZLZfAXT16kzTN+M7PJUhn8BV+rx8ys\nplQGv3v8Zma1pTL4i+7xm5nVlMrg7yq41WNmVks6g99v7pqZ1ZTK4C/kPOM3M6sllcEviWI+5x6/\nmVkVqQx+gK68GPWM38xsivQGfyHnVo+ZWRXpDX63eszMqkpt8J/RlePl4dGZHoaZ2ayT2uBf1F3k\nyPGRmR6Gmdmsk9rgX9hT5MVjwzM9DDOzWSe1wb+op8gRB7+Z2RTpDf7uIi8cOzHTwzAzm3XSG/zz\ni7wyMs7Lw2MzPRQzs1kltcG/uKf83e+e9ZuZnaqp4Je0WtIeSYOSbquyfqGkjZJ+KelhSZc2u2+n\nLOwuB/+Lx3xmj5lZpYbBLykP3AlcB6wEbpS0ctJmnwe2R8TlwCeAO6axb0csnu8Zv5lZNc3M+FcB\ngxGxNyKGgQ3A9ZO2WQn8FCAidgPLJPU1uW9HnJzxH/eZPWZmlZoJ/vOApyvuH0iWVXoc+CiApFXA\nhcDSJvftiMU98wB4YcjBb2ZWqdCmx7kduEPSdmAH8BgwrdNpJK0B1gD09fVRKpVaGsjQ0BClUonx\nCHKC7U8OUhrb39JjzRUTNWdNFuvOYs2Qzbo7WXMzwX8QOL/i/tJk2UkRcRS4CUCSgH3AXuDMRvtW\nPMY6YB1Af39/DAwMNFXAZKVSiYl9F/39/SzoPZeBgctbeqy5orLmLMli3VmsGbJZdydrbqbVsxVY\nIWm5pCJwA7CpcgNJZyfrAD4FbEmeDBru20mL/eldM7MpGs74I2JU0q3AvUAeWB8ROyXdnKxfC7wV\nuFtSADuBT9bbtzOlTLWwp8vBb2Y2SVM9/ojYDGyetGxtxe0HgTc3u+/psrhnHk8+e3QmDm1mNmul\n9pO7UJ7x+wqdZmanSnXwL+qZx29fHmFs3N/EZWY2Id3B391FBPzWH+IyMzsp3cE/v/whLr/Ba2b2\nqlQH/8QVOh38ZmavSnXwT1yvx8FvZvaqVAf/xBU6j7jHb2Z2UqqD/+zuLgCO+EJtZmYnpTr45xXy\nLJhX8IzfzKxCqoMfYKGv12NmdorUB/8iB7+Z2SlSH/y+QqeZ2alSH/xu9ZiZnSr1wT8x44/w9XrM\nzCADwb+wp8iJ0XGOD0/rmyDNzFIr9cG/yJdtMDM7RfqD35dtMDM7RfqD35dtMDM7RfqDf2LG78s2\nmJkBWQj+ZMb/omf8ZmZAk8EvabWkPZIGJd1WZf1Zkv6vpMcl7ZR0U8W6pyTtkLRd0rZ2Dr4ZC+YV\n6MqLF9zjNzMDoNBoA0l54E7gA8ABYKukTRGxq2KzW4BdEfERSecAeyR9JyIm0vaaiHi+3YNvhiQW\ndhfd6jEzSzQz418FDEbE3iTINwDXT9omgAWSBMwHjgCjbR3pa9A7fx7PHH1lpodhZjYrqNEnWiX9\nEbA6Ij6V3P/XwLsi4taKbRYAm4CLgQXAH0fED5N1+4CXgDHgGxGxrsZx1gBrAPr6+t6xYcOGlgoa\nGhpi/vz5pyy7e+cJHvzNKF+7tptCTi097mxWreYsyGLdWawZsln3dGu+5pprHomI/ma2bdjqadKH\ngO3A+4CLgPsk/TwijgJXR8RBSecmy3dHxJbJD5A8IawD6O/vj4GBgZYGUiqVmLzv8cXP8HffeZSF\nF13BOy5c1NLjzmbVas6CLNadxZohm3V3suZmWj0HgfMr7i9NllW6Cfh+lA0C+yjP/omIg8nPw8BG\nyq2j0+rdb1yMBH8/+MLpPrSZ2azTTPBvBVZIWi6pCNxAua1TaT9wLYCkPuAtwF5JPUkbCEk9wAeB\nJ9o1+GYt7CmycsnreGBwRt5fNjObVRoGf0SMArcC9wJPAn8dETsl3Szp5mSzLwLvkbQD+AnwueQs\nnj7gAUmPAw8DP4yIH3WikEbe+6ZeHtv/IseHZ817zmZmM6KpHn9EbAY2T1q2tuL2byjP5ifvtxe4\n4jWOsS3e+6Ze1m3Zy9anXuSfvfmcmR6OmdmMSf0ndye8c9lCuvLi/7ndY2YZl5ng7y4WePsFC93n\nN7PMy0zwA7z3ol52PXOUF335BjPLsEwF/9UrFhMBD+71aZ1mll2ZCv7Ll55NTzHPz3/13EwPxcxs\nxmQq+LvyOd6/so8NW5/muw/tn+nhmJnNiHZdsmHOuP2jl/PSyyN8fuMOjhw7wS3XvInyteXMzLIh\nUzN+gDOLef7HJ/r5g7e9gS/9+B/4z5t2Mjo2PtPDMjM7bTI344dyy+fLH3sbvfPn8c0H9rHn0O/4\n2r+8kt7582Z6aGZmHZe5Gf+EXE78p99fyZc/dgWP7f8tH/nvD/DY/hdnelhmZh2X2eCf8NErl/K9\nf/se8jnxsW88yNqf/Zrx8frfUWBmNpdlPvgBLj3vLH7w767m/W/t4/a/3c3Hv/UQz77kb+wys3Ry\n8CfO7i7yF//qSv7rH17GY/t/ywf+28/41gP7GPEbv2aWMg7+CpL443dewObP/BPefsFCvviDXXzo\nK1u4f9chGn1FpZnZXOHgr2J5bw933/RO1v+bfgj41Le38aGvbOF/bd3PKyNjMz08M7PXJJOnczZD\nEu+7uI+r33QOmx7/Dd96YB+f+94Obv/b3Vx32RJ+77IlvGv5Igp5P3ea2dzi4G+gWMjxR+9Yyh9e\neR4P7n2B7zy0n42PHuS7D+1nUU+Rq964iKveuJhVyxfxpnPm+4nAzGY9B3+TJPGei3p5z0W9vDw8\nRmnPYe7bdYhf7H2BzTueBWBeIcdbXr+Ai1+/gOW987lwcTcXLOpmyVlnsKin6EtDmNms4OBvwZnF\nPNddtoTrLltCRHDgxZd5eN8RnnzmKE8+e5Sf7j7M80MHTtmnmM9xzoJ59M4vsrCnyKLuIq87s4sF\nZxSYP69Az7wC3cU83cU8Z3SV/8wr5JhXyFMs5Cjmc3QVRCGXoysvCvkchZzI5+Q3ns1sWpoKfkmr\ngTuAPPDNiLh90vqzgL8CLkge80sR8ZfN7DvXSeL8Rd2cv6j7lOW/e2WE/UeO8/SR4zz70isc+t0J\nDr30Ci8cG+bIsWF+dWiI370ywtCJUdrxeTH9+IfkJXI5lX+q/OnkXHJbEhOvN5TczwlyyfKJVyMT\nL0okmNijfPvUmk85do077Xh9E8l/Jv+KBBx/+Tjdj5TacJTOqvz91lP5BB4n/3Oq48en1vzq32vz\nv/GImHoMNT/Wdml2HMePHaf70Z+15xgVj1+t0pmcRi3qLvLXN7+748dpGPyS8sCdwAeAA8BWSZsi\nYlfFZrcAuyLiI5LOAfZI+g4w1sS+qbTgjC4uecNZXPKGs+puFxEcGx7j+IlRjg+PcXx4jJdHxjgx\nOsaJkXFOjI4xPBYMj44zMjbO6Ng4I2PByNg4YxGMjQWDe/dxwYUXMjYejEUQQfl28owyHsF4EiqR\nhGhE+dhj4/Hq/Yl/8hVBe/J/mpP3J41/Ui3Vlk/ZYZq5Mvl/0kiOdfjwK5x77utq7tfCoWru08pj\nTew3rSSZ9MQ5OYQPHTq15phyY3rHmjhG1QA+XZoYx+HnXubccxa85mOcfOjKf+9VN5+ZtuyCM05P\nE6aZo6wCBiNiL4CkDcD1QGV4B7BA5X+l84EjwCjwrib2zTRJzJ9Xbve0qpQ/yMDAW9o4qrmhVCox\nMHDlTA/jtMpizZDdujulmVNQzgOerrh/IFlW6WvAW4HfADuAz0TEeJP7mpnZadSu1xUfArYD7wMu\nAu6T9PPpPICkNcAagL6+PkqlUksDGRoaannfuSqLNUM2685izZDNujtZczPBfxA4v+L+0mRZpZuA\n26Pc5B2UtA+4uMl9AYiIdcA6gP7+/hgYGGhm/FOUXxK2tu9clcWaIZt1Z7FmyGbdnay5mVbPVmCF\npOWSisANwKZJ2+wHrgWQ1Ae8Bdjb5L5mZnYaNZzxR8SopFuBeymfkrk+InZKujlZvxb4InCXpB2U\n3zz/XEQ8D1Bt386UYmZmzWiqxx8Rm4HNk5atrbj9G+CDze5rZmYzxxeWMTPLGAe/mVnGaDZe50XS\nc8A/trh7L/B8G4czF2SxZshm3VmsGbJZ93RrvjAizmlmw1kZ/K+FpG0R0T/T4zidslgzZLPuLNYM\n2ay7kzW71WNmljEOfjOzjElj8K+b6QHMgCzWDNmsO4s1Qzbr7ljNqevxm5lZfWmc8ZuZWR2pCX5J\nqyXtkTQo6baZHk+nSDpf0t9J2iVpp6TPJMsXSbpP0q+SnwtneqztJikv6TFJP0juZ6HmsyX9jaTd\nkp6U9O601y3p3yf/tp+QdI+kM9JYs6T1kg5LeqJiWc06Jf1Zkm97JH3otRw7FcFf8S1h1wErgRsl\nrZzZUXXMKPAfImIlcBVwS1LrbcBPImIF8JPkftp8Bniy4n4War4D+FFEXAxcQbn+1NYt6Tzg00B/\nRFxK+RpfN5DOmu8CVk9aVrXO5P/xG4BLkn3+Ism9lqQi+Kn4lrCIGAYmvukrdSLimYh4NLn9O8pB\ncB7leu9ONrsb+IOZGWFnSFoK/B7wzYrFaa/5LOCfAt8CiIjhiPgtKa+b8jXEzpRUALopf8FT6mqO\niC2Uv62wUq06rwc2RMSJiNgHDFLOvZakJfgz+U1fkpYBbwceAvoi4plk1bNA3wwNq1O+AvxHYLxi\nWdprXg48B/xl0uL6pqQeUlx3RBwEvkT5Uu/PAC9FxI9Jcc2T1KqzrRmXluDPHEnzge8Bn42Io5Xr\nki/ESc3pWpJ+HzgcEY/U2iZtNScKwJXA1yPi7cAxJrU40lZ30tO+nvKT3huAHkkfr9wmbTXX0sk6\n0xL8TX/TVxpI6qIc+t+JiO8niw9JWpKsXwIcnqnxdcB7gX8u6SnKbbz3Sfor0l0zlGd1ByLioeT+\n31B+Ikhz3e8H9kXEcxExAnwfeA/prrlSrTrbmnFpCf7MfNOXJFHu+T4ZEV+uWLUJ+JPk9p8A/+d0\nj61TIuLPImJpRCyj/Hf704j4OCmuGSAingWelvSWZNG1wC7SXfd+4CpJ3cm/9Wspv4+V5por1apz\nE3CDpHmSlgMrgIdbPkpEpOIP8GHgH4BfA38+0+PpYJ1XU37590vKX3C/Pal9MeWzAH4F3A8smumx\ndqj+AeAHye3U1wy8DdiW/H3/b2Bh2usG/guwG3gC+J/AvDTWDNxD+X2MEcqv7j5Zr07gz5N82wNc\n91qO7U/umpllTFpaPWZm1iQHv5lZxjj4zcwyxsFvZpYxDn4zs4xx8JuZZYyD38wsYxz8ZmYZ8/8B\nvI5WsFzyyYYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11ac045c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(func_vals)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHVdJREFUeJzt3X+UXWV97/H3xyHgEH6kCowyCSbaNEhFCJkSLFw6USBB\nu0z8sUq4AorSrNyWIreaNtFeq6v1go16q6uRGLlYXSpZrpqkWRoZFTiNij+SGCQkMjU3oGQiRMSk\nDI4wid/7x94DJ8NM5jmH2XNmzv681sqas5/9POd8v2dOznf2s38pIjAzMxvJCxodgJmZTQwuGGZm\nlsQFw8zMkrhgmJlZEhcMMzNL4oJhZmZJXDDMzCyJC4aZmSVxwTAzsyTHNDqA0XTKKafE9OnT6xr7\n5JNPMnny5NENaJwrY85QzrzLmDOUM+9ac962bdtjEXFqSt+mKhjTp09n69atdY2tVCp0dnaObkDj\nXBlzhnLmXcacoZx515qzpJ+l9vWUlJmZJXHBMDOzJC4YZmaWxAXDzMySuGCYmVkSFwwzM0vigmFm\nZkkKLRiSFkjqlrRb0vIh1i+UdJ+keyVtlXRR1bqHJO0YWFdknGZmNrLCTtyT1AKsAi4F9gJbJG2M\niF1V3e4ENkZESHo18GXgzKr18yLisaJiNDOzdEVuYZwP7I6IPRHxNLAWWFjdISJ6IyLyxclAYGZm\n45Ke/b4e5SeW3gosiIjr8uWrgbkRcf2gfm8CbgJOA94QEd/L2x8EDgKHgU9HxJphXmcJsASgra1t\nztq1a+uKt7e3lxNOOKGusRNVGXOGcuZdxpyhnHnXmvO8efO2RURHUueIKOQf8Fbg1qrlq4F/OUr/\ni4FvVS235z9PA34MXDzSa86ZMyfqdffdd9c9dqIqY84R5cy7jDlHlDPvWnMGtkbi93qRU1I9wLSq\n5al525AiYjPwckmn5Ms9+c/9wHqyKS4zM2uQIgvGFmCmpBmSjgUWAxurO0j6fUnKH58HHAf8StJk\nSSfm7ZOBy4D7C4zVzJrIhu09XHjzXezoOciFN9/Fhu3D/q1qNSjsKKmIOCTpeqALaAFui4idkpbm\n61cDbwGukdQP9AFXRERIagPW57XkGOBLEXFHUbGaWfPYsL2HFet20Nd/GKZBz4E+VqzbAcCi2e0N\njm5iK/R+GBGxCdg0qG111eOPAB8ZYtwe4JwiYzOz5rSyqzsrFlX6+g+zsqvbBeN58pneZtZU9h3o\nq6nd0rlgmFlTOX1Ka03tls4Fw8yayrL5s2id1HJEW+ukFpbNn9WgiJpHU93T28xsYD/Fyq5u4Ana\np7SybP4s778YBS4YZtZ0Fs1uZ9HsdiqVCn/1ts5Gh9M0PCVlZmZJXDDMzCyJC4aZmSVxwTAzsyQu\nGGZmlsQFw8zMkrhgmJlZEhcMMzNL4oJhZmZJXDDMzCyJC4aZmSVxwTAzsyQuGGZmlqTQgiFpgaRu\nSbslLR9i/UJJ90m6V9JWSReljjUzs7FVWMGQ1AKsAi4HzgKulHTWoG53AudExLnAO4FbaxhrZmZj\nqMgtjPOB3RGxJyKeBtYCC6s7RERvRES+OBmI1LFmZja2iryBUjvwcNXyXmDu4E6S3gTcBJwGvKGW\nsfn4JcASgLa2NiqVSl3B9vb21j12oipjzlDOvMuYM5Qz7yJzbvgd9yJiPbBe0sXAPwCX1Dh+DbAG\noKOjIzo7O+uKo1KpUO/YiaqMOUM58y5jzlDOvIvMucgpqR5gWtXy1LxtSBGxGXi5pFNqHWtmZsUr\nsmBsAWZKmiHpWGAxsLG6g6Tfl6T88XnAccCvUsaamdnYKmxKKiIOSboe6AJagNsiYqekpfn61cBb\ngGsk9QN9wBX5TvAhxxYVq5mZjazQfRgRsQnYNKhtddXjjwAfSR1rZmaN4zO9zcwsiQuGmZklccEw\nM7MkLhhmZpbEBcPMzJK4YJiZWRIXDDMzS+KCYWZmSVwwzMwsiQuGmZklccEwM7MkLhhmZpbEBcPM\nzJK4YJiZWRIXDDMzS+KCYWZmSVwwzMwsiQuGmZklKbRgSFogqVvSbknLh1j/Nkn3Sdoh6R5J51St\neyhvv1fS1iLjNDOzkRV2T29JLcAq4FJgL7BF0saI2FXV7UHgTyLi15IuB9YAc6vWz4uIx4qK0czM\n0hW5hXE+sDsi9kTE08BaYGF1h4i4JyJ+nS9+H5haYDxmZvY8KCKKeWLprcCCiLguX74amBsR1w/T\n/73AmVX9HwQOAoeBT0fEmmHGLQGWALS1tc1Zu3ZtXfH29vZywgkn1DV2oipjzlDOvMuYM5Qz71pz\nnjdv3raI6EjpW9iUVC0kzQPeBVxU1XxRRPRIOg34pqQHImLz4LF5IVkD0NHREZ2dnXXFUKlUqHfs\nRFXGnKGceZcxZyhn3kXmXOSUVA8wrWp5at52BEmvBm4FFkbErwbaI6In/7kfWE82xWVmZg1SZMHY\nAsyUNEPSscBiYGN1B0lnAOuAqyPiP6vaJ0s6ceAxcBlwf4GxmpnZCAqbkoqIQ5KuB7qAFuC2iNgp\naWm+fjXwAeDFwKckARzK59LagPV52zHAlyLijqJiNTOzkRW6DyMiNgGbBrWtrnp8HXDdEOP2AOcM\nbjczs8bxmd5mZpbEBcPMzJK4YJiZWRIXDDMzS+KCYWZmSVwwzMwsiQuGmZklccEwM7MkLhhmZpbE\nBcPMzJK4YJiZWRIXDDMzS+KCYWZmSZIKhjJXSfpAvnyGJN/QyMysRFK3MD4FvAa4Ml9+AlhVSERm\nZjYupd4PY25EnCdpO0BE/Dq/i56ZmZVE6hZGv6QWIAAknQr8rrCozMxs3EktGJ8E1gOnSfow8B3g\nf480SNICSd2SdktaPsT6t0m6T9IOSfdIOid1rJmZja2kKamI+KKkbcDrAAGLIuInRxuTb5GsAi4F\n9gJbJG2MiF1V3R4E/iSf4rocWAPMTRxrZmZj6KgFQ9KLqhb3A7dXr4uIx48y/Hxgd35/biStBRYC\nz3zpR8Q9Vf2/D0xNHWtmZmNrpC2MbWT7LQScAfw6fzwF+Dkw4yhj24GHq5b3AnOP0v9dwNfrHGtm\nZgU7asGIiBkAkj4DrI+ITfny5cCi0QpC0jyygnFRHWOXAEsA2traqFQqdcXQ29tb99iJqow5Qznz\nLmPOUM68i8w59bDaCyLizwcWIuLrkv5phDE9wLSq5al52xEkvRq4Fbg8In5Vy9g8ljVk+z7o6OiI\nzs7OEcIaWqVSod6xE1UZc4Zy5l3GnKGceReZc+pRUvsk/Z2k6fm/9wP7RhizBZgpaUZ+zsZiYGN1\nB0lnAOuAqyPiP2sZa2ZmYyt1C+NK4O/JDq0F2MyzZ30PKSIOSboe6AJagNsiYqekpfn61cAHgBcD\nn5IEcCgiOoYbW1tqZmY2mlIPq30ceHetT57v89g0qG111ePrgOtSx5qZWeMkFQxJd5Of5V0tIl47\n6hGZmdm4lDol9d6qxy8E3gIcGv1wzMxsvEqdkto2qOm7kn5YQDxmZjZOpU5JVZ/x/QJgDnByIRGZ\nmdm4lDolVX3G9yGya0C9q6igzMxs/EktGK+MiN9WN0g6roB4zMxsnEo9ce+eIdq+N5qBmJnZ+DbS\n1WpfQnYhwFZJs8mmpABOAo4vODYzMxtHRpqSmg+8g+xaTh+van8CeF9BMZmZ2Tg00tVqPwd8TtJb\nIuIrYxSTmZmNQyNNSV0VEV8Apkv668HrI+LjQwwzM7MmNNKU1OT85wlFB2JmVhYbtvewsqubfQf6\nOH1KK8vmz2LR7PZGhzWikaakPp3//NDYhGNm1tw2bO9hxbod9PUfBqDnQB8r1u0AGPdFI/VM71OB\nPwemV4+JiHcWE5aZWXNa2dX9TLEY0Nd/mJVd3c1RMIB/B74NfAs4PEJfMzMbxr4DfTW1jyepBeP4\niPjbQiMxMyuB06e00jNEcTh9SmsDoqlN6pneX5X0+kIjMTMrgWXzZ9E6qeWIttZJLSybP6tBEaVL\n3cJ4N/A+SU8B/WRnfEdEnFRYZGZmTWhgP0XTHSU1ICJOrOfJJS0APkF2X+5bI+LmQevPBD4LnAe8\nPyI+WrXuIbIzyg+T3+u7nhjMzMabRbPbJ0SBGCz1KKnzhmg+CPwsIoa8856kFmAVcCmwF9giaWNE\n7Krq9jhwA7BomJeeFxGPpcRoZmbFSp2S+hTZVsCOfPls4H7gZEn/IyK+McSY84HdEbEHQNJaYCHw\nTMGIiP3AfklvqDN+MzMbI6k7vfcBsyNiTkTMAc4F9pBtPfzTMGPagYerlvfmbakC+JakbZKW1DDO\nzMwKkLqF8QcRsXNgISJ2STozIvZIOtq45+OiiOiRdBrwTUkPRMTmwZ3yYrIEoK2tjUqlUteL9fb2\n1j12oipjzlDOvMuYM5Qz7yJzTi0YOyXdAqzNl68AduV33esfZkwPMK1qeWreliQievKf+yWtJ5vi\nek7BiIg1wBqAjo6O6OzsTH2JI1QqFeodO1GVMWcoZ95lzBnKmXeROadOSb0D2A3cmP/bk7f1A/OG\nGbMFmClphqRjgcXAxpQXkzRZ0okDj4HLyPaZmJlZg6QeVtsHfCz/N1jvMGMOSboe6CI7rPa2iNgp\naWm+fnV+R7+tZHfw+52kG4GzgFOA9fl01zHAlyLijpoyMzOzUZV6WO1M4CayL/MXDrRHxMuPNi4i\nNgGbBrWtrnr8CNlU1WD/BZyTEpuZmY2N1CmpzwK3AIfIpqA+D3yhqKDMzGz8SS0YrRFxJ6CI+FlE\nfBDwuRNmZiWSepTUU5JeAPw03y/Rg+/CZ2ZWKqlbGO8Gjie7jMcc4Grg7UUFZWZm40/qUVJb8oe9\nwLXFhWNmZuPVUQuGpKOeNxERbxzdcMzMbLwaaQvjNWTXg7od+AHZfTDMzKyERioYLyG7wOCVwH8H\nvgbcXn1dKTMzK4ej7vSOiMMRcUdEvB24gOzyIJX8SCkzMyuREXd65xcYfAPZVsZ04JPA+mLDMjOz\n8Waknd6fB15FdnmPD0WELwBoZlZSI21hXAU8SXYexg1V974QEBFxUoGxmZnZOHLUghERqSf2mZlZ\nk3NBMDOzJC4YZmaWxAXDzMySuGCYmVkSFwwzM0tSaMGQtEBSt6TdkpYPsf5MSd+T9JSk99Yy1szM\nxlZhBUNSC7AKuJzsXuBXSjprULfHye6x8dE6xpqZ2RhKveNePc4HdkfEHgBJa4GFwK6BDhGxH9gv\nafDtXkcca2bltGF7Dyu7utl3oI/Tp7SybP4sFs1ub3RYpVDklFQ72aXRB+zN24oea2ZNasP2Hlas\n20HPgT4C6DnQx4p1O9iwvafRoZVCkVsYY0LSEmAJQFtbG5VKpa7n6e3trXvsRFXGnKGceTdLzo8+\n8gR/cebvBrUe4tHuH1E5+NPn9G+WvGtRZM5FFoweYFrV8tS8bVTHRsQaYA1AR0dHdHZ21hwoQKVS\nod6xE1UZc4Zy5t0sOV+7/GvEEBMjAh68ufM57c2Sdy2KzLnIKaktwExJMyQdCywGjnrL11Eaa2ZN\n6vQprTW12+gqrGBExCHgeqAL+Anw5YjYKWmppKUAkl4iaS/w18DfSdor6aThxhYVq5lNDMvmz6J1\nUssRba2TWlg2f1aDIiqXQvdhRMQmsntpVLetrnr8CNl0U9JYMyu3gaOhfJRUY0z4nd5mVi6LZre7\nQDSILw1iZmZJXDDMzCyJC4aZmSVxwTAzsyQuGGZmlsQFw8zMkrhgmJlZEhcMMzNL4hP3zMyG4Xtv\nHMkFw8xsCAP33ujrPww8e+8NoLRFw1NSZmZDWNnV/UyxGNDXf5iVXd0NiqjxXDDMzIaw70BfTe1l\n4IJhZjYE33vjuVwwzMyG4HtvPJd3epuZDcH33nguFwwzs2H43htH8pSUmZklKbRgSFogqVvSbknL\nh1gvSZ/M198n6byqdQ9J2iHpXklbi4zTzMxGVtiUlKQWYBVwKbAX2CJpY0Tsqup2OTAz/zcXuCX/\nOWBeRDxWVIxmZpauyC2M84HdEbEnIp4G1gILB/VZCHw+Mt8Hpkh6aYExmZlZnYosGO3Aw1XLe/O2\n1D4BfEvSNklLCovSzMySjOejpC6KiB5JpwHflPRARGwe3CkvJksA2traqFQqdb1Yb29v3WMnqjLm\nDOXMu4w5QznzLjLnIgtGDzCtanlq3pbUJyIGfu6XtJ5sius5BSMi1gBrADo6OqKzs7OuYCuVCvWO\nnajKmDOUM+8y5gzlzLvInIucktoCzJQ0Q9KxwGJg46A+G4Fr8qOlLgAORsQvJE2WdCKApMnAZcD9\nBcZqZmYjKGwLIyIOSboe6AJagNsiYqekpfn61cAm4PXAbuA3wLX58DZgvaSBGL8UEXcUFauZmY2s\n0H0YEbGJrChUt62uehzAXw4xbg9wTpGxmZlZbXymt5mZJXHBMDOzJC4YZmaWZDyfh2HWNDZs7/Fl\nsm3Cc8EwK9iG7T2sWLfjmftD9xzoY8W6HQAuGjaheErKrGAru7qfKRYD+voPs7Kru0ERmdXHBcOs\nYPsO9NXUbjZeuWCU0IbtPVx4813s6DnIhTffxYbtg6/YYqPp9CmtNbWbjVcuGCUzMJ/ek/91OzCf\n7qJRnGXzZ9E6qeWIttZJLSybP6tBEZnVxzu9S+Zo8+neAVuMgfd1PB0l5aO2rB4uGCXj+fTGWDS7\nfdx8IfuoLauXp6RKxvPp5qO2rF4uGCXj+XTzVqbVy1NSJVM9nw5P0O7569I5fUrrMwc9DG43Oxpv\nYZTQotntfHf5azm7/WS+u/y1LhYl461Mq5e3MMxKZjwetWUTgwuGWQmNp6O2bOLwlJSZmSUptGBI\nWiCpW9JuScuHWC9Jn8zX3yfpvNSxZmY2tgqbkpLUAqwCLgX2AlskbYyIXVXdLgdm5v/mArcAcxPH\njoqBM14XT3uCGz/0DSQ48Jt+Tm6dVOjj06e0Mu/MU7n7gV+y70Bf4a83+LWXzZ/FFI4843e8xFf0\ne/nOV/SN+e+6eh/BcO95kY8H51wd0+Czvsfz7z31vUz5jNf6erU+z9Hey+fzvCn/r4uiiCjmiaXX\nAB+MiPn58gqAiLipqs+ngUpE3J4vdwOdwPSRxg6lo6Mjtm7dmhxj9Rmv7zn7EB/bUZ5dOq2TWvhf\nc4J/2KbnnMTV7Brxu26d1MJNbz4b4IizrMfKUDm3TmrhLXPa+cq2ngn1GUh9L0frM17U+1TE87ZO\nauGmP25h0eWXJo+RtC0iOlL6Fvm/ph14uGp5L9lWxEh92hPHPm9DnfFaFn39h3n8ycP09beM3Nme\nt+ozqcfLZ66v/zC3/+BhDhf0R2NRUt/L0fqMF/U+FfG8ff2HefRg/6g932AT/k9qSUuAJQBtbW1U\nKpXksYunPQHTssdtrdlfYWVyWglzhkb+rp/Ifkwb+1duvs932ntZxs/47x1LTd+DtSiyYPRw5K9z\nat6W0mdSwlgAImINsAayKanOzs7kAN9/813PnPFatikpgPeefZiP7SjfFkajftft+ZnUQ51lXbTh\ncm6RJtwWBqS/l6P1GS/qfSrieVec+zv+rIbvwVoUeZTUFmCmpBmSjgUWAxsH9dkIXJMfLXUBcDAi\nfpE49nkb6ozXsmid1MKLJk8qbf5jbeBM6vH0mWud1MKVc6eNm3hSpb6Xo/UZL+p9KuJ5Wye10Hby\nC0ft+QYr7M+siDgk6XqgC2gBbouInZKW5utXA5uA1wO7gd8A1x5t7GjHOPi6SlPG8ZE9hRwldfCn\n3PTmmaU7SgoOjfnvevCZ1GN9lNTgnKtj6njZiybsUVLDvZcpn/FaX6/6fRrNo6Rqfd6U/9dFKewo\nqUao9SipapVKhVqms5pBGXOGcuZdxpyhnHnXmnMtR0n5TG8zM0vigmFmZklcMMzMLIkLhpmZJXHB\nMDOzJC4YZmaWxAXDzMySuGCYmVkSFwwzM0vigmFmZklcMMzMLIkLhpmZJWmqiw9K+iXwszqHnwI8\nNorhTARlzBnKmXcZc4Zy5l1rzi+LiFNTOjZVwXg+JG1NvWJjsyhjzlDOvMuYM5Qz7yJz9pSUmZkl\nccEwM7MkLhjPWtPoABqgjDlDOfMuY85QzrwLy9n7MMzMLIm3MMzMLEnpC4akBZK6Je2WtLzR8RRF\n0jRJd0vaJWmnpHfn7S+S9E1JP81//l6jYx1tklokbZf01Xy5DDlPkfRvkh6Q9BNJr2n2vCX9z/yz\nfb+k2yW9sBlzlnSbpP2S7q9qGzZPSSvy77duSfOfz2uXumBIagFWAZcDZwFXSjqrsVEV5hDwnog4\nC7gA+Ms81+XAnRExE7gzX2427wZ+UrVchpw/AdwREWcC55Dl37R5S2oHbgA6IuJVQAuwmObM+V+B\nBYPahswz/z++GPjDfMyn8u+9upS6YADnA7sjYk9EPA2sBRY2OKZCRMQvIuJH+eMnyL5A2sny/Vze\n7XPAosZEWAxJU4E3ALdWNTd7zicDFwP/FyAino6IAzR53sAxQKukY4DjgX00Yc4RsRl4fFDzcHku\nBNZGxFMR8SCwm+x7ry5lLxjtwMNVy3vztqYmaTowG/gB0BYRv8hXPQK0NSisovwz8DfA76ramj3n\nGcAvgc/mU3G3SppME+cdET3AR4GfA78ADkbEN2jinAcZLs9R/Y4re8EoHUknAF8BboyI/6peF9kh\nc01z2JykPwX2R8S24fo0W865Y4DzgFsiYjbwJIOmYpot73zOfiFZsTwdmCzpquo+zZbzcIrMs+wF\noweYVrU8NW9rSpImkRWLL0bEurz5UUkvzde/FNjfqPgKcCHwRkkPkU03vlbSF2junCH7K3JvRPwg\nX/43sgLSzHlfAjwYEb+MiH5gHfDHNHfO1YbLc1S/48peMLYAMyXNkHQs2c6hjQ2OqRCSRDan/ZOI\n+HjVqo3A2/PHbwf+faxjK0pErIiIqRExnex3e1dEXEUT5wwQEY8AD0ualTe9DthFc+f9c+ACScfn\nn/XXke2na+acqw2X50ZgsaTjJM0AZgI/rPdFSn/inqTXk81ztwC3RcSHGxxSISRdBHwb2MGz8/nv\nI9uP8WXgDLIr/f5ZRAzeoTbhSeoE3hsRfyrpxTR5zpLOJdvRfyywB7iW7A/Eps1b0oeAK8iOCNwO\nXAecQJPlLOl2oJPsqrSPAn8PbGCYPCW9H3gn2ftyY0R8ve7XLnvBMDOzNGWfkjIzs0QuGGZmlsQF\nw8zMkrhgmJlZEhcMMzNL4oJhpSTpJZLWSvp/krZJ2iTpD+p4nhvyq8F+MT/W/VuS7pV0RX5JjmEv\nZinpjfVeITm/Gu1f1DPWrF4+rNZKJz+x6x7gcxGxOm87BzgpIr5d43M9AFwSEXslXQD8Y0RcMupB\nP/d1pwNfza/MajYmvIVhZTQP6B8oFgAR8WPgO5JW5vdT2CHpioH1kpZJ2iLpvvwEMSStBl4OfF3S\n3wJfAP4o38J4haSKpI687wJJP5L0Y0l35m3vkPQv+eNTJX0lf40tki7M2z+Y3/+gImmPpBvykG4G\nXpG/1kpJL5W0OV++X9J/K/pNtPI5ptEBmDXAq4ChLkj4ZuBcsvtHnAJskbQZOJvskgrnAwI2Sro4\nIpZKWgDMi4jHJP2A/GxygGxDJisGwGeAiyPiQUkvGuK1PwH8n4j4jqQzgC7glfm6M8mK3IlAt6Rb\nyC4m+KqIODd/jfcAXRHx4fx+B8c/nzfIbCguGGbPugi4PSIOk13M7T+APyK7t8RlZJebgOxyEzOB\nzYnPewGwOb8fAcNcmuIS4KyBIgOclF9ZGOBrEfEU8JSk/Qx9ie4twG35BSY3RMS9ibGZJXPBsDLa\nCby1hv4CboqITxcUD2TTwxdExG+PeOGsgDxV1XSYIf7fRsRmSReT3SzqXyV9PCI+X2C8VkLeh2Fl\ndBdwnKQlAw2SXg0cAK5Qdg/wU8m2LH5INj30zoG/+CW1Szqthtf7PnBxfrVQhpmS+gbwV1XxnDvC\ncz5BNkU10P9lwKMR8Rmyiw6eV0N8Zkm8hWGlExEh6U3AP+c7q38LPATcSDbd9GOyG9D8TX6p8Eck\nvRL4Xv4Xfy9wFYn3VoiIX+bFaZ2kF+TjLh3U7QZglaT7yP5fbgaWHuU5fyXpu5LuB74O3A8sk9Sf\nx3dNSmxmtfBhtWZmlsRTUmZmlsQFw8zMkrhgmJlZEhcMMzNL4oJhZmZJXDDMzCyJC4aZmSVxwTAz\nsyT/H1T6TzRInf7bAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11a8da5c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# let's also look at the coefficients of the solution\n",
    "plt.scatter(np.arange(len(theta)), theta)\n",
    "plt.ylabel('Magnitude')\n",
    "plt.xlabel('Coefficients')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: more than one split\n",
    "\n",
    "We will now use ADMM to actually do some distributed computations. For this we can split the first subproblem into two ($N=2$) and using the ADMM algorithm specified at the beginning.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def lasso_ADMM_2(A, b, max_iter=100, lam=1.):\n",
    "\n",
    "    # .. split current dataset ..\n",
    "    A_1, A_2 = np.array_split(A, 2)\n",
    "    b_1, b_2 = np.array_split(b, 2)\n",
    "\n",
    "    # .. initialize variables ..\n",
    "    tau = 1.\n",
    "    n_samples, n_features = A.shape\n",
    "    rho = np.zeros(n_features)\n",
    "    u_1 = np.zeros(n_features)\n",
    "    u_2 = np.zeros(n_features)\n",
    "\n",
    "    # .. to keep track of progress ..\n",
    "    obj_fun_history = []\n",
    "\n",
    "    # .. cache inverse matrix ..\n",
    "    AtA_1_inv = np.linalg.pinv(A_1.T.dot(A_1) / n_samples + tau * np.eye(n_features))\n",
    "    AtA_2_inv = np.linalg.pinv(A_2.T.dot(A_2) / n_samples + tau * np.eye(n_features))\n",
    "\n",
    "    for i in range(max_iter):\n",
    "        theta_1 = AtA_1_inv.dot(A_1.T.dot(b_1) / n_samples + tau * (rho - u_1))\n",
    "        theta_2 = AtA_2_inv.dot(A_2.T.dot(b_2) / n_samples + tau * (rho - u_2))\n",
    "        u = (u_1 + u_2) / 2.\n",
    "        theta = (theta_1 + theta_2) / 2.\n",
    "        rho = np.fmax(theta + u - lam /tau, 0) - np.fmax(-lam/tau - theta - u, 0)\n",
    "        u_1 = u_1 + theta_1 - rho\n",
    "        u_2 = u_2 + theta_2 - rho\n",
    "        obj_fun_history.append(obj_fun(theta_1, A, b, lam))\n",
    "\n",
    "    return theta, obj_fun_history\n",
    "\n",
    "theta_1, func_vals = lasso_ADMM_2(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH25JREFUeJzt3Xt4XXWd7/H3J8lO0vSW1GIsbUrB6YCIUKBy8fYER4/A\noHjBI15QvDw9MurgEYejM3N0dM6Zc46oR3hQe8rIoB4fekQZrAyXYbABFUFoLYVyGcq1rYVSaEPT\ne5Lv+WOtbHbTnWQn7dq7yfq8nmc/Xbe91vdbSr75rd/6/ZYiAjMzM4C6WgdgZmaHDhcFMzMrclEw\nM7MiFwUzMytyUTAzsyIXBTMzK3JRMDOzIhcFMzMrclEwM7OihloHMFozZ86MefPmjem727dvZ/Lk\nyQc3oHEgj3nnMWfIZ955zBlGn/eKFSs2R8RhIx037orCvHnzuO+++8b03a6uLjo7Ow9uQONAHvPO\nY86Qz7zzmDOMPm9JT1dynG8fmZlZkYuCmZkVuSiYmVlRZkVBUrOk30u6X9IaSV8rc4wkXSFpraTV\nkk7KKh4zMxtZlh3Nu4G3RkSPpALwG0k3R8TdJcecBcxPP6cC30//NDOzGsispRCJnnS1kH4Gv9Hn\nXOBH6bF3A62SZmUVk5mZDS/TPgVJ9ZJWAZuA2yLinkGHzAbWlayvT7eZmVkNZDpOISL6gAWSWoF/\nlnRcRDw42vNIWgQsAmhvb6erq2vUsazf1s9vntnBS3uWM61Ro/7+eNbT0zOmv7PxLI85Qz7zzmPO\nkF3eVRm8FhFbJS0HzgRKi8IGoKNkfU66bfD3lwBLABYuXBhjGahy8wMbueW3K7n4vSfzmlnTRv39\n8SyPg3vymDPkM+885gzZ5Z3l00eHpS0EJE0C3g48MuiwZcBH06eQTgO6I2JjFvE0FZJUd/f2Z3F6\nM7MJIcuWwizgh5LqSYrPTyPiRkmfBoiIxcBNwNnAWmAH8PGsgmluqAdg196+rC5hZjbuZVYUImI1\ncGKZ7YtLlgP4TFYxlHJLwcxsZLkZ0dzkloKZ2YhyUxSa3VIwMxtRboqCWwpmZiPLTVFoLiRFYbeL\ngpnZkHJTFNzRbGY2stwUBT+SamY2stwUhUK9EG4pmJkNJzdFQRKFercUzMyGk5uiAFCoc0vBzGw4\nuSoKjXVyS8HMbBj5Kgr1sGuvWwpmZkPJVVFIbh+5pWBmNpR8FYV6uaVgZjaMfBUFtxTMzIaVq6KQ\ndDS7pWBmNpRcFYVCvR9JNTMbTr6KQp0nxDMzG06uikJjvccpmJkNJ1dFwSOazcyGl6ui0FjnuY/M\nzIaTq6JQqJdbCmZmw8hXUaiD3v6gt8+FwcysnHwVheQ9O24tmJkNIVdFobFOgPsVzMyGkquiMNBS\n2OWWgplZWfkqCmlLwQPYzMzKy1VRaEyz9fxHZmbl5aoovNzR7JaCmVk5uSoKL3c0u6VgZlZOropC\nIc3WLQUzs/LyVRQGnj5yS8HMrKzMioKkDknLJT0kaY2ki8sc0ympW9Kq9POVrOKBl28fuaVgZlZe\nQ4bn7gUuiYiVkqYCKyTdFhEPDTru1xFxToZxFBU7mt1SMDMrK7OWQkRsjIiV6fI24GFgdlbXq0Sx\no9ktBTOzsqrSpyBpHnAicE+Z3W+QtFrSzZJem2UcbimYmQ1PEZHtBaQpwB3Af4+I6wftmwb0R0SP\npLOByyNifplzLAIWAbS3t5+8dOnSMcXS/VIPF98l3ju/wLte3Timc4xHPT09TJkypdZhVFUec4Z8\n5p3HnGH0eZ9xxhkrImLhSMdl2aeApALwc+AngwsCQES8VLJ8k6TvSZoZEZsHHbcEWAKwcOHC6Ozs\nHFM8y5cvp047mDVnLp2dx4zpHONRV1cXY/07G6/ymDPkM+885gzZ5Z3l00cCfgA8HBHfHuKYV6XH\nIemUNJ4XMoyJ5kK9bx+ZmQ0hy5bCG4ELgAckrUq3/TUwFyAiFgPnARdJ6gV2AudHxvezmgv17mg2\nMxtCZkUhIn4DaIRjrgSuzCqGcpoa6txSMDMbQq5GNMNAS8FFwcysnNwVhaSl4NtHZmbl5K8ouKVg\nZjak/BWFhjq/o9nMbAi5KwrNhXp2u6VgZlZW7oqC+xTMzIaWu6LgloKZ2dDyVxTcp2BmNqTcFYWm\nQp1bCmZmQ8hdUWhuqHdLwcxsCLkrCk2F5PZR1lOGm5mNR7krCs0N9fQH9Pa7KJiZDZa7otBUSFL2\nLSQzs/3lrig0p+/kdGezmdn+8lcUGpKi4JaCmdn+clcUBm4fuaVgZra//BUFtxTMzIaUv6JQ7Gh2\nS8HMbLBRFQVJbZKOzyqYahjoU9jt9zSbme1nxKIgqUvSNEkzgJXAVZK+nX1o2Sj2KbilYGa2n0pa\nCtMj4iXgvcCPIuJU4G3ZhpUdtxTMzIZWSVFokDQL+I/AjRnHk7lm9ymYmQ2pkqLwdeBW4PGIuFfS\nUcBj2YaVnaaCWwpmZkNpGOmAiLgOuK5k/QngfVkGlaXmBrcUzMyGUklH859Kul3Sg+n68ZL+NvvQ\nsjHQUvA4BTOz/VVy++gq4MvAXoCIWA2cn2VQWRpoKXhEs5nZ/iopCi0R8ftB23qzCKYaGurrqK+T\nWwpmZmVUUhQ2S3o1EACSzgM2ZhpVxpob/EpOM7NyRuxoBj4DLAGOkbQBeBL4SKZRZay54FdympmV\nU8nTR08Ab5M0GaiLiG3Zh5WtJrcUzMzKGrEoSPrKoHUAIuLrGcWUObcUzMzKq6RPYXvJpw84C5g3\n0pckdUhaLukhSWskXVzmGEm6QtJaSaslnTTK+MeksaHO4xTMzMqo5PbRt0rXJX2TZITzSHqBSyJi\npaSpwApJt0XEQyXHnAXMTz+nAt9P/8xUc6HeI5rNzMoYy/sUWoA5Ix0UERsjYmW6vA14GJg96LBz\nSSbZi4i4G2hN51nKVFNDnWdJNTMro5I+hQdIH0cF6oHDSOZDqpikecCJwD2Dds0G1pWsr0+37fPI\nq6RFwCKA9vZ2urq6RnP5op6eHrq6utixbRc9e2PM5xlvBvLOkzzmDPnMO485Q3Z5V/JI6jkly73A\ncxFR8eA1SVOAnwOfT6fgHrWIWELyWCwLFy6Mzs7OsZyGrq4uOjs7uXbdfTy1eQednW8Z03nGm4G8\n8ySPOUM+885jzpBd3kMWhfSlOgCDH0GdJomIeHGkk0sqkBSEn0TE9WUO2QB0lKzPSbdlqqnBfQpm\nZuUM11JYQXLbSGX2BXDUcCdW8uzqD4CHI2KoN7UtAz4raSlJB3N3RGQ+Wrq54KePzMzKGbIoRMSR\nB3juNwIXAA9IWpVu+2tgbnr+xcBNwNnAWmAH8PEDvGZFmhrq2eWWgpnZfirpU0BSG8ljo80D2yLi\nzuG+ExG/oXwro/SYIJlGo6qaC376yMysnEqePvoUcDHJ/f5VwGnA74C3ZhtadgZaChFRHKFtZmaV\njVO4GHg98HREnEHyaOnWTKPKWHOhjgjY2xcjH2xmliOVFIVdEbELQFJTRDwCHJ1tWNlqHnj7mvsV\nzMz2UUmfwnpJrcANwG2StgBPZxtWtpoG3r62t7+kl8TMzCqZ++g96eLfSVoOTAduyTSqjPk9zWZm\n5VXS0XwFsDQi7oqIO6oQU+aKLQXfPjIz20clfQorgL+V9Likb0pamHVQWSv2KfixVDOzfYxYFCLi\nhxFxNskTSI8C/0vSY5lHliG3FMzMyhvN1Nl/AhwDHAE8kk041TGlKblrtm1XxfP6mZnlwohFQdI3\n0pbB14EHgIUR8c7MI8tQa0sjAN0799Y4EjOzQ0slj6Q+DpweEZuzDqZa2loKAGzZvqfGkZiZHVoq\neST1/1QjkGqaPiktCjvcUjAzKzWW13GOew31dUxtbvDtIzOzQXJZFADaWhrZssO3j8zMSlU6dXY9\n0F56fEQ8k1VQ1dDWUvDtIzOzQSoZ0fw54KvAc8DAaK8Ajs8wrsxNb2mk2y0FM7N9VNJSuBg4OiJe\nyDqYamprKfDU5u21DsPM7JBSSZ/COqA760CqzX0KZmb7q6Sl8ATQJelfgN0DGyPi25lFVQXTJxXY\ntquX3r5+Gupz299uZraPSorCM+mnMf1MCAMD2Lp37uUVU5pqHI2Z2aGhksFrXwOQNCVd78k6qGpo\nm5zUty07XBTMzAZUMvfRcZL+AKwB1khaIem12YeWrYFRzd073a9gZjagkpvpS4AvRMQREXEEcAlw\nVbZhZa8tnRRvy3aPVTAzG1BJUZgcEcsHViKiC5icWURVUiwKfgLJzKyooqePJP1X4Mfp+kdInkga\n16aXdDSbmVmikpbCJ4DDgOvTz2HptnFtWnMD9XVyS8HMrEQlTx9tAf6yCrFUlSRaJ3n+IzOzUkMW\nBUnfiYjPS/olyVxH+4iId2UaWRW0thTodlEwMysarqUw0IfwzWoEUgutnurCzGwfQxaFiFiRLi6I\niMtL90m6GLgjy8Cqoa2lwIatu2odhpnZIaOSjuaPldl24UGOoyZaPX22mdk+hutT+CDwIeBISctK\ndk0FXhzpxJKuBs4BNkXEcWX2dwK/AJ5MN10fEV+vPPQD545mM7N9DdencBewEZgJfKtk+zZgdQXn\nvga4EvjRMMf8OiLOqeBcmWib3MjOvX3s2ttHc6G+VmGYmR0yhutTeBp4WtKHgT9GxC4ASZOAOcBT\nw504Iu6UNO+gRZqB1nQA29Yde3nVdBcFMzNF7Pe06b4HSPcBb4iIPel6I/DbiHj9iCdPisKNw9w+\nuh5YD2wAvhgRa4Y4zyJgEUB7e/vJS5cuHenSZfX09DBlypTi+u+f7eV7q3bz92+cRMfUiftOhcF5\n50Eec4Z85p3HnGH0eZ9xxhkrImLhiAdGxLAfYFWZbfeP9L30uHnAg0PsmwZMSZfPBh6r5Jwnn3xy\njNXy5cv3Wf/tY8/HEf/lxrhr7eYxn3M8GJx3HuQx54h85p3HnCNGnzdwX1TwM7aSX4+fl1QcqCbp\nXGBzxeVp6GL0UqTvZoiIm4CCpJkHet7RaE0nxdvqJ5DMzIDKJsT7NPATSVcCInln80cP9MKSXgU8\nFxEh6RSSx2NfONDzjkbb5LRPwZPimZkBlc199Dhw2mjfvCbpWqATmClpPfBVoJCeYzFwHnCRpF5g\nJ3B+2sSpmtZJnj7bzKzUiEVBUhPwPpL+gQZJAMQIYwoi4oMj7L+S5JHVmpnUWE9TQx1bPVbBzAyo\n7PbRL4BuYAWwO9twqq+tpdF9CmZmqUqKwpyIODPzSGqktcWjms3MBlTy9NFdkl6XeSQ10tpScEvB\nzCxVSUvhTcCFkp4kuX0kICLi+Ewjq5K2lkbWbqqo79zMbMKrpCiclXkUNeTbR2ZmL6ukKFT1MdFq\na007miOCgSerzMzyqpKi8C8khUFAM3Ak8Cjw2gzjqpq2lgK9/cH2PX1Maarkr8PMbOKqZPDaPp3M\nkk4C/iKziKpsYKqLLdv3uCiYWe6NemrQiFgJnJpBLDXROunl6bPNzPKukhHNXyhZrQNOAv6YWURV\n1jY5nRRvpx9LNTOr5H7J1JLlXpI+hp9nE071zZzSBMBzL024wdpmZqM23DuafxwRFwBbI+LyKsZU\nVYe3NiPBuhd31DoUM7OaG65P4WRJhwOfkNQmaUbpp1oBZq2poZ5Z05pZt8VFwcxsuNtHi4HbgaNI\nJsMrfYg/0u0TwpwZLW4pmJkxTEshIq6IiNcAV0fEURFxZMlnwhQEgLkzWlj34s5ah2FmVnMjPpIa\nERdVI5Ba6mhr4dmXdrFrb1+tQzEzq6lRj1OYiDpmTAJgw1a3Fsws31wUSG4fgZ9AMjNzUQA6XBTM\nzAAXBQAOm9JEU0Md67b49pGZ5ZuLAlBXJ+a0TeKZF9xSMLN8c1FIdcxo8QA2M8s9F4XU3BktPOM+\nBTPLOReFVEdbC9t29dLtKbTNLMdcFFLFJ5B8C8nMcsxFITUwgM23kMwsz1wUUh6rYGbmolA0rblA\na0vBLQUzyzUXhRIdbS0ewGZmueaiUGKu36tgZjmXWVGQdLWkTZIeHGK/JF0haa2k1ZJOyiqWSs2Z\nMYkNW3bS3x+1DsXMrCaybClcA5w5zP6zgPnpZxHw/QxjqUhHWwt7+vp5btuuWodiZlYTmRWFiLgT\neHGYQ84FfhSJu4FWSbOyiqcSA1Noew4kM8urWvYpzAbWlayvT7fVzEBReNpFwcxyqqHWAVRC0iKS\nW0y0t7fT1dU1pvP09PQM+93+CJrr4ZZ7H+KV2x8f0zUORSPlPRHlMWfIZ955zBmyy7uWRWED0FGy\nPifdtp+IWAIsAVi4cGF0dnaO6YJdXV2M9N0TH7ub53f30tn5pjFd41BUSd4TTR5zhnzmncecIbu8\na3n7aBnw0fQppNOA7ojYWMN4ADhxbisPb3yJXXv7ah2KmVnVZdZSkHQt0AnMlLQe+CpQAIiIxcBN\nwNnAWmAH8PGsYhmNBR2t9PYHD27oZuG8GbUOx8ysqjIrChHxwRH2B/CZrK4/VgvmtgKwat1WFwUz\nyx2PaB7klVObmd06iT+s21rrUMzMqs5FoYwFc1tZ9YyLgpnlj4tCGSd2tLJh6042eWSzmeWMi0IZ\nCzrSfgW3FswsZ1wUyjhu9nQa6sQq9yuYWc64KJTRXKjnNbOmuSiYWe64KAxhQUcr96/bSp+n0Taz\nHHFRGMKCjla27+lj7aaeWodiZlY1LgpDODEdxLbymS01jsTMrHpcFIZw5MzJzJrezL899FytQzEz\nqxoXhSFI4s9fN4s7H3ue7h17ax2OmVlVuCgM45wTDmdvX3DrQ8/WOhQzs6pwURjGCXOm0zFjEjeu\nrvmM3mZmVeGiMAxJnHP84fx27WZe3L6n1uGYmWXORWEE5xw/i77+4OYH3Vows4nPRWEEx86axlEz\nJ3Pj/S4KZjbxuSiMQBLnnHA49zz5gmdNNbMJz0WhAu88fhb9gVsLZjbhuShUYH77VE6a28qSO59g\n196+WodjZpYZF4UKXXrmMTz70i6uueupWodiZpYZF4UKnXbUK+g8+jC+t3ytRzib2YTlojAKl77j\nGLbt7uX7dzxe61DMzDLhojAKxx4+jXcvmM0//fZJnu32k0hmNvG4KIzSF97+p/RH8HfL1tDvF/CY\n2QTjojBKHTNauPQdx3DLmme5/PbHah2OmdlB1VDrAMajT735SP79uW1cfvtjzG+fwjnHH17rkMzM\nDgq3FMZAEv/tPcex8Ig2Lvnp/dy/bmutQzIzOyhcFMaoqaGexReczGFTm/jQVXdz8wMe7Wxm45+L\nwgGYOaWJ6z59OvPbp3LRT1byP29+hD53PpvZOOaicIBmTZ/E//tPp/GhU+ey+I7Hef/iu7jvqRdr\nHZaZ2Zi4KBwETQ31/MN7Xse33n8C67bs5LzFv+OT19zLA+u7iXDLwczGj0yLgqQzJT0qaa2kL5XZ\n3ympW9Kq9POVLOPJ2vtOnsMdf9XJpWcezb1Pvcg7r/wN7/jOnXx3+VqefmG7C4SZHfIyeyRVUj3w\nXeDtwHrgXknLIuKhQYf+OiLOySqOamtpbOAvOv+ED596BMvu/yO/+MMGLrv1US679VFmTW/mlCNn\ncPIRbcx/5VTmt0/hFZMbkVTrsM3MgGzHKZwCrI2IJwAkLQXOBQYXhQlp+qQCF5x2BBecdgTrXtxB\n16ObuOfJF7nr8Rf4xao/7nPcrOnNzJrezKumNzNjciNtLY20tjQytbmBKU0NTG5qYFKhnkmFepoL\ndTQ21FGoT/5sqJOLipkdNFkWhdnAupL19cCpZY57g6TVwAbgixGxJsOYaqJjRgsXnD6PC06fR0Sw\nsXsXazf1sHZTD09u3s7G7p1s7N7FAxu62bJj76ifYKqvE/V1oqFO1EvUpet1SsZU7N2zh+a7bi+u\nA0jpBzFQU0TJ/tILqOziPqpRmEZzhe07djB55R2ZxXKoymPeecr5A6/v4FNvPirTa9R6RPNKYG5E\n9Eg6G7gBmD/4IEmLgEUA7e3tdHV1jeliPT09Y/5uFo4CjmoFWge2FOiPBnb2wva9wc7eYFcv7OwN\n9vTDnr5gTx/09UNvwN7+oK8f+gP6Avojin8GQUSyb8/efhoaeglgoFsjgCCSBYp/JMslK5WUp2r0\nlIy2O2ZKUz/12plNMIewPOadp5w3rXuCrq5ngOx+nmVZFDYAHSXrc9JtRRHxUsnyTZK+J2lmRGwe\ndNwSYAnAwoULo7Ozc0wBdXV1Mdbvjmd5zDuPOUM+885jzpBd3lk+fXQvMF/SkZIagfOBZaUHSHqV\n0vsOkk5J43khw5jMzGwYmbUUIqJX0meBW4F64OqIWCPp0+n+xcB5wEWSeoGdwPnh5zbNzGom0z6F\niLgJuGnQtsUly1cCV2YZg5mZVc4jms3MrMhFwczMilwUzMysyEXBzMyKXBTMzKxI4+0JUEnPA0+P\n8eszgc0jHjXx5DHvPOYM+cw7jznD6PM+IiIOG+mgcVcUDoSk+yJiYa3jqLY85p3HnCGfeecxZ8gu\nb98+MjOzIhcFMzMryltRWFLrAGokj3nnMWfIZ955zBkyyjtXfQpmZja8vLUUzMxsGLkpCpLOlPSo\npLWSvlTreLIgqUPSckkPSVoj6eJ0+wxJt0l6LP2zrdaxHmyS6iX9QdKN6Xoecm6V9DNJj0h6WNLp\nOcn7P6f/vh+UdK2k5omWt6SrJW2S9GDJtiFzlPTl9Gfbo5LecSDXzkVRkFQPfBc4CzgW+KCkY2sb\nVSZ6gUsi4ljgNOAzaZ5fAm6PiPnA7en6RHMx8HDJeh5yvhy4JSKOAU4gyX9C5y1pNvCXwMKIOI5k\nWv7zmXh5XwOcOWhb2RzT/8fPB16bfud76c+8MclFUQBOAdZGxBMRsQdYCpxb45gOuojYGBEr0+Vt\nJD8kZpPk+sP0sB8C765NhNmQNAf4c+AfSzZP9JynA28BfgAQEXsiYisTPO9UAzBJUgPQAvyRCZZ3\nRNwJvDho81A5ngssjYjdEfEksJbkZ96Y5KUozAbWlayvT7dNWJLmAScC9wDtEbEx3fUs0F6jsLLy\nHeBSoL9k20TP+UjgeeCf0ttm/yhpMhM874jYAHwTeAbYCHRHxL8ywfNODZXjQf35lpeikCuSpgA/\nBz5f+h5sgPTNdhPmkTNJ5wCbImLFUMdMtJxTDcBJwPcj4kRgO4NumUzEvNP76OeSFMXDgcmSPlJ6\nzETMe7Asc8xLUdgAdJSsz0m3TTiSCiQF4ScRcX26+TlJs9L9s4BNtYovA28E3iXpKZLbgm+V9H+Z\n2DlD8tvg+oi4J13/GUmRmOh5vw14MiKej4i9wPXAG5j4ecPQOR7Un295KQr3AvMlHSmpkaRTZlmN\nYzroJInkHvPDEfHtkl3LgI+lyx8DflHt2LISEV+OiDkRMY/kv+uvIuIjTOCcASLiWWCdpKPTTX8G\nPMQEz5vkttFpklrSf+9/RtJ3NtHzhqFzXAacL6lJ0pHAfOD3Y75KROTiA5wN/DvwOPA3tY4noxzf\nRNKkXA2sSj9nA68geVrhMeDfgBm1jjWj/DuBG9PlCZ8zsAC4L/3vfQPQlpO8vwY8AjwI/Bhommh5\nA9eS9JnsJWkVfnK4HIG/SX+2PQqcdSDX9ohmMzMrysvtIzMzq4CLgpmZFbkomJlZkYuCmZkVuSiY\nmVmRi4LlgqQuSQf9fbaSLktn7Lxs0PZ3DczGK+ndB3MCRkkLJJ1d7lpmB8qPpFouSOoCvhgR9x3k\n83aTPC/eN8wx15CMn/jZKM7bEBG9Q+y7kGSW0M+OMlyzEbmlYIcMSfPS9wJclf72/a+SJqX7ir/p\nS5qZTmuBpAsl3ZDOL/+UpM9K+kI6SdzdkmaUXOICSavSefhPSb8/OZ27/vfpd/abPVeJy9LvPSDp\nA+n2ZcAUYMXAtpLvXCjpSklvAN4FXJZe+9Xp5xZJKyT9WtIx6XeukbRY0j3ANySdIul3aVx3STo6\nHZH/deAD6fk+MHCtkr/DX0laLel2SXNLzn1Fep4nJJ2Xbp8l6c6Sv5c3H6T/nDZe1Xrknj/+DHyA\neSTvhFiQrv8U+Ei63EXy2zHATOCpdPlCkqmCpwKHAd3Ap9N9/5tkUsCB71+VLr8FeDBd/oeSa7SS\njHqfPCiu9wG3kczd304y1cKsdF/PELlcCFyZLl8DnFey73Zgfrp8KsnUHAPH3QjUp+vTgIZ0+W3A\nzwefu8y1fgl8LF3+BHBDybmvI/lF8FiSqeQBLiEd4Z/mN7XW/w78qe2nYdiKYVZ9T0bEqnR5BUmh\nGMnySN4fsS29nfPLdPsDwPElx10LyVz1kqZJagX+A8mEel9Mj2kG5rLvC3veBFwbyS2i5yTdAbye\nMcyflc5g+wbgumTqHiCZpmHAdfHyrajpwA8lzSeZvqRQwSVOB96bLv8Y+EbJvhsioh94SNLAtMv3\nAlenEyneUPJ3bznlomCHmt0ly33ApHS5l5dvdzYP853+kvV+9v03PrgDLQAB74uIR8ca8CjVAVsj\nYsEQ+7eXLP89ScF7j5L3Y3Qd4LVL/54ExQL5FpKXFF0j6dsR8aMDvI6NY+5TsPHiKeDkdPm8MZ5j\noC/gTSQvZ+kGbgU+l864iaQTy3zv1yT38OslHUZy+2k0s1BuI7m9RSTvt3hS0vvT60nSCUN8bzov\nT4F8YbnzlXEXyWyxAB9OYx+SpCOA5yLiKpI31500bCY24bko2HjxTeAiSX8g6VMYi13p9xeTzDoJ\nyW/jBWC1pDXp+mD/TDIT6f3Ar4BLI5m6ulJLgb9KO4xfTfLD+pOS7gfWMPSrYb8B/I805tIWz3Lg\n2IGO5kHf+RzwcUmrgQtI3l09nE7g/vQaHyB577PlmB9JNTOzIrcUzMysyEXBzMyKXBTMzKzIRcHM\nzIpcFMzMrMhFwczMilwUzMysyEXBzMyK/j/+N/QRFtWIWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116bdc898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# lets plot the objective values of the function\n",
    "# to make sure it has converged\n",
    "plt.plot(func_vals)\n",
    "plt.ylabel('function values')\n",
    "plt.xlabel('numbe of iterations')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color: red; font-weight: bold\">Exercise (10 min): Extend this to ar arbitrary number of splits.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now use the multiprocessing library in the Python standard library to perform the theta_1 and theta_2 steps in parallel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import concurrent.futures\n",
    "\n",
    "def _compute_inverse(A_1, tau, n_features):\n",
    "    return np.linalg.pinv(A_1.T.dot(A_1) / n_samples + tau * np.eye(n_features))\n",
    "\n",
    "def _compute_theta(AtA_inv, A, b, tau, rho, u):\n",
    "    return AtA_inv.dot(A.T.dot(b) / n_samples + tau * (rho - u))\n",
    "\n",
    "def lasso_ADMM_2_MP(A, b, max_iter=100, lam=1.):\n",
    "\n",
    "    # .. split current dataset ..\n",
    "    A_1, A_2 = np.array_split(A, 2)\n",
    "    b_1, b_2 = np.array_split(b, 2)\n",
    "\n",
    "    # .. initialize variables ..\n",
    "    tau = 1.\n",
    "    n_samples, n_features = A.shape\n",
    "    rho = np.zeros(n_features)\n",
    "    u_1 = np.zeros(n_features)\n",
    "    u_2 = np.zeros(n_features)\n",
    "\n",
    "    # .. to keep track of progress ..\n",
    "    obj_fun_history = []\n",
    "\n",
    "    \n",
    "    futures = []\n",
    "    with concurrent.futures.ProcessPoolExecutor() as executor:\n",
    "        futures.append(executor.submit(_compute_inverse, A_1, tau, n_features))\n",
    "        futures.append(executor.submit(_compute_inverse, A_2, tau, n_features))    \n",
    "    concurrent.futures.wait(futures)\n",
    "\n",
    "    AtA_1_inv = futures[0].result()\n",
    "    AtA_2_inv = futures[1].result()\n",
    "\n",
    "\n",
    "    for i in range(max_iter):\n",
    "        futures = []\n",
    "        with concurrent.futures.ProcessPoolExecutor() as executor:\n",
    "            futures.append(executor.submit(_compute_theta, AtA_1_inv, A_1, b_1, tau, rho, u_1))\n",
    "            futures.append(executor.submit(_compute_theta, AtA_2_inv, A_2, b_2, tau, rho, u_2))\n",
    "        concurrent.futures.wait(futures)\n",
    "        theta_1 = futures[0].result()\n",
    "        theta_2 = futures[1].result()\n",
    "        theta = (theta_1 + theta_2) / 2.\n",
    "        u = (u_1 + u_2) / 2.\n",
    "        rho = np.fmax(theta + u - lam /tau, 0) - np.fmax(-lam/tau - theta - u, 0)\n",
    "        u_1 = u_1 + theta_1 - rho\n",
    "        u_2 = u_2 + theta_2 - rho\n",
    "        obj_fun_history.append(obj_fun(theta_1, A, b, lam))\n",
    "\n",
    "    return theta, obj_fun_history\n",
    "\n",
    "theta_1, func_vals = lasso_ADMM_2_MP(A, b)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH25JREFUeJzt3Xt4XXWd7/H3J8lO0vSW1GIsbUrB6YCIUKBy8fYER4/A\noHjBI15QvDw9MurgEYejM3N0dM6Zc46oR3hQe8rIoB4fekQZrAyXYbABFUFoLYVyGcq1rYVSaEPT\ne5Lv+WOtbHbTnWQn7dq7yfq8nmc/Xbe91vdbSr75rd/6/ZYiAjMzM4C6WgdgZmaHDhcFMzMrclEw\nM7MiFwUzMytyUTAzsyIXBTMzK3JRMDOzIhcFMzMrclEwM7OihloHMFozZ86MefPmjem727dvZ/Lk\nyQc3oHEgj3nnMWfIZ955zBlGn/eKFSs2R8RhIx037orCvHnzuO+++8b03a6uLjo7Ow9uQONAHvPO\nY86Qz7zzmDOMPm9JT1dynG8fmZlZkYuCmZkVuSiYmVlRZkVBUrOk30u6X9IaSV8rc4wkXSFpraTV\nkk7KKh4zMxtZlh3Nu4G3RkSPpALwG0k3R8TdJcecBcxPP6cC30//NDOzGsispRCJnnS1kH4Gv9Hn\nXOBH6bF3A62SZmUVk5mZDS/TPgVJ9ZJWAZuA2yLinkGHzAbWlayvT7eZmVkNZDpOISL6gAWSWoF/\nlnRcRDw42vNIWgQsAmhvb6erq2vUsazf1s9vntnBS3uWM61Ro/7+eNbT0zOmv7PxLI85Qz7zzmPO\nkF3eVRm8FhFbJS0HzgRKi8IGoKNkfU66bfD3lwBLABYuXBhjGahy8wMbueW3K7n4vSfzmlnTRv39\n8SyPg3vymDPkM+885gzZ5Z3l00eHpS0EJE0C3g48MuiwZcBH06eQTgO6I2JjFvE0FZJUd/f2Z3F6\nM7MJIcuWwizgh5LqSYrPTyPiRkmfBoiIxcBNwNnAWmAH8PGsgmluqAdg196+rC5hZjbuZVYUImI1\ncGKZ7YtLlgP4TFYxlHJLwcxsZLkZ0dzkloKZ2YhyUxSa3VIwMxtRboqCWwpmZiPLTVFoLiRFYbeL\ngpnZkHJTFNzRbGY2stwUBT+SamY2stwUhUK9EG4pmJkNJzdFQRKFercUzMyGk5uiAFCoc0vBzGw4\nuSoKjXVyS8HMbBj5Kgr1sGuvWwpmZkPJVVFIbh+5pWBmNpR8FYV6uaVgZjaMfBUFtxTMzIaVq6KQ\ndDS7pWBmNpRcFYVCvR9JNTMbTr6KQp0nxDMzG06uikJjvccpmJkNJ1dFwSOazcyGl6ui0FjnuY/M\nzIaTq6JQqJdbCmZmw8hXUaiD3v6gt8+FwcysnHwVheQ9O24tmJkNIVdFobFOgPsVzMyGkquiMNBS\n2OWWgplZWfkqCmlLwQPYzMzKy1VRaEyz9fxHZmbl5aoovNzR7JaCmVk5uSoKL3c0u6VgZlZOropC\nIc3WLQUzs/LyVRQGnj5yS8HMrKzMioKkDknLJT0kaY2ki8sc0ympW9Kq9POVrOKBl28fuaVgZlZe\nQ4bn7gUuiYiVkqYCKyTdFhEPDTru1xFxToZxFBU7mt1SMDMrK7OWQkRsjIiV6fI24GFgdlbXq0Sx\no9ktBTOzsqrSpyBpHnAicE+Z3W+QtFrSzZJem2UcbimYmQ1PEZHtBaQpwB3Af4+I6wftmwb0R0SP\npLOByyNifplzLAIWAbS3t5+8dOnSMcXS/VIPF98l3ju/wLte3Timc4xHPT09TJkypdZhVFUec4Z8\n5p3HnGH0eZ9xxhkrImLhSMdl2aeApALwc+AngwsCQES8VLJ8k6TvSZoZEZsHHbcEWAKwcOHC6Ozs\nHFM8y5cvp047mDVnLp2dx4zpHONRV1cXY/07G6/ymDPkM+885gzZ5Z3l00cCfgA8HBHfHuKYV6XH\nIemUNJ4XMoyJ5kK9bx+ZmQ0hy5bCG4ELgAckrUq3/TUwFyAiFgPnARdJ6gV2AudHxvezmgv17mg2\nMxtCZkUhIn4DaIRjrgSuzCqGcpoa6txSMDMbQq5GNMNAS8FFwcysnNwVhaSl4NtHZmbl5K8ouKVg\nZjak/BWFhjq/o9nMbAi5KwrNhXp2u6VgZlZW7oqC+xTMzIaWu6LgloKZ2dDyVxTcp2BmNqTcFYWm\nQp1bCmZmQ8hdUWhuqHdLwcxsCLkrCk2F5PZR1lOGm5mNR7krCs0N9fQH9Pa7KJiZDZa7otBUSFL2\nLSQzs/3lrig0p+/kdGezmdn+8lcUGpKi4JaCmdn+clcUBm4fuaVgZra//BUFtxTMzIaUv6JQ7Gh2\nS8HMbLBRFQVJbZKOzyqYahjoU9jt9zSbme1nxKIgqUvSNEkzgJXAVZK+nX1o2Sj2KbilYGa2n0pa\nCtMj4iXgvcCPIuJU4G3ZhpUdtxTMzIZWSVFokDQL+I/AjRnHk7lm9ymYmQ2pkqLwdeBW4PGIuFfS\nUcBj2YaVnaaCWwpmZkNpGOmAiLgOuK5k/QngfVkGlaXmBrcUzMyGUklH859Kul3Sg+n68ZL+NvvQ\nsjHQUvA4BTOz/VVy++gq4MvAXoCIWA2cn2VQWRpoKXhEs5nZ/iopCi0R8ftB23qzCKYaGurrqK+T\nWwpmZmVUUhQ2S3o1EACSzgM2ZhpVxpob/EpOM7NyRuxoBj4DLAGOkbQBeBL4SKZRZay54FdympmV\nU8nTR08Ab5M0GaiLiG3Zh5WtJrcUzMzKGrEoSPrKoHUAIuLrGcWUObcUzMzKq6RPYXvJpw84C5g3\n0pckdUhaLukhSWskXVzmGEm6QtJaSaslnTTK+MeksaHO4xTMzMqo5PbRt0rXJX2TZITzSHqBSyJi\npaSpwApJt0XEQyXHnAXMTz+nAt9P/8xUc6HeI5rNzMoYy/sUWoA5Ix0UERsjYmW6vA14GJg96LBz\nSSbZi4i4G2hN51nKVFNDnWdJNTMro5I+hQdIH0cF6oHDSOZDqpikecCJwD2Dds0G1pWsr0+37fPI\nq6RFwCKA9vZ2urq6RnP5op6eHrq6utixbRc9e2PM5xlvBvLOkzzmDPnMO485Q3Z5V/JI6jkly73A\ncxFR8eA1SVOAnwOfT6fgHrWIWELyWCwLFy6Mzs7OsZyGrq4uOjs7uXbdfTy1eQednW8Z03nGm4G8\n8ySPOUM+885jzpBd3kMWhfSlOgCDH0GdJomIeHGkk0sqkBSEn0TE9WUO2QB0lKzPSbdlqqnBfQpm\nZuUM11JYQXLbSGX2BXDUcCdW8uzqD4CHI2KoN7UtAz4raSlJB3N3RGQ+Wrq54KePzMzKGbIoRMSR\nB3juNwIXAA9IWpVu+2tgbnr+xcBNwNnAWmAH8PEDvGZFmhrq2eWWgpnZfirpU0BSG8ljo80D2yLi\nzuG+ExG/oXwro/SYIJlGo6qaC376yMysnEqePvoUcDHJ/f5VwGnA74C3ZhtadgZaChFRHKFtZmaV\njVO4GHg98HREnEHyaOnWTKPKWHOhjgjY2xcjH2xmliOVFIVdEbELQFJTRDwCHJ1tWNlqHnj7mvsV\nzMz2UUmfwnpJrcANwG2StgBPZxtWtpoG3r62t7+kl8TMzCqZ++g96eLfSVoOTAduyTSqjPk9zWZm\n5VXS0XwFsDQi7oqIO6oQU+aKLQXfPjIz20clfQorgL+V9Likb0pamHVQWSv2KfixVDOzfYxYFCLi\nhxFxNskTSI8C/0vSY5lHliG3FMzMyhvN1Nl/AhwDHAE8kk041TGlKblrtm1XxfP6mZnlwohFQdI3\n0pbB14EHgIUR8c7MI8tQa0sjAN0799Y4EjOzQ0slj6Q+DpweEZuzDqZa2loKAGzZvqfGkZiZHVoq\neST1/1QjkGqaPiktCjvcUjAzKzWW13GOew31dUxtbvDtIzOzQXJZFADaWhrZssO3j8zMSlU6dXY9\n0F56fEQ8k1VQ1dDWUvDtIzOzQSoZ0fw54KvAc8DAaK8Ajs8wrsxNb2mk2y0FM7N9VNJSuBg4OiJe\nyDqYamprKfDU5u21DsPM7JBSSZ/COqA760CqzX0KZmb7q6Sl8ATQJelfgN0DGyPi25lFVQXTJxXY\ntquX3r5+Gupz299uZraPSorCM+mnMf1MCAMD2Lp37uUVU5pqHI2Z2aGhksFrXwOQNCVd78k6qGpo\nm5zUty07XBTMzAZUMvfRcZL+AKwB1khaIem12YeWrYFRzd073a9gZjagkpvpS4AvRMQREXEEcAlw\nVbZhZa8tnRRvy3aPVTAzG1BJUZgcEcsHViKiC5icWURVUiwKfgLJzKyooqePJP1X4Mfp+kdInkga\n16aXdDSbmVmikpbCJ4DDgOvTz2HptnFtWnMD9XVyS8HMrEQlTx9tAf6yCrFUlSRaJ3n+IzOzUkMW\nBUnfiYjPS/olyVxH+4iId2UaWRW0thTodlEwMysarqUw0IfwzWoEUgutnurCzGwfQxaFiFiRLi6I\niMtL90m6GLgjy8Cqoa2lwIatu2odhpnZIaOSjuaPldl24UGOoyZaPX22mdk+hutT+CDwIeBISctK\ndk0FXhzpxJKuBs4BNkXEcWX2dwK/AJ5MN10fEV+vPPQD545mM7N9DdencBewEZgJfKtk+zZgdQXn\nvga4EvjRMMf8OiLOqeBcmWib3MjOvX3s2ttHc6G+VmGYmR0yhutTeBp4WtKHgT9GxC4ASZOAOcBT\nw504Iu6UNO+gRZqB1nQA29Yde3nVdBcFMzNF7Pe06b4HSPcBb4iIPel6I/DbiHj9iCdPisKNw9w+\nuh5YD2wAvhgRa4Y4zyJgEUB7e/vJS5cuHenSZfX09DBlypTi+u+f7eV7q3bz92+cRMfUiftOhcF5\n50Eec4Z85p3HnGH0eZ9xxhkrImLhiAdGxLAfYFWZbfeP9L30uHnAg0PsmwZMSZfPBh6r5Jwnn3xy\njNXy5cv3Wf/tY8/HEf/lxrhr7eYxn3M8GJx3HuQx54h85p3HnCNGnzdwX1TwM7aSX4+fl1QcqCbp\nXGBzxeVp6GL0UqTvZoiIm4CCpJkHet7RaE0nxdvqJ5DMzIDKJsT7NPATSVcCInln80cP9MKSXgU8\nFxEh6RSSx2NfONDzjkbb5LRPwZPimZkBlc199Dhw2mjfvCbpWqATmClpPfBVoJCeYzFwHnCRpF5g\nJ3B+2sSpmtZJnj7bzKzUiEVBUhPwPpL+gQZJAMQIYwoi4oMj7L+S5JHVmpnUWE9TQx1bPVbBzAyo\n7PbRL4BuYAWwO9twqq+tpdF9CmZmqUqKwpyIODPzSGqktcWjms3MBlTy9NFdkl6XeSQ10tpScEvB\nzCxVSUvhTcCFkp4kuX0kICLi+Ewjq5K2lkbWbqqo79zMbMKrpCiclXkUNeTbR2ZmL6ukKFT1MdFq\na007miOCgSerzMzyqpKi8C8khUFAM3Ak8Cjw2gzjqpq2lgK9/cH2PX1Maarkr8PMbOKqZPDaPp3M\nkk4C/iKziKpsYKqLLdv3uCiYWe6NemrQiFgJnJpBLDXROunl6bPNzPKukhHNXyhZrQNOAv6YWURV\n1jY5nRRvpx9LNTOr5H7J1JLlXpI+hp9nE071zZzSBMBzL024wdpmZqM23DuafxwRFwBbI+LyKsZU\nVYe3NiPBuhd31DoUM7OaG65P4WRJhwOfkNQmaUbpp1oBZq2poZ5Z05pZt8VFwcxsuNtHi4HbgaNI\nJsMrfYg/0u0TwpwZLW4pmJkxTEshIq6IiNcAV0fEURFxZMlnwhQEgLkzWlj34s5ah2FmVnMjPpIa\nERdVI5Ba6mhr4dmXdrFrb1+tQzEzq6lRj1OYiDpmTAJgw1a3Fsws31wUSG4fgZ9AMjNzUQA6XBTM\nzAAXBQAOm9JEU0Md67b49pGZ5ZuLAlBXJ+a0TeKZF9xSMLN8c1FIdcxo8QA2M8s9F4XU3BktPOM+\nBTPLOReFVEdbC9t29dLtKbTNLMdcFFLFJ5B8C8nMcsxFITUwgM23kMwsz1wUUh6rYGbmolA0rblA\na0vBLQUzyzUXhRIdbS0ewGZmueaiUGKu36tgZjmXWVGQdLWkTZIeHGK/JF0haa2k1ZJOyiqWSs2Z\nMYkNW3bS3x+1DsXMrCaybClcA5w5zP6zgPnpZxHw/QxjqUhHWwt7+vp5btuuWodiZlYTmRWFiLgT\neHGYQ84FfhSJu4FWSbOyiqcSA1Noew4kM8urWvYpzAbWlayvT7fVzEBReNpFwcxyqqHWAVRC0iKS\nW0y0t7fT1dU1pvP09PQM+93+CJrr4ZZ7H+KV2x8f0zUORSPlPRHlMWfIZ955zBmyy7uWRWED0FGy\nPifdtp+IWAIsAVi4cGF0dnaO6YJdXV2M9N0TH7ub53f30tn5pjFd41BUSd4TTR5zhnzmncecIbu8\na3n7aBnw0fQppNOA7ojYWMN4ADhxbisPb3yJXXv7ah2KmVnVZdZSkHQt0AnMlLQe+CpQAIiIxcBN\nwNnAWmAH8PGsYhmNBR2t9PYHD27oZuG8GbUOx8ysqjIrChHxwRH2B/CZrK4/VgvmtgKwat1WFwUz\nyx2PaB7klVObmd06iT+s21rrUMzMqs5FoYwFc1tZ9YyLgpnlj4tCGSd2tLJh6042eWSzmeWMi0IZ\nCzrSfgW3FswsZ1wUyjhu9nQa6sQq9yuYWc64KJTRXKjnNbOmuSiYWe64KAxhQUcr96/bSp+n0Taz\nHHFRGMKCjla27+lj7aaeWodiZlY1LgpDODEdxLbymS01jsTMrHpcFIZw5MzJzJrezL899FytQzEz\nqxoXhSFI4s9fN4s7H3ue7h17ax2OmVlVuCgM45wTDmdvX3DrQ8/WOhQzs6pwURjGCXOm0zFjEjeu\nrvmM3mZmVeGiMAxJnHP84fx27WZe3L6n1uGYmWXORWEE5xw/i77+4OYH3Vows4nPRWEEx86axlEz\nJ3Pj/S4KZjbxuSiMQBLnnHA49zz5gmdNNbMJz0WhAu88fhb9gVsLZjbhuShUYH77VE6a28qSO59g\n196+WodjZpYZF4UKXXrmMTz70i6uueupWodiZpYZF4UKnXbUK+g8+jC+t3ytRzib2YTlojAKl77j\nGLbt7uX7dzxe61DMzDLhojAKxx4+jXcvmM0//fZJnu32k0hmNvG4KIzSF97+p/RH8HfL1tDvF/CY\n2QTjojBKHTNauPQdx3DLmme5/PbHah2OmdlB1VDrAMajT735SP79uW1cfvtjzG+fwjnHH17rkMzM\nDgq3FMZAEv/tPcex8Ig2Lvnp/dy/bmutQzIzOyhcFMaoqaGexReczGFTm/jQVXdz8wMe7Wxm45+L\nwgGYOaWJ6z59OvPbp3LRT1byP29+hD53PpvZOOaicIBmTZ/E//tPp/GhU+ey+I7Hef/iu7jvqRdr\nHZaZ2Zi4KBwETQ31/MN7Xse33n8C67bs5LzFv+OT19zLA+u7iXDLwczGj0yLgqQzJT0qaa2kL5XZ\n3ympW9Kq9POVLOPJ2vtOnsMdf9XJpWcezb1Pvcg7r/wN7/jOnXx3+VqefmG7C4SZHfIyeyRVUj3w\nXeDtwHrgXknLIuKhQYf+OiLOySqOamtpbOAvOv+ED596BMvu/yO/+MMGLrv1US679VFmTW/mlCNn\ncPIRbcx/5VTmt0/hFZMbkVTrsM3MgGzHKZwCrI2IJwAkLQXOBQYXhQlp+qQCF5x2BBecdgTrXtxB\n16ObuOfJF7nr8Rf4xao/7nPcrOnNzJrezKumNzNjciNtLY20tjQytbmBKU0NTG5qYFKhnkmFepoL\ndTQ21FGoT/5sqJOLipkdNFkWhdnAupL19cCpZY57g6TVwAbgixGxJsOYaqJjRgsXnD6PC06fR0Sw\nsXsXazf1sHZTD09u3s7G7p1s7N7FAxu62bJj76ifYKqvE/V1oqFO1EvUpet1SsZU7N2zh+a7bi+u\nA0jpBzFQU0TJ/tILqOziPqpRmEZzhe07djB55R2ZxXKoymPeecr5A6/v4FNvPirTa9R6RPNKYG5E\n9Eg6G7gBmD/4IEmLgEUA7e3tdHV1jeliPT09Y/5uFo4CjmoFWge2FOiPBnb2wva9wc7eYFcv7OwN\n9vTDnr5gTx/09UNvwN7+oK8f+gP6Avojin8GQUSyb8/efhoaeglgoFsjgCCSBYp/JMslK5WUp2r0\nlIy2O2ZKUz/12plNMIewPOadp5w3rXuCrq5ngOx+nmVZFDYAHSXrc9JtRRHxUsnyTZK+J2lmRGwe\ndNwSYAnAwoULo7Ozc0wBdXV1Mdbvjmd5zDuPOUM+885jzpBd3lk+fXQvMF/SkZIagfOBZaUHSHqV\n0vsOkk5J43khw5jMzGwYmbUUIqJX0meBW4F64OqIWCPp0+n+xcB5wEWSeoGdwPnh5zbNzGom0z6F\niLgJuGnQtsUly1cCV2YZg5mZVc4jms3MrMhFwczMilwUzMysyEXBzMyKXBTMzKxI4+0JUEnPA0+P\n8eszgc0jHjXx5DHvPOYM+cw7jznD6PM+IiIOG+mgcVcUDoSk+yJiYa3jqLY85p3HnCGfeecxZ8gu\nb98+MjOzIhcFMzMryltRWFLrAGokj3nnMWfIZ955zBkyyjtXfQpmZja8vLUUzMxsGLkpCpLOlPSo\npLWSvlTreLIgqUPSckkPSVoj6eJ0+wxJt0l6LP2zrdaxHmyS6iX9QdKN6Xoecm6V9DNJj0h6WNLp\nOcn7P6f/vh+UdK2k5omWt6SrJW2S9GDJtiFzlPTl9Gfbo5LecSDXzkVRkFQPfBc4CzgW+KCkY2sb\nVSZ6gUsi4ljgNOAzaZ5fAm6PiPnA7en6RHMx8HDJeh5yvhy4JSKOAU4gyX9C5y1pNvCXwMKIOI5k\nWv7zmXh5XwOcOWhb2RzT/8fPB16bfud76c+8MclFUQBOAdZGxBMRsQdYCpxb45gOuojYGBEr0+Vt\nJD8kZpPk+sP0sB8C765NhNmQNAf4c+AfSzZP9JynA28BfgAQEXsiYisTPO9UAzBJUgPQAvyRCZZ3\nRNwJvDho81A5ngssjYjdEfEksJbkZ96Y5KUozAbWlayvT7dNWJLmAScC9wDtEbEx3fUs0F6jsLLy\nHeBSoL9k20TP+UjgeeCf0ttm/yhpMhM874jYAHwTeAbYCHRHxL8ywfNODZXjQf35lpeikCuSpgA/\nBz5f+h5sgPTNdhPmkTNJ5wCbImLFUMdMtJxTDcBJwPcj4kRgO4NumUzEvNP76OeSFMXDgcmSPlJ6\nzETMe7Asc8xLUdgAdJSsz0m3TTiSCiQF4ScRcX26+TlJs9L9s4BNtYovA28E3iXpKZLbgm+V9H+Z\n2DlD8tvg+oi4J13/GUmRmOh5vw14MiKej4i9wPXAG5j4ecPQOR7Un295KQr3AvMlHSmpkaRTZlmN\nYzroJInkHvPDEfHtkl3LgI+lyx8DflHt2LISEV+OiDkRMY/kv+uvIuIjTOCcASLiWWCdpKPTTX8G\nPMQEz5vkttFpklrSf+9/RtJ3NtHzhqFzXAacL6lJ0pHAfOD3Y75KROTiA5wN/DvwOPA3tY4noxzf\nRNKkXA2sSj9nA68geVrhMeDfgBm1jjWj/DuBG9PlCZ8zsAC4L/3vfQPQlpO8vwY8AjwI/Bhommh5\nA9eS9JnsJWkVfnK4HIG/SX+2PQqcdSDX9ohmMzMrysvtIzMzq4CLgpmZFbkomJlZkYuCmZkVuSiY\nmVmRi4LlgqQuSQf9fbaSLktn7Lxs0PZ3DczGK+ndB3MCRkkLJJ1d7lpmB8qPpFouSOoCvhgR9x3k\n83aTPC/eN8wx15CMn/jZKM7bEBG9Q+y7kGSW0M+OMlyzEbmlYIcMSfPS9wJclf72/a+SJqX7ir/p\nS5qZTmuBpAsl3ZDOL/+UpM9K+kI6SdzdkmaUXOICSavSefhPSb8/OZ27/vfpd/abPVeJy9LvPSDp\nA+n2ZcAUYMXAtpLvXCjpSklvAN4FXJZe+9Xp5xZJKyT9WtIx6XeukbRY0j3ANySdIul3aVx3STo6\nHZH/deAD6fk+MHCtkr/DX0laLel2SXNLzn1Fep4nJJ2Xbp8l6c6Sv5c3H6T/nDZe1Xrknj/+DHyA\neSTvhFiQrv8U+Ei63EXy2zHATOCpdPlCkqmCpwKHAd3Ap9N9/5tkUsCB71+VLr8FeDBd/oeSa7SS\njHqfPCiu9wG3kczd304y1cKsdF/PELlcCFyZLl8DnFey73Zgfrp8KsnUHAPH3QjUp+vTgIZ0+W3A\nzwefu8y1fgl8LF3+BHBDybmvI/lF8FiSqeQBLiEd4Z/mN7XW/w78qe2nYdiKYVZ9T0bEqnR5BUmh\nGMnySN4fsS29nfPLdPsDwPElx10LyVz1kqZJagX+A8mEel9Mj2kG5rLvC3veBFwbyS2i5yTdAbye\nMcyflc5g+wbgumTqHiCZpmHAdfHyrajpwA8lzSeZvqRQwSVOB96bLv8Y+EbJvhsioh94SNLAtMv3\nAlenEyneUPJ3bznlomCHmt0ly33ApHS5l5dvdzYP853+kvV+9v03PrgDLQAB74uIR8ca8CjVAVsj\nYsEQ+7eXLP89ScF7j5L3Y3Qd4LVL/54ExQL5FpKXFF0j6dsR8aMDvI6NY+5TsPHiKeDkdPm8MZ5j\noC/gTSQvZ+kGbgU+l864iaQTy3zv1yT38OslHUZy+2k0s1BuI7m9RSTvt3hS0vvT60nSCUN8bzov\nT4F8YbnzlXEXyWyxAB9OYx+SpCOA5yLiKpI31500bCY24bko2HjxTeAiSX8g6VMYi13p9xeTzDoJ\nyW/jBWC1pDXp+mD/TDIT6f3Ar4BLI5m6ulJLgb9KO4xfTfLD+pOS7gfWMPSrYb8B/I805tIWz3Lg\n2IGO5kHf+RzwcUmrgQtI3l09nE7g/vQaHyB577PlmB9JNTOzIrcUzMysyEXBzMyKXBTMzKzIRcHM\nzIpcFMzMrMhFwczMilwUzMysyEXBzMyK/j/+N/QRFtWIWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11a7c04a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# lets plot the objective values of the function\n",
    "# to make sure it has converged\n",
    "plt.plot(func_vals)\n",
    "plt.ylabel('function values')\n",
    "plt.xlabel('numbe of iterations')\n",
    "\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compressive sensing application: tomography reconstruction with the Lasso\n",
    "\n",
    "\n",
    "This example shows the reconstruction of an image from a set of parallel projections, acquired along different angles. Such a dataset is acquired in computed tomography (CT).\n",
    "\n",
    "Without any prior information on the sample, the number of projections required to reconstruct the image is of the order of the linear size l of the image (in pixels). For simplicity we consider here a sparse image, where only pixels on the boundary of objects have a non-zero value. Such data could correspond for example to a cellular material. Note however that most images are sparse in a different basis, such as the Haar wavelets. Only l/7 projections are acquired, therefore it is necessary to use prior information available on the sample (its sparsity): this is an example of compressive sensing.\n",
    "\n",
    "The tomography projection operation is a linear transformation. In addition to the data-fidelity term corresponding to a linear regression, we penalize the L1 norm of the image to account for its sparsity. **The resulting optimization problem is called the Lasso**.\n",
    "\n",
    "\n",
    "The reconstruction with L1 penalization should give a result with zero error (all pixels are successfully labeled with 0 or 1), even if noise was added to the projections. In comparison, an L2 penalization (sklearn.linear_model.Ridge) produces a large number of labeling errors for the pixels. Important artifacts are observed on the reconstructed image, contrary to the L1 penalization. Note in particular the circular artifact separating the pixels in the corners, that have contributed to fewer projections than the central disk.\n",
    "\n",
    "\n",
    "<span style=\"color: red; font-weight: bold\">Exercise (10 min). I have done the image reconstruction using Ridge regression, do it using the lasso algorithm previously coded.</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/fabian/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:64: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "/Users/fabian/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:68: DeprecationWarning: numpy boolean subtract, the `-` operator, is deprecated, use the bitwise_xor, the `^` operator, or the logical_xor function instead.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import sparse\n",
    "from scipy import ndimage\n",
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.linear_model import Ridge\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def _weights(x, dx=1, orig=0):\n",
    "    x = np.ravel(x)\n",
    "    floor_x = np.floor((x - orig) / dx)\n",
    "    alpha = (x - orig - floor_x * dx) / dx\n",
    "    return np.hstack((floor_x, floor_x + 1)), np.hstack((1 - alpha, alpha))\n",
    "\n",
    "\n",
    "def _generate_center_coordinates(l_x):\n",
    "    X, Y = np.mgrid[:l_x, :l_x].astype(np.float64)\n",
    "    center = l_x / 2.\n",
    "    X += 0.5 - center\n",
    "    Y += 0.5 - center\n",
    "    return X, Y\n",
    "\n",
    "\n",
    "def build_projection_operator(l_x, n_dir):\n",
    "    \"\"\" Compute the tomography design matrix.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "\n",
    "    l_x : int\n",
    "        linear size of image array\n",
    "\n",
    "    n_dir : int\n",
    "        number of angles at which projections are acquired.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    p : sparse matrix of shape (n_dir l_x, l_x**2)\n",
    "    \"\"\"\n",
    "    X, Y = _generate_center_coordinates(l_x)\n",
    "    angles = np.linspace(0, np.pi, n_dir, endpoint=False)\n",
    "    data_inds, weights, camera_inds = [], [], []\n",
    "    data_unravel_indices = np.arange(l_x ** 2)\n",
    "    data_unravel_indices = np.hstack((data_unravel_indices,\n",
    "                                      data_unravel_indices))\n",
    "    for i, angle in enumerate(angles):\n",
    "        Xrot = np.cos(angle) * X - np.sin(angle) * Y\n",
    "        inds, w = _weights(Xrot, dx=1, orig=X.min())\n",
    "        mask = np.logical_and(inds >= 0, inds < l_x)\n",
    "        weights += list(w[mask])\n",
    "        camera_inds += list(inds[mask] + i * l_x)\n",
    "        data_inds += list(data_unravel_indices[mask])\n",
    "    proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))\n",
    "    return proj_operator.toarray()\n",
    "\n",
    "\n",
    "def generate_synthetic_data():\n",
    "    \"\"\" Synthetic binary data \"\"\"\n",
    "    rs = np.random.RandomState(0)\n",
    "    n_pts = 12.\n",
    "    x, y = np.ogrid[0:l, 0:l]\n",
    "    mask_outer = (x - l / 2) ** 2 + (y - l / 2) ** 2 < (l / 2) ** 2\n",
    "    mask = np.zeros((l, l))\n",
    "    points = l * rs.rand(2, n_pts)\n",
    "    mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1\n",
    "    mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)\n",
    "    res = np.logical_and(mask > mask.mean(), mask_outer)\n",
    "    return res - ndimage.binary_erosion(res)\n",
    "\n",
    "\n",
    "# Generate synthetic images, and projections\n",
    "l = 32\n",
    "proj_operator = build_projection_operator(l, l / 5.)\n",
    "data = generate_synthetic_data()\n",
    "y = proj_operator.dot(data.ravel()[:, np.newaxis])\n",
    "y += 0.15 * np.random.randn(*y.shape)\n",
    "\n",
    "# Reconstruction with L2 (Ridge) penalization\n",
    "rgr_ridge = Ridge(alpha=0.2)\n",
    "rgr_ridge.fit(proj_operator, y.ravel())\n",
    "rec_l2 = rgr_ridge.coef_.reshape(l, l)\n",
    "\n",
    "l1_rec, _ = lasso_ADMM(proj_operator, y.ravel(), lam=0.01)\n",
    "l1_rec = l1_rec.reshape(l, l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAADuCAYAAAB4Z3AVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUlfWd5/Hvl71Yi31fVHZkURHRaCLS6EjUqDkaHU80\n6cGOnZOYMU5r8ke6M0lr7EyWzubY5swk0z2TTCSOo8moUWJjRFBAWQRZlaXYd4Fit37zx72kS+r7\nKbhUUfUreL/OyTnmw1P3ee69z731q6fup76eUjIAAAA0rmaNfQAAAABgUQYAAJAFFmUAAAAZYFEG\nAACQARZlAAAAGWBRBgAAkAEWZQ3M3Z9w92/U97YnuZ1B7p7cvYX496XufnVd94Nz08nO0+K5N7gh\njylX9fWaBnB2cv5O2dnP3QeZ2Roza5lSOta4R4Omxt3XmllPM/vQzPab2Ytm9qWU0v5T/PpkZkNS\nSqvP2EECwFmAK2UNyN2bN/YxAKfpxpRSezMbZ2YXmdnXG/l4PsIL6u39rL5vDwBOBW86deTuI9x9\nprvvKf4a8KZq//ZLd/+v7v68u1ea2aRi9vfVtnnI3Te7+yZ3n1b9Vz3Vt3X3q919g7s/6O7bil/z\n+Wq380l3X+Due929wt2/WcJ9WOvuf1H872+6+3R3/5/uvs/d33H3oe7+9eJ+K9z92mpf+3l3X1bc\n9n13/8IJt13b/Wvt7t9z9/XuvrX4q52yUp8DNJyU0hYz+4MVFmdm9tHztPj//6bac/6X1b/e3bu6\n+++K5+k8d/97d59V7d+Hu/vL7r7L3Ve4++3qWIqvu0fc/XUzO2Bm57t7J3f/b8X9byzefvPi9s3d\n/fvuvsPd17j7l6r/Wv80bm+wu7/q7h8Ub/M3xdzd/YfF18ve4mvoQvFY3evuq4v39zl371Pt35K7\n3+fuq4rvLz9zdz+d5w1A08CirA7cvaWZ/c7MXjKzHmb2ZTP7X+4+rNpm/97MHjGzDmY264Sv/3dm\n9lUz+wszG2xmV59kl73MrJOZ9TWz/2BmP3P3zsV/qzSzu82s3Mw+aWZ/7e43n+Zdu9HM/sXMOpvZ\nAit8E25W3O+3zOyfqm27zcxuMLOOZvZ5M/uhu198ivfvMTMbaoVv8IOLt/+3p3nMaADu3s/Mrjez\n8FeRxef8P5nZFDMbYoXnvrqfWeFc7WVm9xT/d/xr25nZy2b2Kyu8nu4ws8fdfWQth/RZM/srK7y+\n1pnZL83smBXOp4vM7Fozm1bc9t7isY8zs4vNLHp9lHJ737bCa7+zmfUzs58U82vN7ONWOLc7mdnt\nZrbzxB25+zVm9p3iv/cu7u9/n7DZDWZ2qZmNKW53XS2PBYAmjkVZ3Uw0s/Zm9lhK6UhK6RUz+72Z\n3Vltm2dTSq+nlKpSSodO+PrbzewXKaWlKaUDZvbNk+zvqJl9K6V0NKX0vBU+3zPMzCylNDOl9E5x\nP4vN7Ndm9onTvF+vpZT+UPz82XQz6168j0et8E1jkLuXF/f7/1JK76WCV63wTeqqk92/4k/8f2Vm\nD6SUdqWU9pnZo1b4Roz8/F9332dmFVZYiP+d2O74c74kpVRpH33Om5vZp83s71JKB1JK75rZ/6j2\ntTeY2dqU0i9SSsdSSgvM7Gkzu62W4/pl8fw6ZmZdzGyqmf3HlFJlSmmbmf3Q/u2cut3MfpRS2pBS\n2m2FHwrqcntHzWygmfVJKR1KKc2qlncws+FW+NzuspTS5mBfd5nZf08pvZ1SOmyFXwlf7oXPgB73\nWEppT0ppvZn9q1W7Qgng7MOirG76mFlFSqmqWrbOCld8jqs42def4rZmZjtP+KD+ASssCs3dL3P3\nf3X37e7+gZndZ2bdTnYHhK3V/vugme1IKX1Y7f9btf1e7+5vFH/9sscK38SO77e2+9fdzNqa2VvF\nX83sscIHyLuf5jHjzLo5pdTBClc7h5s+t058ztdV++/uZtbC9Dkx0MwuO34+FM+Ju6xwVU058etb\nmtnmal//T1a46hYdW/R6K+X2HjIzN7O5Xvjowl+amRV/OPupFa4KbnP3J929Y7CvPlbt8SkWJ3ba\nR98/tlT77z+/3gGcnViU1c0mM+vvH/1A8AAz21jt/9dWb91shV97HNe/DsfyKzN7zsz6p5Q6mdkT\nVviGcca4e2srXMn4npn1TCmVm9nz1fZb2/3bYYUF3qiUUnnxf52KHyZHpopXQ39phec8stk++jwP\nqPbf263wq0B1TlSY2avVzofylFL7lNJf13ZIJ3z9YTPrVu3rO6aURlU7tpO93k759lJKW1JK96aU\n+pjZF6zwq9bBxX/7cUrpEjMbaYVfY/5NsK9NVlj4mdmff33b1T76/gHgHMKirG7etMJPrw+5e0sv\n/K2vG63m50KUp8zs814oC7Q1s7r8/aIOZrYrpXTI3SdY4bNsZ1orM2ttxW+27n69FT5Pc5y8f8Wr\niz+3wmfQepiZuXtfd+czM/n7RzOb4u5jg397ysw+5+4ji8/5n3/NWbza+n/M7Jvu3tbdh1vhc5DH\n/d7Mhrr7Z4uvp5bufqm7jziVgyr+ivAlM/u+u3d092bufoG7H/81/lNm9pXieVZuZg/X5fbc/bbi\nZ+zMzHZbYUFXVTzmy4qfOa00s0NmVhXs4tdWeH2MK/6A86iZvZlSWnsq9xfA2YdFWR2klI5YYRF2\nvRWu/DxuZnenlJaf4te/YGY/tsJnRVab2RvFfzp8GofzRTP7VvFzP39rhW9AZ1Txc2D3F/e12woL\nweeq/fvJ7t/Dx3N332tmM6z4GTnkK6W03cz+2YJSRvE5/0cze8UKz+0rJ2zyJSt8+H2LFcokv7bi\n+VA8n661wme2NhW3+QcrLPxP1d1W+GHhXSuck7+1wofozQo/BLxkZoutUGB53gpX7j6seTOndHuX\nmtmb7r7fCuf9V1JK71uh9PLz4vbrrPAryf9y4g2nlGZY4QeVp61wFe8C4zOVwDmNPx6bkeIVgSVm\n1vps/COvZ/v9Q+nc/R/MrFdK6Z6Tblz/+77ezJ5IKQ086cYA0AC4UtbI3P0WL/y9rs5WuCrwu7Np\nwXK23z+Uxgt/h2xM8W95TbDCn3Z5poH2XebuU929hbv3tcKvVhtk3wBwKliUNb4vWOFPDLxnhV+j\n1Pah5qbobL9/KE0HK3yurNLMfmNm3zezZxto325m/9kKv1ZcYGbLjL+LByAj/PoSAAAgA1wpAwAA\nyACLMgAAgAy0aMiduTu/K8U5KaXkZmZTp04NXwM33xyPKZ09e3aYN2sW/zzVvHnzMO/VK/6j+O+8\n806YHzhwIMzHjBkT5keOHAnzzZuj6UI637JlS5i3bNkyzD/44IMw79ix5h/QHzkyHqG5bt26MF+0\naFGY9+8f/43nIUOGhPnYsdGfczN76623wryiIh7sMXr06DC/+OKLw7yysjLM1UdW2reP/25zVVX0\nJ9bM1q9fH+Y9evQI80cffdT5HoBz1fHvASfDlTIAAIAMsCgDAADIAIsyAACADLAoAwAAyACLMgAA\ngAw0aPsSONeplmKbNm3CfP/+/WHeu3fvMD98OJ5lr5p4qr3YvXv3MFeNO5Wr1uT48ePDfM+ePWG+\nYsWKMN+7d2+Yl5WV1ciiRqaZWatWrcL8vPPOC3P1XO3bty/MlyxZEuYtWsRvvyofODAe0aluv0uX\nLmGunlvVslSPcXl5eZirBi2Ak+NKGQAAQAZYlAEAAGSARRkAAEAGWJQBAABkgEUZAABABmhfAg1I\ntSCPHTsW5mqWpWr6qcbd6tWrw3zEiBFhPm/evDDftm1bmHfr1i3Mhw0bFuabNm0K8wEDBoS5mk+p\nto9mfao25fLly8NctTIPHjwY5qq92K9fvzDv0KFDmK9atSrMVVN269atYa5mWaq5pkePHg3zTp06\nhfn27dvDXDVoAZwcV8oAAAAywKIMAAAgAyzKAAAAMsCiDAAAIAMsygAAADJA+xLIwKFDh8JcNQBV\ng65169Zhrppy7dq1O4Wj+zc9evQI80GDBoW5mpuo2qNt27YNczVDc8yYMWE+ZMiQGtnatWvDbdUs\ny5RSmKtW44cffhjmqmnatWvXMP/c5z4X5qqBumbNmjBX59TgwYPDXN2vnj17lrRf9bih4UTntDof\nkBeulAEAAGSARRkAAEAGWJQBAABkgEUZAABABliUAQAAZID2JdCAOnfuHObNmsU/H6l5hKq9OGPG\njDC/+uqrw1w1D/v06RPmaq6hah5u2bIlzFWLU81Z7NixY5irVmnU1lTNTjX/s6KiIszVDEr12Kjn\n6vzzzw/zaG6nmVlZWVmYb968Ocznz58f5itWrAhz1frcv39/mKv7pRrAaDhNuWmpmtaqKX62zVrl\nShkAAEAGWJQBAABkgEUZAABABliUAQAAZIBFGQAAQAZoXwINaPjw4WG+a9euMF+1alWYq8adaiip\npp9qEm7dujXMVWtSzXccMGBAmHfp0iXMFTVPUTUqlyxZUiObMGFCuO3BgwfDXM2CrKqqCvNSG2/d\nu3cP840bN4a5Ov7Zs2eHuWqyqrba8uXLw3z9+vVh7u5hrmZrAqdCNc779esX5rQvAQAAUO9YlAEA\nAGSARRkAAEAGWJQBAABkgEUZAABABmhfAg1IzRFUVCvz8OHDYf6JT3wizDds2BDmqkGnWpmqxanm\nHarb37FjR0m3rx43NXMzytU+VbtQPcaqNTlw4MAwb9WqVZiruaDl5eVhrtqUQ4cODfO33347zNX8\n1SNHjoR5165dw3zbtm1hrmYXAnXRosW5sVzhShkAAEAGWJQBAABkgEUZAABABliUAQAAZIBFGQAA\nQAZczZQ7Iztzb7idARlJKbmZ2XnnnRe+Bp588snw6+64444wLysrC/MxY8aE+YEDB8K8W7duYa4a\ngKrFqWZQqhmdvXv3DvOxY8eGuZpvt3DhwjCPZoCqWZa7d+8OczWDUjVN1WMwadKkMFf3Sd2+mn2p\nnpOZM2eGuZrRqdqmF154YZirc0q1PhcuXOjqe4Bq6Tbk9yfkTb1XqVZ1bo5/DzgZrpQBAABkgEUZ\nAABABliUAQAAZIBFGQAAQAZYlAEAAGTg3BgmBWRCNYhUA1DNe1NNQjWzsn///mGuWpyqoadmUKq5\niep2qqqqwlzNslyxYkWYqwZg3759a2RqBmVFRUWYHzt2LMzVc6Wao+vWrQvz5s2bh3nUHK1te9Vc\nvO2228L81VdfDfP33nsvzF977bWS9qvOqdrQssTJnM551RRxpQwAACADLMoAAAAywKIMAAAgAyzK\nAAAAMsCiDAAAIAO0L5GVM93CUo2xhjJixIgw37t3b5g/8sgjYa4acWr+YteuXcNczbicO3dumKs5\niGqO46hRo8L84MGDYa506NAhzFWrtEuXLjWyqJFpptuUHTt2DPM2bdqEubqv27ZtC3N1PKp92aNH\njzC/+uqrwzx6DMzMli5dGuYrV64M85EjR4b51q1bw1yda7W58sorw/z1118P8/p6n+jevXuYb9++\nvV5uv7FETUX1WlHP+5mmjkc9Jw888ECYT548OczVe/3ll18e5mvXrg3zhsaVMgAAgAywKAMAAMgA\nizIAAIAMsCgDAADIAIsyAACADGTdvlQNm8Zu0KHuGuu5LbW1Vd/Hs2bNmjCfOHFimE+YMCHMe/bs\nGebPPPNMmKuZlWrWpJopqZqBO3fuDHM1E1M1r9ScyKNHj4a5aq1GcyJVe1G1LFUrU7UsN27cGOaq\nTVZfs/wWLVoU5qpNptqUqhGrZmKqmaGtW7cO89qoWag/+tGPwly1j6dPn17Sfpt6y1KJnsvGalkq\nahateo98++23w3zhwoVh/sILL4T5LbfcEubz5s0L81mzZoX5mcKVMgAAgAywKAMAAMgAizIAAIAM\nsCgDAADIAIsyAACADGTRviy1iUcrs+nI7bkq9Zyqb2qGo2r0qdl/kyZNCvOLL744zFXT6f333w9z\n1cpcsmRJmA8cODDMVWtSNRtVS1Q9P4cPHw7zY8eO1cjUfVLty27duoW5mhfaokX8dnreeeeFuZqh\nuWrVqjBXLbPly5eHefv27cNcPVejR48Oc/W4qZbc6by2v/rVr4b5Y489Fua33nprmM+ePTvMe/Xq\nFeZvvfXWKRwdqlOvix07doS5eq9S70l/+MMfwvzpp58O8/Xr14e5ombOXnXVVWH+5ptvhrl6b6sr\nrpQBAABkgEUZAABABliUAQAAZIBFGQAAQAZYlAEAAGQgi/ZlqWhl5ofH/tSohl7UFjQz27BhQ5gv\nWLAgzKdOnRrmqtHUqVOnMFcNw127doW5Ov62bduG+bJly8J88ODBYa4ehz59+oR51Gw8dOhQuK1q\nF06ePDnMP/jggzBXzVf1mKlWl2qxDRo0KMzV/EY1+7JZs/hn8eHDh4e5mmV5xRVXhHnv3r3DvDZb\nt24N8/nz54f5nj17wvzRRx8N85/85CdhruaPqjmg55LrrrsuzNV5NWXKlDCfO3dumKuZqop6rZfa\nvnz55ZfDvEuXLiXl6pwdP358ScdzIq6UAQAAZIBFGQAAQAZYlAEAAGSARRkAAEAGWJQBAABkoEm2\nL9F4ztaWZUM1elW7ULUgN27cGOaqQXfJJZeEuZqzqNqganvVMNy5c2eYb9myJczV7MvmzZuXdDyq\neRi1R1WLql+/fmGu5iWq52TdunVhruaFqlbm+eefH+abN28Oc/UYqAahek7UjE7VJlu8eHGYq5Zr\nbdRzo6xevTrM1WNUVVUV5mo+6LnUvvzYxz4W5hMmTAhz1Q5WrUa1vaLeA9SM1/qiGt633XZbmE+f\nPj3M1ezaU8WVMgAAgAywKAMAAMgAizIAAIAMsCgDAADIAIsyAACADJxV7UtmYpZOPTYKj1ndqJZZ\n165dw1w18VTD7emnnw7zm266KcyfffbZMFfPc8eOHUvaXs0oHDVqVJirFmfLli3DXLVHO3fuXCNT\nj+WIESPCfObMmWGuGn6qEajaZGouaGVlZZirBuGRI0fCvH///mGumqwqV41edX/ff//9MK/NlVde\nWfLXRL72ta+F+eOPPx7mjzzySL3stymYOHFimD/22GNh/r3vfS/M1fOuZsKWSs2iVe8l9UXN71Uz\nPWfMmBHmdW2JcqUMAAAgAyzKAAAAMsCiDAAAIAMsygAAADLAogwAACADZ1X7EhoN1DyoBpFq7Ki5\njGr+4r59+8L8vffeC/NSW5Oqoae2V40m1TZVjUQ1J7KUWZlqluVzzz0X5uvXrw/zHj16hLl6LNUs\nPNWOXLp0aZirxqpqlZaVlYW5amtu3bo1zDdt2hTmqhGr5rKa6bbsH//4xzAfOXJkmL/77rtyHxH1\nmJbaPm/K7rrrrjD/+te/HuZz5swJ84suuijM58+ff3oHlgl1/qs2sTqX64orZQAAABlgUQYAAJAB\nFmUAAAAZYFEGAACQgXPig/6MX0Iu1IdJ1Yev1WidVq1ahbkaUaI+bK5GiCxatCjM16xZE+aqYNCt\nW7cwb926dZirDxGr8UtK9OHcd955J9x2x44dYd6pU6cwVx98P3jwYJirD8Sr50qVKcrLy8NcUeWI\nAwcOhLm6X+pD9erxUWOrzHTxQ91nNdarVG+88UaYf+c73wnze+65p172W6px48aF+cKFC+t82+p5\nXLduXZir74O5faBflYZKfc9QxZ65c+eGuSrq1LU8wpUyAACADLAoAwAAyACLMgAAgAywKAMAAMgA\nizIAAIAMnBPty3NN1P6gUXp66ru5q1qH7du3D/PFixeHuRr1c+jQoTBX44gqKyvDfOrUqWH+05/+\nNMzV/VYNJdWq69mzZ5ir1ueKFSvC/FOf+lSYR1QLcvfu3WGuWo3qsR87dmyYq1FZ6jFQ7UvVrN2+\nfXuYq/urzkF1nGPGjAlz1QyuTffu3cNctSZLtXLlyjBXbdDevXuH+ebNm+vleJSjR4/Wy+3cd999\nNTJ1XysqKupln42l1Jaloprub731VphPmDAhzAcNGlSn4+BKGQAAQAZYlAEAAGSARRkAAEAGWJQB\nAABkgEUZAABABmhfAqeh1FbmcaqZ1rlz5zBXcwrVDM0OHTqEuWrcLV26NMxHjhwZ5hMnTgxzpays\nLMzVvD3VBGvbtm2Yq3Zn1CRU7Sp1G2r+5+HDh8N89OjRYa6OXT0G6hxSz6FqWapWmpoVqKjmrmpl\n9urVq6TbN6u/GZdKs2bx9YcFCxaE+bBhw8L8TLcv1euxVFGbdcmSJeG29dU0bdEiXk5ccMEFYa6a\n07lRc4nXrl0b5qU0vyNcKQMAAMgAizIAAIAMsCgDAADIAIsyAACADLAoAwAAyADty7NQ1Aw83VmN\nqF9qPqLKVZNNNaPatGkT5qqVqVpp8+fPD/M77rgjzJ955pkwVw0rNeNPNbi2bNkS5tdcc02Yv/vu\nuzWyyy67LNxWtQhHjBgR5qqZqrzzzjthrmYR7t+/P8w3bdoU5qp92aVLlzBXTUfV9FXzWl977bUw\n/8Y3vhHmtTnTrUbVMFRN1JkzZ4b58OHDw3z58uWndVx11alTpzD/8Y9/XCP7zGc+E26rmtx9+/YN\nc/UaVS3mWbNmhXlToV4X6n49/PDDddofV8oAAAAywKIMAAAgAyzKAAAAMsCiDAAAIAMsygAAADJw\nTrQvaR5qJ5vVeKacq4/9+eefH+ZDhgwJ80suuSTM1YzIQYMGhfmYMWPCXLUjV61aFeZq/mK3bt3C\nXLUsVQOwXbt2YX7w4MEwj2ZcmsUzRidMmBBu27NnzzBXbUfVylyzZk2Yq9l5qjmnmqaDBw8Oc9WO\nVDMrVetTnYPqPULN+qytndqyZcswv/HGG8P8iiuuCPPZs2eH+eTJk8P8Bz/4QZg/+eSTYa6o862x\nqLmq0Tnx4IMPhttefPHFYT59+vQw/+IXvxjmaparOsaoIW1mtmHDhjA/09RcX3X+79u3L8wHDBhQ\np+PgShkAAEAGWJQBAABkgEUZAABABliUAQAAZIBFGQAAQAa8Idt37h7u7Ey3I2lf5udsfU5qeT25\nmdm4cePCDe6777562U95eXmYDxs2LMyfeuqpMH/hhRfCXDUVP/7xj4f5b3/72zAvKysLc9Xo69y5\nc5irduell15aI1NzQSsrK8NczZRUTdPVq1eHuZqtqeYxqu3VTMzDhw+H+QcffFDS9qpV+v7774e5\nmoOqZiNOmzbN1feAL3/5y+HXzJkzJ8zVY6HmNX76058Oc9UkbCrU++U999xTI1Pnw3vvvRfmixcv\nDvNx48aF+bJly8JcnW9XXnllmDeVWZkXXnhhmKv26O7du0/pmxtXygAAADLAogwAACADLMoAAAAy\nwKIMAAAgAyzKAAAAMpDF7EvVIKmvZuhZ3Oirs8Z6bEp9znN7Dks9zuPb9+rVK/x3lav5izfffHNJ\n+3/xxRfDfOjQoWG+Z8+eMFczKKNZk2ZmkyZNCnM19+6aa64Jc9W+VMaPH18jU03NV155JcxVk0/N\nHVXbq0apmgOpWpBqhuDbb78d5kuWLAlz1dBVx3PnnXeGubpfqs1aGzVTcteuXWH+7W9/O8x79OgR\n5lOnTg1zNR9UzXg901SL+aabbgpz1dS9++67a2Rqpuq0adPCvEOHDmG+cOHCMFezWdX80oEDB4a5\nOq/Ue+4TTzwR5vXllltuCXPVwlaN4VPFlTIAAIAMsCgDAADIAIsyAACADLAoAwAAyACLMgAAgAxk\n0b5UcmvcnWlnsnnYkDNO66KptzJPRjXcdu7cGea///3vw3zVqlVhrtpkR44cCXPV+ty6dWuYr127\nNsxXrlwZ5tdff32Yd+nSJczVfL7169eHuZrFOWPGjBqZauZt3rw5zNVjpmYCtm7dOsxHjhwZ5ocO\nHSppvxdddFGYr1mzJsyrqqrCXD1Xu3fvDvPPfvazYf7aa6+FuTqXp0yZEuZmegarmrt51113ydtq\nytQ81OnTp4e5ajA+8sgjNbILLrgg3FY9X+r8VNTsV/UeM3v27DBXze/777+/pOOpL7feemuYq9dF\nXXGlDAAAIAMsygAAADLAogwAACADLMoAAAAywKIMAAAgA1m3L1F/zvZW45lS34+PmvG3b9++MFez\nL/fu3Rvm7dq1C3M1x07Nw+vUqVOYd+/ePcxV61M1uL7yla+EuZrjqGZuqscn2l49Z2rW5KJFi8Jc\ntdhUa0w1VtW5VVZWFuaKmompnsM2bdqE+e233x7mqlW3YsWKMB8+fHiY12bdunUlf01T9sADD4S5\nOlfUc6AsW7asRvbJT34y3Pbee+8N85///OdhXmor8/XXXw/zUv8igHoPU+e5anIrN9xwQ5jPmjWr\npNupK66UAQAAZIBFGQAAQAZYlAEAAGSARRkAAEAGWJQBAABkgPZlI6DxWLqzpT0ataLMzPr27Rvm\nam6cammpxpFqfR47dizMW7SI3xrU/EVFNaPUjD81E3PHjh1hXlFREeatWrWqkanHYP/+/WGu2pTq\nGNXsSzVHtLKyMsxV01TNgVTt0Y0bN4a5amuOHz8+zLdt2xbml19+eZirWZy1+fDDD0v+mqZsyJAh\nYa7axGo+pcojzZrF12DUea5alur8j15zZrrhXapXXnklzK+77rowV/NUlSuuuCLM33zzzTCfOHFi\nmL/xxhsl7fdEXCkDAADIAIsyAACADLAoAwAAyACLMgAAgAywKAMAAMgA7ctG0BhNwqbWUjzRmT7+\nhnp8Ro8eHebXXHNNmL/77rthruYOqramauhVVVWF+dKlS8NcNfrU7MuhQ4eG+Zw5c8JcNQAXLFgQ\n5qrx2L9//xqZuq8LFy4Mc9VkVbejjuXAgQNh3rlz5zBXsylVu03NQZ03b16Y33nnnSUdj2rnqVZp\ny5Ytw7w2AwcODHPVMm4qevToEea/+tWvwvzBBx8M89/85jcl7bdnz541MtWaVDMuld27d4d5eXl5\nSbejqCZ6x44dw1zN1lSuvfbaMO/Xr1+YP/vss2GuWpl1xZUyAACADLAoAwAAyACLMgAAgAywKAMA\nAMgAizIAAIAM0L5sAkptBqrt0fgtVDV/UTWLVCPo+uuvD3PVjFL7Vbd/+PDhMFfzFy+55JIwnzVr\nVkn7veGF/lKMAAAEXUlEQVSGG8JcNaZmzJgR5gMGDKiRqYaoakeqGZGlzqxU8z9V41DNgVTb/+lP\nfwrzqVOnhvmoUaPCXDViVZtSzVLs2rVrmNdGzWBt6tTcUJVPnz49zB966KEw/+53vxvmU6ZMqZGp\nJqiacamo99B9+/aFuTpPVOO8e/fuYa5mTaoW5GWXXRbmkydPDvP7778/zEs9N6Pmdym4UgYAAJAB\nFmUAAAAZYFEGAACQARZlAAAAGWBRBgAAkAHalxkptQHY2E3CnOX62Kh5gevXrw/z7du3h/mwYcPC\nvKysLMy3bNkS5mvWrAlz1bI8cuRImKv5iGp2oWoYPv7442E+duzYMFdzIqP7q/apmnBq1mSpj8Hw\n4cPDXLUdVVttzJgxYa6oOaIzZ84MczVndfDgwWGuzk3VJK7Nhg0bSv6as9GkSZPCXDUGVaNy5cqV\nNbJVq1ad/oGdAtVS7NOnT5ir5rRy8803l3T7EyZMCPOHH364pP2WqqKiok5fz5UyAACADLAoAwAA\nyACLMgAAgAywKAMAAMgAizIAAIAM0L5swhq7SZiDXFuWimr0VVVVhXnnzp1Lun01x07NL1THo5pU\nao6japWOHDkyzFUb9Pnnnw/z5cuXh7lqoUbz8EaMGBFuq+6TOkY1m0+1L5cuXRrmqg2qWpZz584N\n83bt2oW5ek5Ue7Rv375hrtqpqrl7Ou1LFKjzXOWqYRudK+PGjQu37dWrV5ir879UmzZtKmn7QYMG\nhfmLL74Y5uvWrQtzdf7njitlAAAAGWBRBgAAkAEWZQAAABlgUQYAAJABFmUAAAAZoH2JJiHXNmWp\nVPNt7969Ya7alKrpp+YmqnanatZdddVVYd6tW7cwf/nll8NczdBs3rx5mO/atSvMN27cGOZq9uW8\nefNqZPv37w+37devX5irc27BggVhrtqLqslaXl4e5q+++mqYt23bNsynTZsW5qqtefTo0TAvtVWq\nnsP6au3h5ObPn3/K2y5cuPAMHkn9UfNyf/GLX4S5aqirVmbuuFIGAACQARZlAAAAGWBRBgAAkAEW\nZQAAABlgUQYAAJAB2pdAA1JNNmXfvn1hfvDgwTBXsyxVs05R8wvVDM0pU6aE+eLFi8NczcNTTcX2\n7duHedSyNItnYqrH8qWXXgrzFi3it0c1z3PgwIFhPmfOnDCvrKwMc9VAnTp1apirGZQ7duwIc9Wa\nVM+5mo24ffv2MFctVOBM2L17d2MfQr3iShkAAEAGWJQBAABkgEUZAABABliUAQAAZIBFGQAAQAY8\npdTYxwAAAHDO40oZAABABliUAQAAZIBFGQAAQAZYlAEAAGSARRkAAEAGWJQBAABkgEUZAABABliU\nAQAAZIBFGQAAQAZYlAEAAGSARRkAAEAGWJQBAABkgEUZAABABliUAQAAZIBFGQAAQAZYlAEAAGSA\nRRkAAEAGWJQBAABkgEUZAABABliUAQAAZIBFGQAAQAZYlAEAAGSARRkAAEAG/j9TFAnmsjQLdwAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11ada2978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.figure(figsize=(8, 3.3))\n",
    "plt.subplot(131)\n",
    "plt.imshow(data, cmap=plt.cm.gray, interpolation='nearest')\n",
    "plt.axis('off')\n",
    "plt.title('original image')\n",
    "plt.subplot(132)\n",
    "plt.imshow(rec_l2, cmap=plt.cm.gray, interpolation='nearest')\n",
    "plt.title('Ridge regression')\n",
    "plt.axis('off')\n",
    "plt.subplot(133)\n",
    "plt.imshow(l1_rec, cmap=plt.cm.gray, interpolation='nearest')\n",
    "\n",
    "plt.axis('off')\n",
    "\n",
    "plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,\n",
    "                    right=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}