ukmlv/Compare-ISTA-FISTA.ipynb

## Compare-ISTA-FISTA.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing matrixmodel.py\n"
     ]
    }
   ],
   "source": [
    "%%file matrixmodel.py\n",
    "\n",
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "This file defines a class for matrix forward operators that are useful for \n",
    "iterative reconstruction algorithms.\n",
    "\n",
    "Created on Tue Jul 26 2016\n",
    "\n",
    "@author: U. S. Kamilov, 2016.\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "class MatrixOperator:\n",
    "    \"\"\"\n",
    "    Class for all matrix forward operators\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, H_matrix, sigSize):\n",
    "        \"\"\"\n",
    "        Class constructor\n",
    "        \"\"\"\n",
    "        \n",
    "        self.H_matrix = H_matrix\n",
    "        self.sigSize = sigSize\n",
    "        self.L = np.linalg.norm(H_matrix, 2) ** 2\n",
    "\n",
    "\n",
    "    def apply(self, f):\n",
    "        \"\"\"\n",
    "        Apply the forward operator\n",
    "        \"\"\"\n",
    "\n",
    "        fvec = np.reshape(f, (f.size, 1))\n",
    "        z = np.dot(self.H_matrix, fvec)\n",
    "        return z\n",
    "\n",
    "\n",
    "    def applyAdjoint(self, z):\n",
    "        \"\"\"\n",
    "        Apply the adjoint of the forward operator\n",
    "        \"\"\"\n",
    "\n",
    "        fvec = np.dot(np.transpose(self.H_matrix), z)\n",
    "        f = np.reshape(fvec, self.sigSize)\n",
    "        return f\n",
    "\n",
    "\n",
    "    def getLipschitzConstant(self):\n",
    "        \"\"\"\n",
    "        Return the largest eigen value of HTH\n",
    "        \"\"\"\n",
    "        \n",
    "        return (self.L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing algorithm.py\n"
     ]
    }
   ],
   "source": [
    "%%file algorithm.py\n",
    "\n",
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "This file implements fast iterative shrinkage/thresholding algorithm (FISTA)\n",
    "\n",
    "For more details on this code see corresponding paper:\n",
    "\n",
    "U. S. Kamilov, \"Minimizing Isotropic Total Variation without Subiterations,\" \n",
    "Proc. 3rd International Traveling Workshop on Interactions between Sparse models \n",
    "and Technology (iTWIST 2016) (Aalborg, Denmark, August 24-26)\n",
    "\n",
    "Created on Tue Jul 16 07:08:59 2016\n",
    "\n",
    "@author: U. S. Kamilov, 2016.\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def fistaEst(y, forwardObj, tau, numIter = 100, accelerate = True):\n",
    "    \"\"\"\n",
    "    Iterative reconstruction with FISTA\n",
    "    \"\"\"\n",
    "    \n",
    "    # Lipschitz constant is used for obtaining the step-size\n",
    "    stepSize = 1/forwardObj.getLipschitzConstant()\n",
    "    \n",
    "    # Size of the reconstructed image\n",
    "    sigSize = forwardObj.sigSize\n",
    "    \n",
    "    # Define iterates\n",
    "    fhat = np.zeros(sigSize)\n",
    "    s = fhat\n",
    "    q = 1\n",
    "    \n",
    "    # Tracks evolution of the energy functional\n",
    "    cost = np.zeros((numIter, 1))\n",
    "\n",
    "    \n",
    "    # Iterate\n",
    "    for iter in range(numIter):\n",
    "        \n",
    "        # Gradient step\n",
    "        fhatnext = s - stepSize*forwardObj.applyAdjoint(forwardObj.apply(s)-y)\n",
    "        \n",
    "        # Proximal step\n",
    "        fhatnext = softThresh(fhatnext, stepSize*tau)\n",
    "        \n",
    "        # Update the relaxation parameter\n",
    "        if accelerate:\n",
    "            qnext = 0.5*(1+np.sqrt(1+4*(q**2)))\n",
    "        else:\n",
    "            qnext = 1\n",
    "        \n",
    "        # Relaxation step\n",
    "        s = fhatnext + ((q-1)/qnext)*(fhatnext-fhat)\n",
    "        \n",
    "        # Compute cost\n",
    "        cost[iter] = evaluateCost(y, forwardObj, fhat, tau)\n",
    "        \n",
    "        # Update variables\n",
    "        q = qnext\n",
    "        fhat = fhatnext\n",
    "        \n",
    "    return fhat, cost\n",
    "\n",
    "\n",
    "def softThresh(w, lam):\n",
    "    \"\"\"\n",
    "    Soft-thresholding function\n",
    "    \"\"\"\n",
    "\n",
    "    # norm along axis 3 of differences\n",
    "    abs_w = np.abs(w)\n",
    "\n",
    "    # compute shrinkage\n",
    "    shrinkFac = np.maximum(abs_w-lam, 0)\n",
    "    abs_w[abs_w <= 0] = 1  # to avoid division by zero\n",
    "    shrinkFac = shrinkFac/abs_w\n",
    "\n",
    "    # save output\n",
    "    what = shrinkFac * w\n",
    "\n",
    "    return what\n",
    "\n",
    "def evaluateCost(y, forwardObj, f, tau):\n",
    "    \"\"\"\n",
    "    Evaluate the energy functional\n",
    "    \"\"\"\n",
    "    \n",
    "    cost = 0.5*np.sum((y-forwardObj.apply(f))**2) + tau*np.sum(np.abs(f))\n",
    "    cost = cost/(0.5*np.sum(y ** 2))\n",
    "    return cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%reset -f\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "from scipy import misc\n",
    "from matrixmodel import MatrixOperator\n",
    "from algorithm import fistaEst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAEKCAYAAACYBHl/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFzlJREFUeJzt3Xu03lV95/H3JwlJBIRcwAQChFsMLqJEJkYpaTWIyG0K\nldYqKogzrSMFFcZVrHZGnZmu6cyyhRZKNY5F0LpYERgu2kBkhUWJ5ZaShCQkGEg9hJAEcgNCSMjl\nO388vzAPh7P3OXlu55ydz2uts9Z5ft/97N9OVviwf89vP7+tiMDMrCRD+nsAZmat5mAzs+I42Mys\nOA42MyuOg83MiuNgM7PiONisZSQdLekVSerAufZIOr7d57HBycHWjyTNkPQrSVskbZD0kKR/19/j\nypE0QdJtkl6StFnSk5IuAYiI1RFxSHRmcaQXYFrSsP4ewP5K0juBe4AvAj8DhgO/Dexow7mGRsTu\nFnX3Y2AhcDTwBvBeYHyL+t4XbZ8V2uDlGVv/eTcQETE7anZExP0RsRRA0qWS5ku6vprRPSXpjL1v\nlvT56tgrkp6R9Md1tQ9LWi3pTyWtBf5B0lhJ91SzrI2SHqxrf0Q1C3tR0rOSrsyM+wPAzRGxPSL2\nRMTiiLiv6mdidYk4pHp9rKQHJb0saa6kGyT9uFvbSyR1Vef+Rt2YPiDpX6rxrqn+Hvw/YusTB1v/\n+TWwW9KPJJ0taVQPbT4IrATGAt8G7qhrtx44NyIOAS4DrpU0te6944FRwDHAHwP/GVhd9fUu4BsA\n1edh91CbhR0BfBT4iqSPJcb9MHCjpD+UdHQP9fpLxJ8Cj1Tn/A7wOd5+CXk6MAk4E/ivkiZXx3cD\nXwXGAKcBZwCXJ8Zk9hYOtn4SEa8CM4A9wCzgRUl3STq8rtn6iPjbiNgdEbOBp4HzqvfPiYjfVL8/\nBMyldim7127gWxGxMyJ2ADupBddxVX+/qtp9ADgsIv6iOv4b4P8An0oM/Q+Afwb+HFgl6QlJ07o3\nknQMMK0aw67qfHd3/2sAvh0Rb0TEk8Bi4JTqz/RERDxWzWafq/6OPpz8CzWr42DrRxHxdER8ISKO\nAaYARwLX1TVZ0+0tXVUbJJ0j6eHqsnIzcA5wWF3blyJiZ93r/w08C8ytLl2vqY5PBCZI2lT9bAb+\njNqsrqcxvxwR34iI9wLjqIXR/+2h6RHApojYXndsdQ/t1tf9vg04uPrzTaounddK2gL8Rbc/n1mS\ng22AiIhfAz+iFnB7TejW7BjgBUnDgduohdXhETEamMNbP1B/yyVfRLwWEV+LiBOA3wWuljSTWtis\niogx1c/oiDg0Iv59H8a8CfgucKSk0d3Ka4ExkkbWHevp0jXl74HlwAkRMQr4Jr5hYH3kYOsnkiZL\nulrShOr10cCnqX2Gtde7JF0paZikPwBOAn5B7Q7qcGBDROyRdA5wVi/nO0/SCdXLV4Fd1C6DHwNe\nrW40jJQ0VNLJPV1eVv38ZVUfWt3ZvRx4JiI2720CUF0+LgC+LekASacB3cMyF1TvBF6JiG2STgK+\nlPvzmdVzsPWfV6ndHHhU0qvAvwBPAl+ra/MotQ/WNwD/HbgoIrZExFbgy8DPJG2i9nnYXb2cbxJw\nf3WuXwF/FxEPRsQe4HxgKvBvwIvAD4BDEv0cSO3SczPwDLVZ2O/W1etnip8Bfqsa/38DbuWty1m6\n30iof/014DOSXgG+X7031dbsLeQHTQ5Mki4F/kNE/E5/j6VVJN0KLI+I7/T3WKxsnrFZ20iaJul4\n1ZxNbWZ3Z3+Py8rnBY/WTuOBO6itRXse+E8Rsbh/h2T7A1+KmllxfClqZsVp+6XoVVdd1dCUcPz4\n9Peq161bl6xdeOGFydqdd/rjnX21fPnyZO0973lPsvbhD6e/JPDggw8mawPJ/Pnzk7UZM2a0/HzX\nXXddsvbVr341WRs3blyy9vWvf72ptX/HHntsdHV19bV5V0Qc28z5WsWfsZlZUldXF3v27OlT2yFD\nhkxs83D6zMFmZlmD8XN4B5uZZTnYzKw4fb0UHUgcbGaWNRhnbG1fx9boXdGB5PDDD0/WXnrppQ6O\nJO/cc89N1jZt2pSsPfLII+0Yjg0A1157bVN3RSXF9u3be28IjBw5kogYEE9g8To2M8uKiD799ETS\nUZLmSVomaYmkL6fOUz0OfqekTzQ7Zl+KmllWk1d1u4CrI2KRpIOBf5U0NyJW1Deq9sn4S+C+Zk62\nl2dsZpbVzIwtItZFxKLq963UHh7a/QGqAFdSe3jqi60Ys2dsZpbVqs/hJR1L7bl/j3Y7fiRwYUTM\nlDS9FedysJlZViuWe1SXobcBX6lmbvWuA66pb97s+RxsZpaVmrE99NBD2e/T7lXtB3sb8OOI6OlJ\nz9OAW6utIA8DzpG0MyK672rWZ17u0Qc7dqQ3Zx8xYkQHR9IeTzzxRLJ26qmndnAkA8u1116brF11\n1VUtP9+VV6b3qb7++usb6rMVyz22bNnSp7ajRo3qcbmHpFuo7c9xdR/OdxNwT0Tcsc+DreMZm5ll\nNTP5kXQ6tb0vlkhaSG2vim9Q2/YxImJW99M1fLI6DjYzy2om2KqNsofuQ/svNHyyOg42M8sajF+p\ncrCZWZaDzcyK46d7mFlxBuOMzcs9BqD3v//9ydrChQs7OJLGjR49OlnbvHlzB0eyf2vFco/cHiP1\nxo8fP2Ce7uEZm5llDcYZm4PNzLIcbGZWHAebmRXHwWZmxfFyDzMrjmdsA8CUKVOStaVLlzbU50UX\nXZSs3X777Q31mZNb0jFp0qRkbeXKlS0fS6O8pKMcDjYzK46DzcyK42Azs+I42MysOL4rambF8YzN\nzIrjYNtHuR1uZsyY0VCfjS7pyMkt6Zg9e3ay9slPfjJZGzt2bLK2cePGZG0gLemwno0fPz5Z6+uT\nMgYSB5uZFWcwBtuQ/h6AmQ1sEdGnn55I+qGk9ZKeTNQPkXS3pEWSlkj6fCvG7GAzs6xmgg24Cfh4\npvs/AZZFxFRgJvBX1QbLTfGlqJllNbPcIyLmS5qYawK8s/r9ncDGiNjV8AkrDjYzy2rzZ2w3AHdL\negE4GPjDVnTqYDOzrFSwPf744yxYsKDZ7j8OLIyIMySdAPxS0vsiYmsznfZrsDW6pGMgyS3pyMkt\n6RhInn322WTthBNO6OBI2mPnzp3J2gEHHNBQn4NxSUdOKtimTZvGtGnT3nz9ve99r5HuLwP+Z3We\nZyX9G3AS0FRi+uaBmWU1efMAQNVPT7qAMwEkjQPeDaxqdsy+FDWzrGY+Y5P0U+AjwFhJzwHfAobX\nuo1ZwP8AflS3HORPI2JTcyN2sJlZL5oJtoi4uJf6WvLLQRriYDOzLD/dw8yKMxi/UuVgM7MsB1uH\nnHPOOcnanDlzOjiS9hhIG7aUsKQjZ/r06clablOdTmvHk3D6ysFmZsVxsJlZcRxsZlYcB5uZFcfL\nPcysOJ6xmVlxHGwdUsKSjpxt27a1vM9ly5YlawceeGCydtxxx7V8LA888ECyNnPmzIb6fP3115O1\nd7zjHcnaQFrSkdOfT8JxsJlZcRxsZlYcB5uZFcfBZmbF8XIPMyuOZ2xmVhwHm7XEmjVrWt7nySef\n3PI+G9Xoko6c3JKOnNzSodxTZPYnDjYzK85gDDbvUmVmWc3sUiXph5LW123W0r1+saTF1c98Se9t\nxZgdbGaWtWfPnj79JNxEfrOWVcDvRMQp1Has+kErxuxLUTPLanKXqvmSJmbqj9S9fASY0PDJ6jjY\nzCyrg5+x/UegJV8Ed7CZWVYq2JYsWcLSpUtbcg5JM4HLgJZ829/BNsicdtppydrDDz/cwZG0Rzue\n/JFTwpKOds+oUv1PmTKFKVOmvPn61ltvbah/Se8DZgFnR8TmhjrpxsFmZlktCE5VP28vSMcAtwOf\ni4hnmz3RXg42M8tqJtgk/RT4CDBW0nPAt4DhtW5jFvBfgDHAjZIE7IyI9J6IfeRgM7OsZr4EHxEX\n91L/I+CPGj5BgoPNzLIG4zcPHGxmluVgM7PiONgGgKlTpyZrixYt6uBI2qOEJR057VjSMVgMHz48\nWXvjjTeStdpn7u3jYDOz4jjYzKw4DjYzK473PDCz4njGZmbFcbCZWXEcbD0477zzkrVf/OIXLT9f\nCUs6urq6krWJE5PP7CvCzTffnKxdeumlHRxJe1xyySXJ2i233NJQn+vWrWt0OH3iYDOz4jjYzKw4\nDjYzK46Xe5hZcTxjM7PiONjMrDgOth60Y0lH6Q477LD+HkKf5DZCmTOnsV3USljSkdPoko6cyZMn\nt7zPeg42MyvOYAy2If09ADMb2CKiTz8pks6WtELSryVdk2jzEUkLJS2VlN6DsY88YzOzrGaWe0ga\nAtwAfBR4AXhc0l0RsaKuzaHA3wFnRcQaSU1/FuMZm5llNTljmw6sjIiuiNgJ3Apc0K3NxcDtEbGm\nOt+GZsfsYDOzrCaDbQKwuu7189Wxeu8Gxkh6QNLjkj7X7Jh9KWpmWR24eTAMOBU4AzgIeFjSwxHx\nTDMdDjorVqxI1k466aQOjgQ2btyYrG3dujVZy31ucdxxxzU0lgULFiRr06ZNa6jPnHnz5rW8T9t3\nL7/8clv7TwXbypUreeaZXrNnDXBM3eujqmP1ngc2RMR2YLukfwZOAfavYDOzzkkF24knnsiJJ574\n5uv77ruvp2aPAydKmgisBT4FfLpbm7uA6yUNBUYAHwT+upkxO9jMLKuZu6IRsVvSFcBcap/p/zAi\nlkv6Yq0csyJihaT7gCeB3cCsiHiqmTE72Mwsq9nP2CLiXmByt2Pf7/b6u8B3mzpRHQebmWUNxm8e\nONjMLMvBZmbFcbB1SKeXdOQcfPDBydqXvvSlZG327NktH0s7lnTk7Nixo6Pnu+OOO5K1T3ziEx0c\nyf7FwWZmxXGwmVlxvOeBmRXHMzYzK46DzcyK42Azs+I42HqwcuXKZG3SpEntPn3bjRgxIllrx5KO\nRq1duzZZO+KIIzo4ksZ1eknH6NGjk7XNmzd3cCT9y8FmZsVxsJlZcbzcw8yK4xmbmRXHwWZmxXGw\nmVlxHGw96PSSjsceeyxZmz59egdHMrAMliUdA8lAWtJxww03JGtXXHFFW8/tYDOz4jjYzKw4g3G5\nh3eCN7OsJneCR9LZklZI+rWkaxJt/lbSSkmLJE1tdswONjPLaibYJA0BbgA+DpwMfFrSSd3anAOc\nEBGTgC8C32t2zA42M8tqcsY2HVgZEV0RsRO4FbigW5sLgFuqcz0KHCppXDNjdrCZWVaTwTYBWF33\n+vnqWK7Nmh7a7JO23zwYOXJksrZ9+/aWn29/XtJhA8OwYen/rHbt2pWs3Xjjjclau5d05KRCa/Xq\n1Tz//PMdHk3f+K6omWWlgu2oo47iqKOOevP1o48+2lOzNcAx9W+rjnVvc3QvbfaJL0XNLGvPnj19\n+kl4HDhR0kRJw4FPAXd3a3M3cAmApA8BWyJifTNj9ozNzLKaWaAbEbslXQHMpTaR+mFELJf0xVo5\nZkXEP0k6V9IzwGvAZc2O2cFmZlnNfvMgIu4FJnc79v1ur1v6IaKDzcyy/JUqMyuOg60H7VjSYftu\n06ZNydqYMWM6OBJLufzyy/t7CD1ysJlZcQbjl+AdbGaW5RmbmRXHwWZmxXGwmVlxHGxmVhwH235o\ny5YtydqoUaM6OJK8dizpeOGFF5K1I488suXn67T7778/WTvzzDOTtcmTJydry5Yta2gsZ5xxRrI2\nb968hvrsKwebmRXHyz3MrDiesZlZcRxsZlYcB5uZFcfBZmbFcbDthwbSko5OK2FJx9SpTe/N+zaN\nLunIafeSjhwHm5kVx8s9zKw4g3HG5l2qzCyryQ2TkySNljRX0tOS7pN0aKLdoZJ+Jmm5pGWSPthb\n3w42M8tqV7ABXwfuj4jJwDzgzxLt/gb4p4h4D3AKsLy3jh1sZpbVxmC7ALi5+v1m4MLuDSQdAvx2\nRNxUjWVXRLzSW8cONjPLamOwvWvvxsgRsQ54Vw9tjgM2SLpJ0hOSZkl6R28d++ZBk3bs2JGsjRgx\nooMj6bxTTjklWdu2bVuytnLlynYMpyGLFi3q7yEMeKnQeumll9iwYUP2vZJ+CYyrPwQE8Oc9naqH\nY8OAU4E/iYgFkq6jdgn7rdx5HWxmlpVa7jF27FjGjh375uunn376bW0i4mOpfiWtlzQuItZLGg+8\n2EOz54HVEbGgen0bcE1vY/alqJlltfFS9G7g89XvlwJ39XDu9cBqSe+uDn0UeKq3jj1jM7OsNq5j\n+1/AbElfALqATwJIOgL4QUScX7X7MvCPkg4AVgGX9daxg83MstoVbBGxCXjbo4gjYi1wft3rxcAH\n9qVvB5uZZQ3Gbx442Mwsy8G2H2rHko758+cnazNmzEjWhg4dmqzt3r27qTH1ZPHixS3vs9Muuuii\nZO3222/v4EgGLgebmRXHT/cws+J4xmZmxXGwmVlxHGxmVhwHm5kVx8FmLZFb0pHTjiUdnXbvvfcm\na2effXbLz5db0iEpWVuxYkWyNnny5IbGcvzxxydrq1ataqjPVvBdUTMrjmdsZlYcB5uZFcfBZmbF\ncbCZWXEcbGZWHAdbCx1wwAHJ2s6dOzs4kjIccsghydorr/S6m1lL5Z5C0uiSjlGjRiVrW7ZsaajP\n3H/QjS7pWLNmTUPv609e7mFmxfGMzcyKMxiDzbtUmVlWu3apkvT7kpZK2i3p1ESboyTNk7RM0hJJ\nX+5L356xmVlWG2dsS4DfA76fabMLuDoiFkk6GPhXSXMjIv2dNhxsZtaLNu5S9TSAMl/KjYh1wLrq\n962SlgMTAAebmTVuoHzGJulYYCrwaG9t2x5sP/nJT5K1z372s8mal3S0VqeXdOQ0+hSSYcPS/1xf\nf/31ZK3RJ4bs2LEjWWt0E58JEyY09L7+lFrusXXrVrZu3Zp9r6RfAuPqDwEBfDMi7unrGKrL0NuA\nr0RE/qR4xmZmvUjN2A466CAOOuigN1+vX7++p/d+rNnzSxpGLdR+HBF39eU9DjYzy+rQpWj64Xfw\nD8BTEfE3fe3Myz3MLKuNyz0ulLQa+BDwc0lzquNHSPp59fvpwGeAMyQtlPSEpF6/nuIZm5lltfGu\n6J3AnT0cXwucX/3+KyD9HbwEB5uZZQ2Uu6L7wsFmZlkOth7klnS0Qztu0dvAsGvXroZqV199dbL2\n1FNPJWv+91Ljp3uYWXE8YzOz4jjYzKw4DjYzK46DzcyK42Azs+I42AaAdtyi37hxY7I2duzYhvqc\nOHFistbV1dVQnzlz585N1s4666yWn28gyS3pGDNmTLK2adOmdgxn0PFyDzMrjmdsZlYcB5uZFcfB\nZmbFcbCZWXEcbGZWHAdboRpd0pHTjiUdM2fObHmfpfOSjt55uYeZFcczNjMrzmAMNm/mYmZZbdzM\n5fclLZW0W9KpmXZXVe2elPSPkob31reDzcyy2hVswBLg94AHUw0kHQlcCZwaEe+jdpX5qd469qWo\nmWW1cZeqpwEk5fYUhdouVQdJ2gMcCLzQW98ONjPL6s+7ohHxgqS/Ap4DtgFzI+L+3t7nYCvIAw88\n0N9DsAI1M2OT9EtgXP0hIIBvRsQ9fXj/KOACYCLwMnCbpIsj4qe59znYzCwrFWy7du3K7g5Wvfdj\nTZ7+TGBVRGwCkHQH8FuAg83MGpcKtqFDhzJ06P/fpP2NN95o5jSpz9meAz4kaSSwA/go8Hhvnfmu\nqJlltXG5x4WSVgMfAn4uaU51/AhJP6/O/RhwG7AQWEwtAGf11rdnbGaW1ca7oncCd/ZwfC1wft3r\n7wDf2Ze+HWxmljUYv3ngYDOzLH8J3t5i9+7dyVr9h66D1emnn56svfbaa8naokWL2jEcaxPP2Mys\nOA42MyuOg83MiuNgM7PiONjMrDgONjMrzmBc7qHBmMZm1hmSfkPtyRp90RURx7ZvNH3nYDOz4vhL\n8GZWHAebmRXHwWZmxXGwmVlxHGxmVhwHm5kVx8FmZsVxsJlZcRxsZlYcB5uZFcfBZmbFcbCZWXEc\nbGZWHAebmRXHwWZmxXGwmVlxHGxmVhwHm5kVx8FmZsX5f2OtTQYf7EMuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c11f828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Size of the signal\n",
    "sigSize = (32, 32)\n",
    "\n",
    "# Read and resize the Shepp-Logan phantom\n",
    "f = np.random.binomial(1, 0.1, sigSize) * np.random.randn(sigSize[0], sigSize[1])\n",
    "\n",
    "# Total number of pixels\n",
    "numPixels = np.size(f)\n",
    "\n",
    "# Plot image\n",
    "plt.figure()\n",
    "imgplot = plt.imshow(f, interpolation='nearest')\n",
    "imgplot.set_cmap('gray')\n",
    "plt.colorbar()\n",
    "plt.axis(\"off\")\n",
    "plt.title(\"Sparse Signal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define the measurement model and its transpose\n",
    "\n",
    "# Number of measurements\n",
    "numMeasurements = np.ceil(2*numPixels/3).astype(int)\n",
    "\n",
    "# Measurement matrix\n",
    "H_matrix = np.random.randn(numMeasurements, numPixels)\n",
    "\n",
    "# Define a class object for easier manipulation\n",
    "forwardObj = MatrixOperator(H_matrix, sigSize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Noise standard deviation\n",
    "stdDev = 0.01\n",
    "\n",
    "# Generate noise\n",
    "noise = stdDev*np.random.randn(numMeasurements, 1)\n",
    "\n",
    "# Generate measurements\n",
    "y = forwardObj.apply(f) + noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Benchmark reconstruction\n",
    "\n",
    "# Regularization parameter\n",
    "tau = 1\n",
    "\n",
    "# Total number of iterations\n",
    "numIter = 100\n",
    "\n",
    "# Reconstruct with ISTA\n",
    "[fhatISTA, costISTA] = fistaEst(y, forwardObj, tau, numIter, accelerate=False)\n",
    "\n",
    "# Reconstruct with FISTA \n",
    "[fhatFISTA, costFISTA] = fistaEst(y, forwardObj, tau, numIter, accelerate=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEZCAYAAABvpam5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVPW9//HXZ2kCghRBUapRVIyCGBETy1oSYzT2Bmo0\nptwbb6KJJmrM/YXlJteSGxNN9GqMBgV7bxivdUmIDRsKSlEEKQpILwIL+/n98T3DFLbN7syc2Zn3\n8/E4j5lzZuac73zE/cy3nO/X3B0REZFsVMRdABERaX2UPEREJGtKHiIikjUlDxERyZqSh4iIZE3J\nQ0REsqbkIa2amdWa2W7N/OwhZvZBrsvUhOsONrO3zWyVmf24iZ9p9vcUyQclDykIM5trZuvNbLWZ\nrYke/5SDUzf5RqXMP8DuPtnd985BGbJ1GfCiu+/g7jdmvmhmL5nZBRmHc3ZDlpm1M7MqM5sV/beY\nY2a3mVn/FpxzQBRf/U0pE/oPLYXiwHHu3tXdu0SPF+XgvJZlGYrBAGB6lp/J5ns25mHgeOAsYAdg\nKPAGcFQLzmmE+OaynFLElDykkLb5w2Jm7c1shZkNSTm2Y1RL2THa/4GZzTazz83sMTPrU+fJM36x\nm9l5ZvbP6Pmk6PrvRrWe083scDObn/L+vaJzrDCz98zs2ymvjTOzG83sqejzr5jZoHq/qNkJZjbN\nzJab2Ytmtmd0/AXgCOCm6Dy7Z3zut8ChwI111M6+HtUWlpvZjRmfu8DM3jezZWb29/pqEWZ2NCFJ\nnODub7l7rbuvcfdb3H1c9J4+ZvZ4dK5ZZvb9lM8faGZToia3T83s99FLk6LHlVG5D6ovNlIalDwk\nVu6+ifBLeFTK4TOAanf/3MyOBK4CTgP6AJ8A92Vzieg6h0f7+0a1ngdTXzeztsCTwDNAL+Ai4G4z\n2yPlXGcCY4BuwEfAf9d1QTMbDNwTnaMX8HfgKTNr6+5HAf8E/iMqx4cZ8fjP6PUf11E7Ow44gFBT\nOMPMvhFd70TgCuCk6Hr/BO6tJx5HAa+7+6J6Xge4nxDnnYHTgavMrDJ67QbgenffAfgS8EB0/LDo\nsWtU7tcaOL+UACUPKaTHol/NK6LH70XH7yU9eYwG7k55fru7T3X3GuCXwMEtaJ+vr1nlYKCzu1/r\n7pvd/SXgqYxyPerub7p7bVS+YfWc6wzgKXd/0d23AL8HOgJfbWaZE66OagnzgZdSrv9v0WuzorJd\nAwwzs351nKMn8Gl9FzCzvoRYXO7uNe4+FbgN+E70lhpgdzPr6e7r3f31zFM0+9tJq6LkIYV0orv3\ncPfu0ePt0fGXgI5Rk8gAwi/rR6PXdgHmJU7g7uuAZcCuOS5bH2B+xrF5Gdf5LOX5emD7es6VWWaP\nzt3SMi+u5/oDgBuihLycEB+v53rLCN+1PrsAy919fcqx1DhcAOwJzDCz18zsuOy/hpSCtnEXQMpK\nnb9K3b3WzB4g1DIWE361J/54LSL8cQwnMOtM+PW8oI5TrQM6pezvnEXZFgGZv9T7AzOzOEfqub6c\ncawfdZe5Ltl27M8Hfuvu9TVVpXoeuMjMdqmn6WoR0MPMOkeJGkIcFgK4+0eE/06Y2anAQ2bWoxll\nllZONQ8pFvcS+hRGE/oLUo9/18z2M7MOhP6PV6Omm0zvAKeYWceoI/p7Ga9/BtR3r8RrwHozu8zM\n2kZt/MdTf99BQx4AjjOzI6Jz/RzYALzSxM8vbqCcdbkFuDIx6MDMdjCz0+p6o7u/ADwHPGpmw82s\njZltb2b/Zmbnu/sC4GXgajPrYGb7EeI4ITr32YmBDMAqQtKoBZZGj1/KotzSihV18jCzTmZ2h5n9\nxcxGx10eabEno5E4ie3hxAtR2/k6QpPK31OOvwD8P+ARwq/fQYQhplvfkvL8j4Q2+c+AccBdGdev\nAsZHzTtpf1yj/pRvA98CPgduBM5199l1XKdB7j4LOCc6x1JCR/e33X1zE891A3B6NNrp+no+s3Xf\n3R8j9HPcZ2YrgXeBbzZw/tOApwkd4yuB9wgd8c9Hr48ixHkRYTDD/4v6gIjOO93MVhPifaa7b3T3\nLwgDCP4VxXdEI99RWjkr5sWgzOwcYIW7TzSz+9z9rEY/JCIieVfQmoeZ3W5mi83s3Yzj3zSzGdGY\n8stTXupLshNzS8EKKiIiDSp0s9U44JjUAxamM7gxOr4PMMrM9openk9IIKAhgCIiRaOgycPdJwMr\nMg6PAGa7+7yo3fk+4MTotUeB08zsJsINXCIiUgSKYajurqSPr19ASChEwzUzJ4gTEZGYFUPyaBEz\nK94efxGRIubuze4OKIahugsJNyEl9I2ONZm7l+02ZsyY2MtQTJvioVgoHk3bWiqO5GGkd35PIcyV\nM8DM2hPG8D8RQ7lapblz58ZdhKKieCQpFukUj9wq9FDdewh3rw42s0/M7LseJo77CfAsYY2D+9w9\nq9XdqqqqqK6uznl5RURKTXV1NVVVVS0+T1HfJNgUZuat/Tu0RHV1NZWVlXEXo2goHkmKRTrFI52Z\n4S3o81DyEBEpQy1NHsXQYd5imzdsbvxNJUrNdekUj6RSjMXAgQMxM21ZbAMHDszLf4tWP1QX4IpL\nfsnxZxynKqlIiZs3b15ORgqVE7P0ykV1dXVOfliURLPVuw/NZN9TB8ddFBHJs6ipJe5itCr1xUzN\nVsDSuesaf5OIiORMSSSPJfO+iLsIsSnFdu2WUDySFAvJp5Lo87jz+VvZuXqT+jxERBqRqz6Pkqh5\nfKXn+WWbOMr1e9dH8UhSLApv0KBBvPjii9TU1HDppZfSr18/unbtym677cYll1wCQJcuXejatStd\nu3alTZs2dOrUaeuxe+9Nrnp8xx13UFFRwYMPPpjTMlZWVubkJsGSSB5LPi+JryEiJeLqq6/mrbfe\n4o033mD16tVUV1czfPhwANasWcPq1atZvXo1AwYMYOLEiVuPjRo1aus5xo8fT8+ePRk/fnxcX6NB\nJfFXd+mq9nEXITZq106neCQpFvGZMmUKJ598MjvttBMA/fv355xzztnmffVNUjhv3jz+8Y9/cOut\nt/LMM8+wZMmSvJc5WyWRPJas7RR3EUREtho5ciTXXXcdN998M9OmTcv68+PHj+crX/kKJ598Mnvv\nvTd33313HkrZMiWRPD5Yd3vZ/spSu3Y6xSOpbGNh1vKtha688kquuOIK7rnnHg488ED69u2bVfPT\nhAkTOPvsswEYPXp0TpuucjUxYuxzyudgTnrvZitcREpf+JNVvAYOHOgvvPBC2rENGzb4TTfd5G3a\ntPEZM2Y0+v7Jkyd727ZtffHixe7uPnfuXK+oqPCpU6c2q0z1xSw63uy/vSVR81jrndm0dlPcxYhF\nuda46qN4JCkWxaFDhw5ceOGFdO/enffff7/R9995550ADBs2jD59+jBy5EjMbOvxYlESyWPHiuV8\nPmt53MUQEQHghhtuYNKkSWzYsIEtW7Zw5513snbtWvbff/8GP7dx40YefPBB/vrXv/LOO+8wdepU\npk6dyp/+9CfuvvtuamtrC/QNGlcSyaNX+1Us/XBV3MWIRdm2a9dD8UhSLAovMQlhp06duPTSS+nT\npw+9evXi5ptv5pFHHtlmhtvMSQsfe+wxOnXqxLnnnkvv3r23bhdccAFbtmzhmWeeKdRXaVRJTIx4\nVPc3ufwy5+tXHBB3cUQkjzQxYvY0MWIDPttyE5NefTnuYsRC7drpFI8kxULqkqvRViWRPI4a9B12\ntP3iLoaISNHL1fQkJdFs9ZujXmL9erjq5cq4iyMieaRmq+yp2aoBvXeuYOnykvgqIiKtQkn8xe21\na3uWrOoQdzFioXbtdIpHkmIh+VQS63n0HtiJpWtbPqWAiIg0TUn0ecx69mOO/ZbxYc2AuIsjInmk\nPo/s5avPoyRqHrdPvIVPNx8CKHmIiDQkVysJlkTNo3ZLLR3a1LBmldOha3n1fVRXV+tO4hSKR1Ip\nxkI1j+xptFUDrMLoVbGMpTM1v5WIxGfgwIF06tSJrl27bl1a9pVXXqGiomLrvFQLFy7ktNNOo1ev\nXnTv3p399tuP8ePHM3ny5K2f2X777amoqEg7z4IFC7Ze5/zzz6ddu3YsXrw4rq9aGskDoHeHVSyZ\ntTLuYhRcqf2ybCnFI0mxKDwzY+LEiaxevXrr0rK77LJL2hxW5557LgMGDGD+/PksW7aMCRMmsNNO\nO3HIIYds/cz06dMxM1atWrX1WN++fQFYv349jzzyCN26deOuu+6K66uWTvLo1XEdS+eui7sYIlLm\nGmtWmzJlCueddx7bbbcdFRUVDB06lGOOOabJ53rooYfo3r07v/71r7njjjtyUeRmKZnk0bvrBpbM\n3xB3MQpOY/nTKR5JikVxOvjgg7nwwgu5//77mT9/ftafHz9+PKNHj+bMM89kxowZvP3223koZeNK\nJ3n0qGHJws1xF0NEYhb3KrQnnXQSPXr0oEePHpxyyinbvP7ggw9y2GGH8dvf/pbddtuN4cOH88Yb\nbzTp3J988gkvvfQSo0ePpnfv3hx99NE5XaI2GyWTPHrtCEuXxl2KwlO7djrFI6lcY+He8q0lHn/8\ncZYvX87y5ct55JFHtnl9hx124KqrruK9995j8eLFDB06lJNPPrlJ554wYQJDhgxh3333BWDUqFHc\nc889bNmypWWFboaSSR69+7RhyfI2cRdDRMpcNkOJe/Towc9//nMWLVrEihUrGn3/hAkTmDNnDn36\n9KFPnz5ceumlfP755zz99NMtKXKzlEzy6LVre5aW4fxWatdOp3gkKRbFIzWhXHHFFUyfPp0tW7aw\nZs0a/vd//5fdd9+d7t271/sZgFdeeYU5c+YwZcqUrcvTTp8+nVGjRsWyvnlJ3GFeVVVF700DWLxu\nRNxFEZEylrmsbF3H169fz8knn8xnn31Gx44dOeigg3jiiScaPdf48eM56aSTGDJkSNrxiy++mMMO\nO4yVK1fSrVu3RsuoO8wjZubuzuJpS9lnvwo+r+0Zd5FEJE90h3n2dId5I3oP2ZHN3pZls3WXuYhI\nvpVM8rAKY3DnBcx6aWHcRSkotWunUzySFAvJp5JJHgCDe61k1pRVcRdDRKTklUyfB8BvjnyJDRuN\n//5XZbyFEpG8UJ9H9tTn0QSD9+3ArHnlN1xXRKTQSit5jOzBzGXlNdpK7drpFI8kxULyqSTu80jY\n44i+fLihgtrNtVS0Lam8KCLAgAED6r2XQuo2YEB+VlgtqT4PgL5tPuXlybX0P3jXGEslIlLc1OeR\nYXDXz5g56bO4iyEiUtJKL3n0WcOst9bGXYyCUbt2OsUjSbFIp3jkVskljz0H1zJrVtylEBEpbSWR\nPKqqqrb+qhg8rDOzFnaKt0AFVK5rNtRH8UhSLNIpHkF1dTVVVVUtPk/JdZjPfm4ux3yrgjk1/WMs\nlYhIcVOHeYZBh+zKos292bh6Y9xFKQi146ZTPJIUi3SKR26VXPJo27EdA9sv4qPq7BeWFxGRpim5\nZiuAE3Z+je+e55x87ciYSiUiUtxa2mxVUneYJ+zZ7wtmvRt3KURESlfJNVsBDN67glkftYm7GAWh\ndtx0ikeSYpFO8cit0kweB3Rl5mc7xF0MEZGSVZJ9Hp9NXcx++1ewpLZXTKUSESluGqpbh5327c0G\n78CKj1fGXRQRkZJUksnDKow9Oy9g5rPz4i5K3qkdN53ikaRYpFM8cqskkwfAkL5rmP6CZtcVEcmH\nkuzzAPjdWW/y2Vuf8odZx8dQKhGR4qY+j3rsc+TOTJ/fJe5iiIiUpNJNHl/fhekb94BPP427KHml\ndtx0ikeSYpFO8citkk0e/QcYqyq6sfLFt+IuiohIySnZPg+AEbsu5PojHuerd11Y4FKJiBS3ku3z\nMLNBZnabmT3Q3HPss4/z/pR1uSyWiIhQxMnD3T929++35BxDvtaD6R93gi1bclWsoqN23HSKR5Ji\nkU7xyK28Jw8zu93MFpvZuxnHv2lmM8xslpldno9r73NgJ6a3GwrTp+fj9CIiZSvvfR5mdgiwFhjv\n7vtFxyqAWcBRwCJgCnCWu88ws3OB/YH/cfdPzexBdz+9gfPX2+cxbx58de/lLLzhYfjBD3L7xURE\nWrGi7/Nw98nAiozDI4DZ7j7P3WuA+4ATo/dPcPdLgI1mdjMwrLk1k/79YXXt9qycNLUF30BERDLF\ntRjUrkDqOrELCAllK3dfDvyoKSc7//zzGThwIADdunVj2LBhVFZWYgZ9+7zAhBdm8pPovYl2z8rK\nypLYv/7667d+32IoT9z7ikdyP7WNvxjKE/d+ucejurqaO+64A2Dr38uWKMhQXTMbADyZ0mx1KnCM\nu/8w2j8HGOHuFzXj3PU2WwFccH4tI++9iB8uvQq6dm3eFyhi1dXVW/+hiOKRSrFIp3ikK/pmq3os\nBPqn7PeNjuXcPvtWML3n4TBlSj5OHzv9z5BO8UhSLNIpHrlVqORh0ZYwBdjdzAaYWXvgLOCJ5p68\nqqoqrUqaap99CCOuXn21uacXESkZ1dXVVFVVtfg8hRhtdQ9QCfQEFgNj3H2cmR0LXE9IYLe7+zXN\nPH+DzVbz58OIoRv4dP/j4IUXmnOJoqaqeDrFI0mxSKd4pGtps1XeO8zdfXQ9x/8O/D3f1+/bF9Zt\n7sDy12bTY+1a2H77fF9SRKTkNVrzMLNr3f3yxo7FpbGaB8DIkfA/NT/l0DFHwgknFKhkIiLFqxAd\n5l+v49ixzb1gPjTU5wFRv8eAb8EzzxSuUCIiRSjvfR5m9iPgQmA34KOUl7oA/3L3c1p89RxoSs3j\nj3+Ej15fxo2vfgXmzAFrdrItOmrHTad4JCkW6RSPdPns87iH0CdxNXBFyvE10Q18rcawYfDwwz1g\n82aYNQv23DPuIomItGpN6fP4ErDA3TeaWSWwH2GeqpUFKF+jmlLzWLEiTFWy6swfUvHlIfDTnxao\ndCIixakQfR4PA1vMbHfgVqAfoVbSanTvDjvuCLOHngZ/z/sALxGRkteU5FHr7puBU4A/u/svgD75\nLVZ2GuswBxg+HN7ufAi8/DKsX1+YghVAY9+73CgeSYpFOsUjyFWHeVOSR42ZjQK+AzwVHWvX4ivn\nUFVVVaMdYfvvD2/P7BSyiP4RiUiZqqysLMwd5mY2BPh34BV3v9fMBgFnuPu1Lb56DjSlzwNg4kS4\n/np47qhrYNEi+NOfClA6EZHi1NI+jyZNTxLNPzU42p0ZrcFRFJqaPBYtgv32g6XPvYOdfhrMnl1S\nQ3ZFRLKR9w7zaITVbOAm4H+BWWZ2WHMvGJc+faBtW1jQcyjU1sI778RdpJxQO246xSNJsUineORW\nU/o8rgO+4e6Hu/thwDHAH/NbrOw0pcPcLOr3eMdg1Ci4p1UNGBMRyYmCzaprZu8mFnFq6Fhcmtps\nBXDlldC+PVSdPh2OOQY++QQq4lrSREQkPoW4z+MNM7vNzCqj7TbgjeZeME777w9vv02Y7KpnT/jn\nP+MukohIq9SU5PEj4H3gomibRhPXFi82w4dHyQPg7LNLoulK7bjpFI8kxSKd4pFb9SYPM+tlZkPc\nfaO7/8HdT3H3U4DngFa5GPigQbBqFXz+OXDWWfDww7BpU9zFEhFpdRqqefwZ2LGO4z2AG/JTnOZp\nSoc5hO6NYcOi2kf//jBkCPzf/+W9fPmkWULTKR5JikU6xSMoxJTsb7j7V+p5bZq7f7nFV8+BbDrM\nIcyJuMsucNllwC23wKRJcO+9+SugiEgRymeHeZcGXiuq6Umykdbvcdpp8PTTsGZNrGVqCbXjplM8\nkhSLdIpHbjWUPD40s29lHjSzY4E5+StSfu2/P7z1VrSz445w6KHw0EOxlklEpLVpqNlqD2Ai8DLw\nZnT4K8DBwPHuPqsgJWxEts1WmzeHKdrnz4du3YCnnoKxY2HKlPwVUkSkyOSt2crdZwP7ApOAgdE2\nCdivWBJHc7RtCwccAK+9Fh049lhYtgxefz3WcomItCYN3ucRDdMd5+6XRtvf3H1DoQqXLwcfHJb1\nAKBNG7jwQrjxxljL1Fxqx02neCQpFukUj9wqy7k5vvpVeOWVlAMXXABPPglLlsRWJhGR1qRJU7IX\nMzPzMWPGUFlZ2eRx3EuXwh57wPLlKVNbfe978KUvhQmwRERKVHV1NdXV1YwdO7Yg63l0BPq7+8zm\nXihfsu0wT9hjD3j0Ufhy4m6Vt96CE0+Ejz8OHSMiIiWsEOt5fBt4B3gm2h9mZk8094LFIq3fA8IN\nIP37h+arVkTtuOkUjyTFIp3ikVtN6fOoAkYAKwHc/R1gUB7LVBDb9HsA/PjH8MeiWqpERKQoNWU9\nj1fdfaSZve3u+0fHWuV6HqmmToUzz4QZM1IObt4Me+8Nt90Ghx+eu0KKiBSZQqznMd3MRgNtzGwP\nM/sz4cbBVu3LXw7rmi9blnKwbVv41a/CTYMiIlKvpiSPnwD7ABuBe4BVwE/zWahCaNMGRoyAV1/N\neOHss2Hu3FazUJTacdMpHkmKRTrFI7eakjz2cvdfufuB0fafpXCjIIRO8236Pdq1U+1DRKQRTenz\neAnYGXgIuN/dpxWiYE3V3D4PCBPq/v738OKLGS/U1MDgwXDXXfC1r7W8kCIiRSbvfR7ufgRwBLAU\n+IuZvWdm/9ncC+ZDUxeDyjRyJLzxRugnT9OuXbhZ8L/+KyflExEpFnlfDKrON5vtC1wGnOnu7Vt8\n9RxoSc0DwuCqe+8NKwym2bQJ9twTxo2DIl6BrLq6WiukpVA8khSLdIpHukLcJLi3mVWZ2XuEpWlf\nBvo294LF5pBD4B//qOOF9u3h2mvD0oNbthS8XCIixawpfR6vAPcDD7j7ooKUKgstrXncd1+oeTz+\neB0vusNhh8F558H3v9/8QoqIFJmW1jxKYmLElnyHJUtC3/jnn9czpdWbb8Lxx8PMmdC1a/MLKiJS\nRPLWbGVmD0SP75nZuynbe2b2bnMvWGx694aBAxtYSPCAA+Cb34T//u9CFqvJNHY9neKRpFikUzxy\nq6HpYy+OHo8vREHidNRR8Pzz4b6POl11Vbgl/Yc/DNO2i4iUuab0eVzr7pc3diwuLW22Avj73+Ga\na2DSpAbedM018NJL8MwzYM2u6YmIFIVCzG319TqOHdvcCxajQw8NXRvr1jXwpksvDR0kd95ZsHKJ\niBSrhvo8fhQNz90zo8/jY6Bk+jwAtt8+dG00OJ1Vu3bhno/LLoNPPy1Y2Rqjdtx0ikeSYpFO8cit\nhmoe9wDfBp6IHhPbAe5+TgHKVlCJfo8GDRsW+j1+9KMwjFdEpEw1eaiumfUGtkvsu/sn+SpUNnLR\n5wFhVcELL4R33mnkjRs3wv77w5gxYUEQEZFWKO/3eUTL0P4B2AVYAgwAPnD3fZp70VzKVfKoqYEd\nd4QPP4RevRp582uvhfXO33gD+pbMzfYiUkYK0WH+W2AkMMvdBwFHAZmrYMSquRMjpmrXLiweuM0M\nu3U56CD4yU9g9Og6ZlUsLLXjplM8khSLdIpHkKuJEZuSPGrcfRlQYWYV7v4S8JUWXzmHqqqqcjLh\n2dFHw3PPNfHNV1wBHTpo3Q8RaVUqKysLM6uumT0PnARcDexIaLo60N2/2uKr50Cumq0AZs8OU1kt\nXAgVTUmrixfD8OFh+O7RR+ekDCIihVCIZqsTgS+AnwHPAB8RRl2VnD32CP0e2yxNW5+ddoLx4+E7\n3wkLoouIlImmLAa1zt23uPtmd7/T3f8UNWOVpFNOgYcfzuIDRx0F//EfcNJJsH593spVH7XjplM8\nkhSLdIpHbjV0k+AaM1udsq1JfSxkIQvplFPgkUeyvI3jyivD1Lznnw+1tfkqmohI0Sj7KdkzuYfm\nqwcfDLdzNNmGDXDkkaHvQ8vXikiRK0SfB2Z2iJl9N3q+o5kNau4Fi51ZsvaRle22g8cegwkT4K67\n8lI2EZFi0ZRlaMcAlwO/jA61B0r6r2OzkgeExUGeegp+/vPwWABqx02neCQpFukUj9xqSs3jZOAE\nYB1AtBRtl3wWKm4jRsDKlTBjRjM+vM8+8MQTcMEFYQp3EZES1JT7PF539xFm9pa7DzezzsAr7r5f\nYYrYsFz3eST85CfQp0/oC2+WSZPg9NNDIhk5MqdlExFpqUL0eTxgZn8BupnZD4Dngduae8HWotlN\nVwmHHx5uHjzxxDAXlohICWnKfR6/Bx4CHgb2BH7t7n/Kd8Hiduih8MknYaLEZjv2WPjb3+Db385b\nE5bacdMpHkmKRTrFI7eaNNrK3Z9z91+4+8+BF8zs7DyXK3Zt24Ybx29raR3ruOPCuN8zzyxYJ7qI\nSL7V2+dhZl2B/wB2JSwI9Vy0/3NgqrufWKhCNiRffR4AM2eG1qdPPoH27Vt4stdfhxNOgN/9LmQl\nEZEY5W09DzN7HFgBvEKYhr03YMDF7t7YkkkFk8/kAXDEEWGRqNNPz8HJ3n8fjj8eRo2C3/ymibMv\niojkXj47zHdz9/Pd/S/AKGAIcEwxJY5C+Ld/g7/8JUcnGzIkdJ5PmgRnnJGTubDUjptO8UhSLNIp\nHrnVUPKoSTxx9y3AAnffkP8iBWZ2opndamb3mtnXC3XdTCefDO++28KO81S9esELL0CnTqFXfs6c\nHJ1YRKRwGmq22kJ0YyChuaojsD567u7etSAFNOsG/I+7/6Ce1/PabAXwi1+EFqZrr83hSd3hz3+G\n3/4Wbr01zMorIlIgeV/DvKXM7HbgeGBx6o2FZvZN4HpC7ed2d6/zT7OZ/R64q77mskIkj1mzQiVh\n/vwcdJxneu21MBLr1FPh6qvzcAERkW0VZGLEFhoHHJN6wMwqgBuj4/sAo8xsr+i1c83sD2a2i5ld\nAzwddz/L4MFh1pEW3TRYn4MOgrfeCu1iI0aENrIsqB03neKRpFikUzxyK+/Jw90nE0ZtpRoBzHb3\nee5eA9xHWLEQd5/g7pcApxJGeZ1mZj/Mdzkbc8klcNVVeVquo0ePMCPvT38aFpe6+mrYvDkPFxIR\nyY2CrOdhZgOAJxPNVmZ2KmHk1g+j/XOAEe5+UTPO7eeddx4DBw4EoFu3bgwbNozKykog+WujpfuH\nH17JyJGC0oIUAAAQ30lEQVRwzDHVHHlky89X7/7998PvfkelGdxyC9Vr1+b2/NrXvvbLcr+6upo7\n7rgDgIEDBzJ27Nji7vOA/CePQi1o9eyzcNFFMG1auAM9b9zDuiCXXQannRY61bt1y+MFRaTctIY+\nj7osBPqn7PeNjhW1r389jLS95548X8gs3IX+/vtQUwN77x1GZNXRlJX4ZSGB4pGkWKRTPHKrUMnD\noi1hCrC7mQ0ws/bAWYQpUIqaWagEjB0b/qbnXY8e4Q7Fp54KGWvoUHj66SwXWBcRyb1CDNW9B6gE\negKLgTHuPs7MjiV9qO41zTy/jxkzhsrKyq3tfPl29NFhdO0P6rzzJE/cQxL5xS/CioVjx4a5U0RE\nslBdXU11dXXr6PPIp0L2eSS88kqY62ratBi6IjZvDrWQ3/wGdt0VxoyByspQLRIRaaLW2ufRqh18\ncFii45JLYrh4Yq74Dz6A88+n+jvfCfeHPPCAhveidu1UikU6xSO3lDya6Xe/g+pqmDgxpgK0bQvn\nnx9WK/zVr+CGG8LdjNddB8uXx1QoESkXJdFsVeg+j4TqajjnnHBTeI8eBb103V59FW66KfSNnHJK\nmBL4wAPVpCUiW6nPIxJHn0eqiy6CFSvCbRlFY8kSuP32sHXoAN/9Lpx7Luy0U9wlE5EioT6PmF19\ndZjb8M4747l+ne24vXvDL38Js2fDzTfD9Omw555wzDFwxx2walWhi1kwatdOUizSKR65peTRQp07\nw+OPw+WXhzvQi4oZHHYYjBsHCxfCBReEObT69w9L4o4bB8uWxV1KEWmFSqLZKq4+j1STJ4duhmef\nhWHDYitG06xcGW42fOQReO45GD4cjjsubHvtpT4SkRKmPo9I3H0eqR5+GC6+GP71LxgwIO7SNNH6\n9fDii2HY2MSJYRTXN74RtiOP1JxaIiVKfR5F5NRT4Yorwj17H3xQmGu2uB23Uyc4/vjQNzJvXmiD\nGzwY/vrX0Lw1YkSYoPHpp2H16pyUOZ/Urp2kWKRTPHIrn3PDlqUf/xi6dAkJ5IEH4PDD4y5RFsxg\n333DdsklsHEjvP56GJN83XVwxhmw++7wta+FbeRIGDRIzVwiZUjNVnnywgswahT88Y9w9tlxlyZH\nNm2Ct98O7XL/+le4r2TTplA7GTEi9J0ccADsskvcJRWRRhT9Gub5Viwd5nWZNi0Majr66PDDvUuX\nuEuUBwsXhrHKU6aE5XTffDP0mwwdGkYODB0aajJ77qn12UWKgDrMI8Va80hYvRp+9jN46aVwL8ih\nh+b2/NXV1cWVNN1h/nyYOjW5TZsGc+fCbrvBkCFhfZLEtsceYbxzjhRdPGKkWKRTPNK1tOahPo88\n69o13Oj95JNhGvcTTwyzqffuHXfJ8sQsdLT37x9mj0zYuBFmzgw3LH7wQRgm/MEH8NFH0LNn6KTf\nY4/Qp7L77vClL4X+lJKsrom0fqp5FNCyZWExqQkTQm3kZz8Lg53K2pYtoaYyc2ZIJB9+mHycOxc6\ndgxJZODAMP45sfXrF7aePdVhL9IM6vNoRckj4aOPwuwh//xnGJ317/8e/gZKBvcwT9ecOWEY8bx5\nIaF88klIOPPnhxrNrrtC377hcdddoU+f0Gnfp0/Ydt4Ztt8+7m8jUlSUPFph8kiYNg3+8Ad49FE4\n66ywMuH++2f3Q7rs23HXrAmd9gsXwoIFVE+eTGXHjrBoUdgWL4ZPP4WKijAxZGLr3TtsvXqlbzvu\nGDL5dtvF/c1arOz/bWRQPNKpzwOoqqoqytFWjfnyl+Fvf4OrroJbboHTTgutNOecE/pHdtst7hK2\nAl26hClV9tor7A8YEG6ySeUekszixclt6dJQq5k9Oww7XroUPv88PC5bFkaGJRJJjx7JrXv35GP3\n7uEO/MTjDjuErU2bgodBpKkSo61aSjWPIuIOL78c+kQeeyz8TTruODj22LB6Ydn3jxRKItksWxa2\n5cvDtmxZmH8/dVu5Mvl81arwuc6dw0iJRDLZYYewn9i6dEk+Jrbtt9/2sVMnJSLJGzVblVDySFVb\nG26bmDgRnnkmLDi1775hqO+BB4b78XbbLbTGSBGprQ0JZNWq5LZ6ddgSx9esSR5buzbsJ7Z165LH\n1q8PzWedO4dt++2Tzzt12vYxc+vYMfmYuW23XfJ5hw4adFCGlDxKNHlkWr8+zBQyeTK88Ua40Xvl\nSujXr5qDDqpkzz3DfXiJgUldu8Zd4niUVLt2bS188UVIJuvWhX8E69Ylt8T++vXJbd268Jn166me\nO5fKLl227vPFF+nbhg3hsaYmJJDttksmlcTz1OMdOiS3zP3UrX37bZ+nPjZla9cu5wmtpP5t5ID6\nPMpEp06hKT/13/6yZaGJq1MnmDEDJk0Kg5Hmzg3/b/frFwYdJbZEH3Hv3qEpP7F16BDPd5JGVFQk\naxrNUV29bf9PXWprw6i1RFLZuDEklkRy2bgxeSz1MXNbuzb5fNOmsNW1n3heUxP2a2qSryWe19SE\nmQratUsmk8SWup/5WmJLfDZ1f8mSMDqlbdvk63U9Zj5PbG3abHusodcSx+p6zHzeCpsQVPMoQe6h\n73fBgjAIKTHwaMmSsC1enGzCX748/H+S2eebaLKvq3k+tWk+sXXtGhKWWj8kJ9xDAklsiYTS0PPE\ntnlz/fubNyf36zteUxPuP0rsZz6v67XNm9OfZ743832Zr5mlJ5PM5FLX1tjrmVtFRdq+3X+/mq1a\n+3eIk3v4wbhqVbLvN9Ecn2ieT2yJZvv6ts2b0/uFE1tqMurWreGtY0clICkz7qH2V19iaWhryntS\nt9rarc9t9Gg1W7XWobq50NJ2XLNk7aFv35aVZdOm9CSTSECZ24IFyUSVmrRWrgz/rhtLMIktNRkl\nnk+ZUs0RRzQ/HqVEbfzpijYeZskaQQHkaqhuySQPiV/79sl+lObauDEklEQyyUwyK1eGG8tT91Of\nb9jQeIJJbZ7LPN61a2gNEClViR/aY8eObdF51GwlJaWmpv7Ekqj5JJJRXe9bvToMQKgv8dT3mJqQ\n1PQmrYGG6ip5SA7V1oY+oNSE0lgyytzPpukt8yb17t01+k0KQ8mjzJNH0bbjxqQY4rFhw7a1nMwb\n0VP7eVKb51asCE3fTU00ma9165ZsOi+GWBQTxSOd7vMQKTKJe+p22in7z7qH5JPZ55OaZBYvDjPY\npyajxHtWrw63hXTvHvpu+vVLTsPV2NatWxi2LdIUqnmIlJDa2pBAUmsyy5dvOx1XXcdWrgxJr6Hk\nkjonZOamVYZbFzVbKXmI5ERiPsi6kktj28qVoa+mvsSSmXRSJyju1k0j3OKg5FHmyUPtuOkUj6RC\nxiJxs2ldiaWuRJQ4tnx5aHZLNLU1lHDqqw019fYI/dtIpz4PyvsmQZFikHqzaf/+2X02MRFxakLJ\nTDgff1x3Ylq9OkyV05Q+nfnz05NUNomnlGg9j0i51zxEytmWLSGB1LfESkM1oNTBBZnrfNX1mLp1\n7tz67+VRs5WSh4g0Q21taDJLTSqZz+uqDS1fHm5GTSSSzMUmM48nZl0otqSj5FHmyUPtuOkUjyTF\nIl0u47FhQ3oyyVxwMrHoZOZilJs3pyeb1ATT0GPHjjkpdhr1eYiIFNh22yXXycnGF18kk0pq0kkk\nmQ8/TE9CiceKim1rMalr8tS1de+e3z4d1TxERIqYe1gIMjXJpNZoMo8ltlWrwkSf9SWaX/9aNQ8R\nkZJlllxQsl+/pn9uy5YweKC+ZNPicrX2X+3lXvNQu3Y6xSNJsUineKRraZ9H61s4V0REYqeah4hI\nGVLNQ0RECk7Jo5XLxTQDpUTxSFIs0ikeuaXkISIiWSuJPo8xY8ZoYkQRkSZITIw4duxYTU/S2r+D\niEihqcO8zKkdN53ikaRYpFM8ckvJQ0REsqZmKxGRMqRmKxERKTglj1ZO7bjpFI8kxSKd4pFbSh4i\nIpI19XmIiJQh9XmIiEjBKXm0cmrHTad4JCkW6RSP3FLyEBGRrKnPQ0SkDKnPQ0RECk7Jo5VTO246\nxSNJsUineOSWkoeIiGRNfR4iImVIfR4iIlJwRZs8zGwvM7vZzB4ws3+PuzzFSu246RSPJMUineKR\nW0WbPNx9hrv/CDgT+Grc5SlW77zzTtxFKCqKR5JikU7xyK28Jw8zu93MFpvZuxnHv2lmM8xslpld\nXs9nvw08BTyd73K2VitXroy7CEVF8UhSLNIpHrlViJrHOOCY1ANmVgHcGB3fBxhlZntFr51rZn8w\nsz7u/qS7HwecU4ByiohIE7XN9wXcfbKZDcg4PAKY7e7zAMzsPuBEYIa7TwAmmNnhZnYF0AGYmO9y\ntlZz586NuwhFRfFIUizSKR65VZChulHyeNLd94v2TwWOcfcfRvvnACPc/aJmnFvjdEVEmqElQ3Xz\nXvPIt5Z8eRERaZ64RlstBPqn7PeNjomISCtQqORh0ZYwBdjdzAaYWXvgLOCJApVFRERaqBBDde8B\nXgYGm9knZvZdd98C/AR4FpgO3OfuH+S7LK2dmfU1sxfNbLqZvWdmF0XHu5vZs2Y208z+z8x2iLus\nhWJmFWb2lpk9Ee2Xcyx2MLMHzeyD6N/IQeUaDzP7mZlNM7N3zexuM2tfTrGo6xaJhr6/mf3SzGZH\n/3a+0aRraF6o1sPMdgZ2dvd3zGx74E3CKLXvAsvc/XfRPTPd3f2KOMtaKGb2M+AAoKu7n2Bm11K+\nsbgDmOTu48ysLdAZuJIyi4eZ7QJMBvZy901mdj/hXrEhlEkszOwQYC0wPmWgUp3/b5jZEOBu4EBC\nF8LzwB6NTRpYtHeYy7bc/TN3fyd6vhb4gPAf+0TgzuhtdwInxVPCwjKzvsC3gNtSDpdrLLoCh7r7\nOAB33+zuqyjTeABtgM5REu1I6FMtm1i4+2RgRcbh+r7/CYTWn83uPheYTbidokFKHq2UmQ0EhgGv\nAju5+2IICQboHV/JCuqPwC+A1F9I5RqLQcDnZjYuasa71cw6UYbxcPdFwHXAJ4Skscrdn6cMY5Gh\ndz3ff1dgfsr7FkbHGqTk0QpFTVYPARdHNZDM6mXJt0Wa2XHA4qgm1tBw7ZKPRaQtMBy4yd2HA+uA\nKyjPfxvdCL+yBwC7EGogZ1OGsWhEi76/kkcrE1XDHwImuPvj0eHFZrZT9PrOwJK4yldAXwNOMLM5\nwL3AkWY2AfisDGMBsACY7+5vRPsPE5JJOf7bOBqY4+7Lo8E5jxImVy3HWKSq7/svBPqlvK9Jt04o\nebQ+fwPed/cbUo49AZwfPT8PeDzzQ6XG3a909/7uvhthqPeL7n4u8CRlFguAqDlivpkNjg4dRRjJ\nWHb/NgjNVSPNbDszM0Is3qf8YpF5i0R93/8J4KxoRNogYHfg9UZPrtFWrYeZfQ34B/AeocrphNE0\nrwMPEH49zAPOcPeymULUzA4HLo1GW/WgTGNhZkMJgwfaAXMIo/DaUIbxMLMxhB8VNcDbwPeBLpRJ\nLKJbJCqBnsBiYAzwGPAgdXx/M/sl8D1CvC5292cbvYaSh4iIZEvNViIikjUlDxERyZqSh4iIZE3J\nQ0REsqbkISIiWVPyEBGRrCl5iORBND36j+Iuh0i+KHmI5Ed34MK4CyGSL0oeIvlxNbBbNMPttXEX\nRiTXdIe5SB6Y2QDgycRCPCKlRjUPERHJmpKHiIhkTclDJD/WEGZxFSlJSh4ieeDuy4F/mdm76jCX\nUqQOcxERyZpqHiIikjUlDxERyZqSh4iIZE3JQ0REsqbkISIiWVPyEBGRrCl5iIhI1pQ8REQka/8f\nXxe/hqfz1gcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e55cbe0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAADDCAYAAABJYEAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHHpJREFUeJztnXmMVdWWxtcWJ6CKwSqgoKCYZVIRUZEIxRQRRUSNxIgt\nT6Mt2iFMTRw6arQ1rZH3BI2+jv+8aGxBoyioaBSUoBgwLUopPECKoQARGYpRwfH0HxQvddb6qu/y\nUreum/f9EvPeXu6zz7777LNyXF+ttUOSJEIIISQOTsn3BAghhPih0yaEkIig0yaEkIig0yaEkIig\n0yaEkIig0yaEkIig0yaEkIig084xIYTNIYThIYTTQgh/CSFsCyEcDCFsCiE8WdPnUI3tYAjh1xDC\nD7VsN9Ya65YQwm8hhHH5+0Xkn5UQwpaavXmw1v4cWLMnT6npUxpCeC2EsDuEsC+E8GUIYUIIYVCt\naw7XXFN7nPa17vN8COHnEEKb/P3aPy502g3HfSJygYhcmCRJMxEZKiKfi4gkSVKYJEmzGnuViIyu\nZZtba4wJIrK35n8JaWgSObY3mx3fnyKyo8Z+nBfl2B7uICJFInKziHyXJMmyWtf0qbmmea19vl1E\nJITQRESuE5H9IvIvDfbLIoJOu+G4SETeSJLkOxGRJEm2JknyP6BfqPknbQyho4iUi8gdIjIqhNA6\nl5MlpA7M3lRcJCIvJElyNEmS35IkqUiS5L3fMdb1IrJPRP5TRG7JfponL3TaDccKEfn3EMJdIYRz\nsrh+goh8liTJGyKyVkRuqtfZEVI/LBeRv4YQbgghdMji+gkiMkdEXhGRniGEfvU6u5MAOu2G479E\n5HERGS8i/xtC2B5C+D1hjptF5KWa/z9HGCIh+WF+CKG65p/Xwb8fJyIficj9IrIphPB5COFCz8Ah\nhDIRGSYic5Ik2SUii4X73ECn3UAkx/jvJEkGi0gLOebE/xZC6JHp2hDCpSLSWY59fYiIzBWR80II\n5+VswoRgxiZJclbNP9fpf5kkyYEkSf4jSZJzRaSNiFSIyBvOsW8Wkb8nSfJVTXuuiIwPITSql5mf\nJNBp54EkSX5MkuSvcix219txyZ9q/ndVCOFbORZqSWrZCWkoMsW0/0GSJNUi8mcRaRdCaOm45GYR\n6RJC+LZmn/9FRIpF5MqsZnqSQqfdQIQQpoQQhoQQzgwhNAoh/ElECkTkiwzXnSHH/pPzX0XkfBHp\nW/PPZBG56fifWhGSR/7hyEMIj4cQ+tTs8UIR+TcRqUySZF9d19RcN1BEusgxIfP4Hu8jx762+XFS\ni1PzPYF/Ao7/OdQPcuzLoWuN7WsRuS5Jki119D/ONTXXvpgkya/HjSGEv4nIwyIySkTeqf9pE2Ko\nq/h+bXsTORYOKRGRIyLyqYhc7RhrgojMT5Lk77WNIYSnROSjEEKLJEn2ZzXrk4zAQxAIISQe+J/W\nhBASEXTahBASEXTahBASEXTahBASETn/65EpU6YYpbO0tNT0a9MmXdBrz549ps9FF11kbEuXLjW2\nH374IdVu1Mj+bf5vv/1mbLt37za27t27G9v+/WkRu6CgwDX+6aefbmyHDx9OtT///HPTZ/DgwcbW\nt29fY5s3b56x9erVK9XWa1OXrXnz5sZWVVVlbBUVFan2kCFDTJ+ff/7Z2Jo0aWJsjz/+eKp9zz33\nmD5o70ycONH9t8P1ydSpU83ePuUU+x3Utm3bVHvHjh2mz7XXXmts8+fPNzbPHw6EYJfj119/NTb0\nXujxve8OQo+1du1a06d3b5umUF5ebmxLliwxNj037x9VoPVB1y5btizVHjRokGt8xOzZs1PtKVOm\nmD7t2rUztrvvvttMll/ahBASEXTahBASEXTahBASETlPrnniiSfMDfbu3Wv66TgZisGdeeaZxoZi\nbt9//32qXVxcbProWDKagwiOdel+p55qpQEUo/z222+NbdWqVRnnheK/6J5IB/jpp59SbRR/R/Pq\n3LmzsZ1xxhnGpqmurjY2FB/XzwiBfiOKj8+cOTMvMe3p06ebzYH2i7ahPigWjvrpeGxRUZHpg94v\nz1hexowZY2woTr9y5cqsxvfG5DVoDX/55RdjQ/vKsxbeNfT0Q74GjTVr1izGtAkhJGbotAkhJCLo\ntAkhJCLotAkhJCJynlyDEjeQoLVz585Uu6yszPRBAgsK3mtBorKy0vQpKSnJaiwRkX370qWBdWKQ\niMjLL79sbK1b27N433zzzVR7+PDhrnkhIePQoUPG1qVLl1QbCTPnnGOPrESCKOL+++9PtadPn276\noN+NRNM77rgj1X7ggQdMn7POOss1r4bAK+JrEQ0JYQiPOLZ9+3ZjQ4K9F/2b0BzeeustY0P9VqxY\nkWoPGDDAdZ13XT3JNdmKjiIiTz75ZKo9bdo013Vo/EmTJqXaTz31lOmDfA2CX9qEEBIRdNqEEBIR\ndNqEEBIROU+uue+++8wNUDz2yJEjqTZKmkGJBLt27TI2nUzTrFkz00cXfRLBMSWdnCIi0qdPn1Qb\nJRGg+XuSBlAiCopVo+IySD/Ytm1bqt2iRQvT5+jRo8aG4qKon94/KAlCawAiIueee66x6cQctF5o\nXR999NG8JNdMmzbN9fLo33EiSRoatB7oGXjv2a9fv1T7iy/sEabZzt+bQIT26IEDB4xNJ1qh+LW3\nUJYHNBbSiNC743mWqM/s2bOZXEMIITFDp00IIRFBp00IIRFBp00IIRGR8+QaLTCKYKGhZcuWqTYS\nADt16mRsOilHxAqPl156qemDTrxB90RJJvoUja5du5o+3333nbEhoUTPFc0BiY5IAEE2fS06faZ9\n+/bGhoQfz6ksW7dudY3/448/GpsWqJFgjX7jyQASZpEIqPfQ2LFjTZ8FCxYYGxLRkHioRXX0zqFn\n7BEnvQIm2nvo2tNOOy3VRhUgTySRSe93tB9RomC2oqM36Ydf2oQQEhF02oQQEhF02oQQEhF02oQQ\nEhE5FyLRUVmoSthVV12VaqOg/IYNG4zt7LPPNjZ9fNbChQtNn5kzZxrbbbfdZmy9evUyNp1liI7Y\nQgKsR6BAQlurVq2MDVU8REKhFkq04CuCRUFUTQ+JmKeffnqqjZ73xRdfbGwff/yxsXmyQ5FQmy/Q\nHkUCd3l5ecbrKioqjA2JXPr3v/7666bPnDlzjG38+PHGhvaC3stbtmwxfbI9YgvhPVIOoe95IqIj\neu905iR6vzp27Ghs6D3RIqY3OxTBL21CCIkIOm1CCIkIOm1CCIkIOm1CCImIvBw3dsUVVxibFliQ\nWNC4cWNjQ6KFLmWKrps8ebKxoXuuW7fO2HQWI8rwRKICKgmps7qWLVvmug6Vm0ViihYZ0VjIhkre\nInFS3xOVpURriARFbUOiZmlpqbHlC/SMtegoYgVWb9le9Dy1gIWuQ6IjAh35psdHc81WdFy/fr2x\n9ejRwzWW5ygulPWJQJmNnvHRvHbs2OEaX+PJ8KwLfmkTQkhE0GkTQkhE0GkTQkhE5DymjSrUbd++\n3dh07GzcuHGmD0pcQFX+dIwZxZhQJbGSkhJjQ1X+OnTokGrrZBsRnLCCYm46lta/f3/TByUNoMQL\nFAPV90Tzatq0qbHppBkRkS5duhibTnhC6+WpPihij4BDVe5Q9cR8kW2SxsiRI02f999/P6s5eOKn\nIngPobiqrli5cePGrOaFQPFrb2zXm3Cj8VbO00cIioisXr0643VoXp7EmU8//dT0QUloCH5pE0JI\nRNBpE0JIRNBpE0JIRNBpE0JIRORciESJJ0i400kZb7/9tumDRAsk7mmRsaCgwPQpLi52jYWSC7RI\n9/3335s+SCjctGmTsemKhIMGDTJ9kNiEjnxCCTGXXHLJ/3s/Eb94WFlZaWxa6EHJNajiIbqnfk4X\nXHCB6VNYWGhs+QIJZp4Ej0WLFtXbHLxCm/e4MU+SG3oP165da2x6fbp165bVHESwcDdw4MCM80Lv\nJkraWrNmjWtuHtBaayFYv5e/53780iaEkIig0yaEkIig0yaEkIig0yaEkIjIuRCJKokhtPCFgvmo\nyhwK3mvxEwmRR48eNbaDBw8aG7pWVytDxzah6l/oaDR9lBgS2tBaILEGiau7d+/OOBaytW3b1tie\nffZZY7vzzjtTbSQwtmnTxtiQaKozyVD2JnpG+QI9A8+RUV7B23v8lAa9E55sXBH8RwKesdCxfJ55\neSsGomxBz/og0RGN7zkC0fs8sn2WPG6MEEJOQui0CSEkIui0CSEkIui0CSEkInIuRCJhCgmKuiwn\nOqYMZduh4L0uNYrKeaKSih07djQ2JMxoMezyyy83fd577z3XXIuKilJtlMGFyqQi0PosX7481e7c\nubPpg4TU6upqY7vpppuMTa8Fmisqn4uOgNPiMHpGSED+I+EtgZotWuRCoiC6H9obSEjV4w0bNsz0\nQSWSPYJifa6DiMiSJUtS7cGDB5s+SGxFNi06nggeIdgryiL4pU0IIRFBp00IIRFBp00IIRFBp00I\nIRGRcyHSW75SZy6hoDwq1Tlv3ryM46Myo0gMRedGouwmXSp14cKFGecgIlJaWmpsOmMRCRRInETj\no9KXWhhDa4GyMNH6I+HKU04SiZNorVu3bp1qozKy+nzOPxpo3fR+QUJY3759jQ2VI0UlfzXomWR7\nxiISHdH4SIDV+wX9boT33EgtknqzK09EBPSMhdDPDb1L3vXhlzYhhEQEnTYhhEQEnTYhhEREzmPa\nQ4cONTYUJ9OxY50gIyLywQcfGFv37t2NTSf0oMQCFCdGcTNdhU9EZMOGDal2WVmZ6YNix3v37jU2\nXckOxXpR9T6UnLJx40Zje+GFF1LtGTNmmD7e+BqqmKaTXVD8GsVTJ02aZGwvvvhiqo32AEr6yRej\nR482NnRMnt7b6Pi4lStXGhuKX+vYK4rFemOjSK+pqqpKtTt16pRxDiJ4v+t5oD4oFu6tUvj888+n\n2rfccovpg8h1/HrChAnGpt9DNNa+fftc4/NLmxBCIoJOmxBCIoJOmxBCIoJOmxBCIiLUd+Utzb33\n3mtugIQSXdWvXbt2ps/mzZuNrXnz5samx0eiIxJ5kCiCRDSdmKMTZOoaHwlrem5IcBkwYICx6ep9\nIlic1OOhtUeCFAId/6XHR6Ivuif6nTq5ZsuWLaYPSgR6+OGHfQpRPTN16lSztz1ilTcJxIO7Mhx4\nBkgQ1WKhVyjMNmFl+PDhxvbhhx9mvM5LfSbSILzPUgv0SPxHRxs+8sgjZjB+aRNCSETQaRNCSETQ\naRNCSETQaRNCSETkPCMSZUEhwUxn2y1evNj06dmzp7GhDEItYiJRcNOmTcaGBAR9nJYIFuQ069at\nM7b+/fsbm16fI0eOmD4VFRXGhgRStNb6GDd0tBj63ehYr7Zt2xqbzuJCwhUSi9ERcLqqnydjL59k\nK8ihvdGjR4+M14nYNUH327Nnj7Gh54n2kD6ODs3hs88+M7YLL7zQ2DRorI8++ijjdSJYNNWcSBVB\nz7XoOm8GMBIeNfpdrQt+aRNCSETQaRNCSETQaRNCSETQaRNCSETkXIhEWYBIyNOCHzp+CWU26iw6\nNL7OthTB4uTtt99ubHPnzjU2LT4gEeP88883tkOHDhmbFpLQerVo0cLYkEiHSsvquaHytiNGjDA2\nlNmIRBctsKDykij7EQmWenx0P2/2Zr7wHDfWq1cv13VeYU2DhP677rrL2F599dWMYyGhE4mOnszA\nE8lORHtbC+/z5883fa655pqM86oLTxlcj8Aoguef6X51wS9tQgiJCDptQgiJCDptQgiJCDptQgiJ\niJyXZh0/fry5QdeuXU0/LR56M/KQqKCzCpEQhsQ9lEmGRFMt0ulSrSI4K8qTSYmEDSS+ebPx9Dw8\nYmVdNpSVptdx586dpg8Si7dv325s+rejZ4TWcObMmXkpzXr11Ve79rYH9IzRenveV6/gh/ao3h9o\nLLQ30Pmn6Plpcl061Ts+6qczeVHGond87YO8gvqsWbNYmpUQQmKGTpsQQiKCTpsQQiIi58k13bp1\nMzZPTAklmaDEkFGjRhlbdXV1qo2SWnbs2GFsKOGjffv2xqZjXSgOjRJ60Dx0fLBVq1amD6o0iOKK\naF11QhLqg9Ya/SZ0T11lEVVdRMke6J46xvrYY4+ZPjfeeKOx5Yts49eITz75xNgGDhxobPr5Ib0G\nPScUv0YJZvo9RLFXtDeaNWtmbBqkp6B5IdC+1TbvO+GNQ2s/4k2AQuPrdXz66adNn4kTJxobgl/a\nhBASEXTahBASEXTahBASEXTahBASETkXIpGQgar16X4ouea8884zNnTElha5kACCBEbUT4sRIlbE\nLCwsdM0LiXQ6MWf9+vWmT0lJibEh0RTNVScCoeQU9DwKCgqMDaEFV5RohJJEkJj1zDPPpNoPPvig\n6YOEt3yBhDxPAhXikksuMTaPsOat9ojm6qlQh/axN/FK29BRemgfZ1uFD83BK0565o/GQvNHe1QL\nj5MnTzZ9vPBLmxBCIoJOmxBCIoJOmxBCIoJOmxBCIiLnQiQSFFHwXosUSIzwintaYEFz2Lp1q7Gh\nzD0kUOhjsFB2JcoEraqqyjgWug7hyUATsWuBjvBCoCp8TZo0yWhDld3QPdFz08IjypZDzyhfeEVH\ntIc8eLL50NjoPfFWhdR4j4pDVf50dq/3+DQkrnoFUQ/eLEk9D08fEbxmHuHRK8DyS5sQQiKCTpsQ\nQiKCTpsQQiKCTpsQQiIi50IkysRCNp3FhcRDlKWHsr+0UILuh8qkIgEBiU1aANFZhyJ4/j179jQ2\nLcihLDWUsVhUVGRsaH20UIXKvKL16d+/v7Gha3W53BEjRpg+nTp1MrbKykpj06LmDTfcYPq88sor\nxpYvvKU6tVjlFcI8ICHs8OHDxoZK4Xru6c0oREfKea5DIj46VtAzD++6ImEcXbt48eJUe9iwYaYP\n+sOBNWvWGJtm5MiRxrZo0aKM14nwS5sQQqKCTpsQQiKCTpsQQiIi5zFt9EfxKOam454ouWP37t3G\nhmJ6utIc6oPiv94Ype6HKtvt2bPH2FCyiI4FokpoKM6NfpMn0QIlRqDxd+7caWx9+vQxtvLy8lQb\nzX/58uXGVlxcnHEeCxYsMH28yUENgffYKr3mSDtBzyDbCnjeCo1ov+h95U0eQe85mr+mXbt2GfuI\n4HXVeOPv6D3s27evsenficZC8WvPM0LxaybXEELISQidNiGERASdNiGERASdNiGERETOhciWLVsa\n27Zt24ztm2++SbVRRTl9tJUIFqa0ALJ582bTp0OHDsbmOX5JxAoeSADZv3+/saEkEy30oGQDlLCC\nxCAk3mpRCq0XEniR6LJu3Tpj06IOEnlQRUKUtKQrBCIRDCVF/dHRIh0SnFCSmEd8QwIgAomCHuHO\nKzqi+WvQM0dJXOg937hxo7Hp+XuTa1C/VatWZeyXbaVEBLrOWw2SX9qEEBIRdNqEEBIRdNqEEBIR\ndNqEEBIRORciUcZQWVmZsWmRAglmSJxEx4Zp0aW0tDTj/URwFbuSkhJjW7p0aap98cUXmz6o+te+\nffuMTYuHKOvwyy+/NDZUWRDNv7Cw0Ng0KEMVVSREguhrr72WaiPxEImy6FnqZ4L6eMS5hgKJrp4K\nct7f4MmSRAIjGt+bsaizVwcOHGj6eEVTvT5oz6K97T1GTAt36Dok7qH5jx071tjmzZuXcQ7ZVmxE\n74l3X/BLmxBCIoJOmxBCIoJOmxBCIoJOmxBCIiLnQmSbNm2MDYmAWhxAR2yhkpNIKNRlUVEmIhIL\n0FhIPBwwYECqjcQOVJoVlS3V90S/G9mQiITW+qWXXkq1keCCBCI0/+eeey7j3NBc+/XrZ2wo602L\nn1999ZXpg8S/fJFtFqMXdAwcylDUIEHOm7mHhEcPSOjM9pg171zfeeedVPvKK690XYfmikRH/V5v\n2LDB9EGCPfqd3bt3zziWF35pE0JIRNBpE0JIRNBpE0JIRNBpE0JIRORciGzatKmxoQw8nUn29ddf\nu8ZH2Y763DlUzhOJSEePHjU2lFGohQYkGKHxkZCqxUkkNLVt29bYkFCIMgiHDBmSaiPRFAmRSDRF\nZ2Hq0rtIYERjVVVVZeyH1gI9o3zhLaWp1xf9BiSwon2FMuk0uc4aReOjEsyeMrre0qnofRo1alTG\n8dE7h8RyNA8t6GoxUQTPFZVX1nj3DoJf2oQQEhF02oQQEhF02oQQEhEh1/GvRo0amRvMnDnT9NPV\n7lA8D/1RPIoF6jg0Ok5LH20lgquqoXnoe6LKbggUc9YV8FCsC123a9cuY0Pz0JoCSrxAscd3333X\n2K6//vqM16L4OALFQPWRc2hvovj+Qw89lN2ZTydIy5YtzQRvvfVW0y+X7xjaL94qefWJ5/gstA7e\nuXoSc1DCGWL+/PnGNmbMGGPTGot3b6O10PqM93fPmjXLDMYvbUIIiQg6bUIIiQg6bUIIiQg6bUII\niYicJ9cgYQYJg1q0QAIjEiJRZTudgICSeZDoiIQ8lFCiRUCUUOIRZkSsyNi4cWPTB61X586dXf30\nOiJBp3Xr1sZ29913u8bX64iSP1BiBBKN9Fqj9UIVIvPFhAkTjA2tryeRAl3nEdZORHREz8pTRRDN\nNVuhEN3Pm3CjQf4BMWPGDGNbv369sXmER09SDuJExGl+aRNCSETQaRNCSETQaRNCSETQaRNCSETk\nXIhEVfhQNpwW81Am4saNG40NHQdWVFSUaqOMvw4dOmScgwgW6Q4cOJBqo6ptqLrhwoULjW3cuHGp\ndnV1temD1gL9biTEaMEDiU+o6tmmTZuMrXnz5samhRgkFB48eNA1Vw363X8kvMd6aRsSvJEw7hH3\ndEatCH6/kOCHjuHbsmVLqo3+IAAJy4sWLTK2yy67LOMcvIIcWmst8KLxkfi5bt06Y2vWrJmxab+B\nBGVk8/ymbMVWEX5pE0JIVNBpE0JIRNBpE0JIRNBpE0JIRORciEQCCzpKrH379qk2Ok6rrKzMNZYW\nH1BWIzruCokWa9asMTadtYiEn6FDh7rG1wILEnm8x06h8bW4UVlZafr07t3b2JC4imx6/uh5I9FF\nHwknYo9pQvdDGaN/JJAYptcIiXuecp4i9rkjcRg9AyTkbd68OeM80H4cNmyYsSH0WLkuzYrEbc87\nIYLFeC0yItERPTdPpqlnn9QFv7QJISQi6LQJISQi6LQJISQich7TRnEmFKPV8UyU1IKO2EJJMjrB\no6CgwPQpKSkxNhRHRzExHY/dtm2b6bNgwQJj00k/IiI7d+5MtVFlMZRIg45L88TJUCwZxUBRkgyy\n6Rhz165dTR9UHRBVXtQJDmhe3kpuDYE3RqtjoSjm6a3yp/t5dJK6xkfxWJ34g57BkiVLXGN5jtg6\nkdiupj6PA0P9vFUKPbFv7zNC8EubEEIigk6bEEIigk6bEEIigk6bEEIiIudCJBLfVq9ebWw9evRI\ntVGAHwmKSOTSCQco6I+qo6E/sEeihUcMQ5Xctm7damzFxcWp9pEjR0yf0aNHGxv63StWrDA2Lbh6\nxRS01mgdtaiGKjEWFhYam66UKGLXAq0zSvbIF2g90Jz1GnkFJ5SEo5+ft8ocGgsJop7j3Lz39Pzu\n8vJyY0PvZkVFRcbx0by8oqan6l59Ho3mnQOCX9qEEBIRdNqEEBIRdNqEEBIRdNqEEBIRIdugOSGE\nkIaHX9qEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIR\ndNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqEEBIRdNqE\nEBIRdNqEEBIR/wcC9/bxDu4hWAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1113876d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.linspace(1, numIter, numIter)\n",
    "\n",
    "plt.figure()\n",
    "plt.semilogy(t, costISTA, 'r-', t, costFISTA, 'b-')\n",
    "plt.xlim(1, numIter)\n",
    "plt.grid(True)\n",
    "plt.legend([\"ISTA\", \"FISTA\"])\n",
    "plt.title(\"Evolution of the Cost\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"Relative Cost\")\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "# Plot original\n",
    "plt.subplot(1, 2, 1)\n",
    "imgplot = plt.imshow(fhatISTA, interpolation='nearest')\n",
    "imgplot.set_cmap('gray')\n",
    "plt.axis(\"off\")\n",
    "plt.title(\"ISTA\")\n",
    "\n",
    "# Plot reconstruction\n",
    "plt.subplot(1, 2, 2)\n",
    "imgplot = plt.imshow(fhatFISTA, interpolation='nearest')\n",
    "imgplot.set_cmap('gray')\n",
    "plt.axis(\"off\")\n",
    "plt.title(\"FISTA\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}