DouglasRachfordL1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Douglas Rachford Proximal Splitting\n",
    "===================================\n",
    "\n",
    "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n",
    "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n",
    "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n",
    "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n",
    "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n",
    "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n",
    "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n",
    "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n",
    "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n",
    "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$\n",
    "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n",
    "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n",
    "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n",
    "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n",
    "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n",
    "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n",
    "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n",
    "$\\newcommand{\\CC}{\\mathbb{C}}$\n",
    "$\\newcommand{\\RR}{\\mathbb{R}}$\n",
    "$\\newcommand{\\EE}{\\mathbb{E}}$\n",
    "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n",
    "$\\newcommand{\\Ww}{\\mathcal{W}}$\n",
    "$\\newcommand{\\Vv}{\\mathcal{V}}$\n",
    "$\\newcommand{\\Nn}{\\mathcal{N}}$\n",
    "$\\newcommand{\\NN}{\\mathcal{N}}$\n",
    "$\\newcommand{\\Hh}{\\mathcal{H}}$\n",
    "$\\newcommand{\\Bb}{\\mathcal{B}}$\n",
    "$\\newcommand{\\Ee}{\\mathcal{E}}$\n",
    "$\\newcommand{\\Cc}{\\mathcal{C}}$\n",
    "$\\newcommand{\\Gg}{\\mathcal{G}}$\n",
    "$\\newcommand{\\Ss}{\\mathcal{S}}$\n",
    "$\\newcommand{\\Pp}{\\mathcal{P}}$\n",
    "$\\newcommand{\\Ff}{\\mathcal{F}}$\n",
    "$\\newcommand{\\Xx}{\\mathcal{X}}$\n",
    "$\\newcommand{\\Mm}{\\mathcal{M}}$\n",
    "$\\newcommand{\\Ii}{\\mathcal{I}}$\n",
    "$\\newcommand{\\Dd}{\\mathcal{D}}$\n",
    "$\\newcommand{\\Ll}{\\mathcal{L}}$\n",
    "$\\newcommand{\\Tt}{\\mathcal{T}}$\n",
    "$\\newcommand{\\si}{\\sigma}$\n",
    "$\\newcommand{\\al}{\\alpha}$\n",
    "$\\newcommand{\\la}{\\lambda}$\n",
    "$\\newcommand{\\ga}{\\gamma}$\n",
    "$\\newcommand{\\Ga}{\\Gamma}$\n",
    "$\\newcommand{\\La}{\\Lambda}$\n",
    "$\\newcommand{\\si}{\\sigma}$\n",
    "$\\newcommand{\\Si}{\\Sigma}$\n",
    "$\\newcommand{\\be}{\\beta}$\n",
    "$\\newcommand{\\de}{\\delta}$\n",
    "$\\newcommand{\\De}{\\Delta}$\n",
    "$\\newcommand{\\phi}{\\varphi}$\n",
    "$\\newcommand{\\th}{\\theta}$\n",
    "$\\newcommand{\\om}{\\omega}$\n",
    "$\\newcommand{\\Om}{\\Omega}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This numerical tour presents the Douglas-Rachford (DR) algorithm to\n",
    "minimize the sum of two simple functions. It shows an\n",
    "application to\n",
    "reconstruction of exactly sparse signal from noiseless measurement using\n",
    "$\\ell^1$ minimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "%pylab inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Douglas-Rachford Algorithm\n",
    "--------------------------\n",
    "The Douglas-Rachford (DR) algorithm is an iterative scheme to minimize\n",
    "functionals of the form\n",
    "$$ \\umin{x} f(x) + g(x) $$\n",
    "where $f$ and $g$ are convex functions, of which one is able to\n",
    "compute the proximity operators.\n",
    "\n",
    "This algorithm was introduced in\n",
    "\n",
    "P. L. Lions and B. Mercier\n",
    "\"Splitting Algorithms for the Sum of Two Nonlinear Operators,\"\n",
    "_SIAM Journal on Numerical Analysis_\n",
    "vol. 16, no. 6, 1979,\n",
    "\n",
    "\n",
    "as a generalization of an algorithm introduced by Douglas and Rachford in\n",
    "the case of quadratic minimization (which corresponds to solving\n",
    "a positive definite linear system).\n",
    "\n",
    "\n",
    "To learn more about this algorithm, you can read:\n",
    "\n",
    "\n",
    "Patrick L. Combettes and Jean-Christophe Pesquet,\n",
    "\"Proximal Splitting Methods in Signal Processing,\"\n",
    "in: _Fixed-Point Algorithms for Inverse\n",
    "Problems in Science and Engineering_, New York: Springer-Verlag, 2010.\n",
    "\n",
    "\n",
    "\n",
    "The Douglas-Rachford algorithm takes an arbitrary element $s^{(0)}$, a parameter $\\ga>0$, a relaxation parameter $0<\\rho<2$ and iterates, for $k=1,2,\\ldots$\n",
    "$$\n",
    "\\left|\\begin{array}{l}\n",
    "x^{(k)} = \\mathrm{prox}_{\\gamma f} (s^{(k-1)} )\\\\\n",
    "s^{(k)} = s^{(k-1)}+\\rho\\big(\\text{prox}_{\\ga g}( 2x^{(k)}-s^{(k-1)})-x^{(k)}\\big).\n",
    "\\end{array}\\right.\n",
    "$$\n",
    "\n",
    "It is of course possible to inter-change the roles of $f$ and $g$,\n",
    "which defines a different algorithm.\n",
    "\n",
    "The iterates $x^{(k)}$ converge to a solution $x^\\star$ of the problem, i.e. a minimizer of $f+g$.\n",
    "\n",
    "Compressed Sensing Acquisition\n",
    "------------------------------\n",
    "Compressed sensing acquisition corresponds to a random projection\n",
    "$y=Ax^\\sharp$ of a signal $x^\\sharp$ on a\n",
    "few linear vectors (the rows of the matrix $A$). For the recovery of $x^\\sharp$ to be possible, this vector is supposed\n",
    "to be sparse in some basis. Here, we suppose $x^\\sharp$ itself is sparse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We initialize the random number generator for reproducibility."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dimension of the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 400"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = round(N/4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a random Gaussian measurement matrix $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = randn(P,N) / sqrt(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparsity of the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = 17"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by generating a $S$-sparse signal $x^\\sharp$ with $S$ randomized values.\n",
    "Since the measurement matrix is random, one does not care about the sign\n",
    "of the nonzero elements, so we set values equal to one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sel = random.permutation(N)\n",
    "sel = sel[0:S]   # indices of the nonzero elements of xsharp\n",
    "xsharp = zeros(N)\n",
    "xsharp[sel] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We perform random measurements $y=Ax^\\sharp$ without noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = A.dot(xsharp) # matrix-vector multiplication in Python, more precisely dot product of a 2-D array with a 1-D array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compressed Sensing Recovery with the Douglas-Rachford algorithm\n",
    "------------------------------------\n",
    "Compressed sensing recovery corresponds\n",
    "to solving the inverse problem $y=A x^\\sharp$, which is ill posed because\n",
    "$x^\\sharp$ is\n",
    "higher dimensional than $y$.\n",
    "\n",
    "\n",
    "The reconstruction can be performed by $\\ell^1$ minimization,\n",
    "which regularizes the problem by exploiting the prior knowledge that the solution is sparse.\n",
    "$$ x^\\star \\in \\arg\\min_x \\norm{x}_1 \\quad\\mbox{s.t.}\\quad Ax=y$$\n",
    "where the $\\ell^1$ norm is defined as\n",
    "$$ \\norm{x}_1 = \\sum_{n=1}^N \\abs{x_n}. $$\n",
    "\n",
    "\n",
    "This is the minimization of a non-smooth function under affine\n",
    "constraints. This can be shown to be equivalent to a linear programming\n",
    "problem, for wich various algorithms can be used (simplex, interior\n",
    "points). We propose here to use the Douglas-Rachford algorithm.\n",
    "\n",
    "\n",
    "It is possible to recast this problem as the minimization of $f+g$\n",
    "where $g(x) = \\norm{x}_1$ and $f(x)=\\iota_{\\Omega}$ where $\\Omega =\n",
    "\\enscond{x}{Ax=y}$ is an affine space, and $\\iota_\\Omega$ is the indicator\n",
    "function\n",
    "$$ \\iota_\\Omega(x) = \\choice{ 0 \\qifq x \\in \\Omega, \\\\ +\\infty \\qifq x \\notin \\Omega. } $$\n",
    "\n",
    "\n",
    "The proximal operator of the $\\ell^1$ norm is soft thresholding:\n",
    "$$ \\text{prox}_{\\gamma \\norm{\\cdot}_1}(x_n) = \\max\\pa{ 0, 1-\\frac{\\ga}{\\abs{x_n}} } x_n. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prox_gamma_g (x, gamma) :\n",
    "    return x - x/maximum(abs(x)/gamma,1) # soft-thresholding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the 1-D curve of the thresholding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.0, 1.0, -0.80000000000000004, 0.80000000000000004)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAFwCAYAAAAYFxnDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF/5JREFUeJzt3X2wbXd91/HP13sBi4IxE+eGJJeJ2iDgSAXHNFJnOApp\nL0ES6ig01qJY2ozTaNsZbUJRe/9QpzDjqAyWyaSRiWM1ImAHpjzdUq44TAtkCE/lXpJUM+YBLg8t\nqIAziXz94+zQw8k5++xz1n7er9dMJnud9dt7/ZKddfc7v73uutXdAQCYlz+w6AkAAJtFfAAAcyU+\nAIC5Eh8AwFyJDwBgrsQHADBXg+Ojqk5V1fmquq+qbtlj/yVV9b6q+mRVfbaq/vbQYwIAq6uG3Oej\nqo4l+XySlyR5OMnHk9zY3ed2jDmd5Cnd/bqqumQ0/kR3PzZk4gDAahq68nF1kvu7+4HufjTJXUlu\n2DXmC0mePnr89CRfFR4AsLmOD3z+5Uke3LH9UJLv3zXm9iS/UVWPJHlaklcOPCYAsMKGrnxM8p3N\nzyf5ZHdfluTPJvk3VfW0gccFAFbU0JWPh5Oc3LF9MturHzu9MMk/S5Lu/p2q+h9J/lSSu3cOqip/\nyAwArJHurr1+PjQ+7k5yVVVdmeSRJK9KcuOuMeezfUHqR6rqRLbD478fZpLLpKpOd/fpRc+DvXl/\nlp/3aPl5j5bbqrw/4xYVBsVHdz9WVTcneX+SY0nu6O5zVXXTaP9tSf55krdW1aey/TXPz3X37w45\nLgCwuoaufKS735vkvbt+dtuOx19J8vKhxwEA1oM7nB7e2UVPgLHOLnoCHOjsoifAgc4uegKMdXbR\nExhq0E3GpqmqehWu+QAADjbuc93KBwAwV+IDAJgr8QEAzJX4AADmSnwAAHMlPgCAuRIfAMBciQ8A\nYK7EBwAwV+IDAJgr8QEAzJX4AADmSnwAAHMlPgCAuRIfALAiqvK8qjx70fMYSnwAwAqoyvOSvD/J\ncxc9l6HEBwAsuR3h8TPdeeei5zOU+ACAJbYrPP7TouczDeIDAJbUOoZHIj4AYCmta3gk4gMAls46\nh0ciPgBgqax7eCTiAwCWxiaERyI+AGApbEp4JOIDABZuk8IjER8AsFCbFh6J+ACAhdnE8EjEBwAs\nxKaGRyI+AGDuNjk8EvEBAHO16eGRiA8AmBvhsW1wfFTVqao6X1X3VdUt+4zZqqp7quqzVXV26DEB\nYNUIj99X3X30J1cdS/L5JC9J8nCSjye5sbvP7RhzUZKPJPmh7n6oqi7p7q/s8Vrd3XXkyQDAktrE\n8Bj3uT505ePqJPd39wPd/WiSu5LcsGvM30jyju5+KEn2Cg8AWFebGB4HGRoflyd5cMf2Q6Of7XRV\nkour6kNVdXdV/djAYwLAShAeezs+8PmTfGfzpCQvSPLiJE9N8ptV9Vvdfd/AYwPA0hIe+xsaHw8n\nOblj+2S2Vz92ejDJV7r7W0m+VVUfTvJ9SZ4QH1V1esfm2e4+O3B+ADB3mxgeVbWVZGuisQMvOD2e\n7QtOX5zkkSQfyxMvOH12kjcn+aEkT0ny0SSv6u7P7XotF5wCsPI2MTz2Mu5zfdDKR3c/VlU3Z/tf\n8rEkd3T3uaq6abT/tu4+X1XvS/LpJN9Ocvvu8ACAdSA8JjNo5WOarHwAsMqEx3eb5W+1BYCNJzwO\nR3wAwADC4/DEBwAckfA4GvEBAEcgPI5OfADAIQmPYcQHAByC8BhOfADAhITHdIgPAJiA8Jge8QEA\nBxAe0yU+AGAM4TF94gMA9iE8ZkN8AMAehMfsiA8A2EV4zJb4AIAdhMfsiQ8AGBEe8yE+ACDCY57E\nBwAbT3jMl/gAYKMJj/kTHwBsLOGxGOIDgI0kPBZHfACwcYTHYokPADaK8Fg88QHAxhAey0F8ALAR\nhMfyEB8ArD3hsVzEBwBrTXgsH/EBwNoSHstJfACwloTH8hIfAKwd4bHcxAcAa0V4LD/xAcDaEB6r\nQXwAsBaEx+oQHwCsPOGxWgbHR1WdqqrzVXVfVd0yZtyfr6rHquqvDj0mADxOeKyeQfFRVceSvDnJ\nqSTPTXJjVT1nn3FvSPK+JDXkmADwOOGxmoaufFyd5P7ufqC7H01yV5Ib9hj395K8PcmXBx4PAJII\nj1U2ND4uT/Lgju2HRj/7jqq6PNtB8pbRj3rgMQHYcMJjtQ2Nj0lC4l8lubW7O9tfufjaBYAjEx6r\n7/jA5z+c5OSO7ZPZXv3Y6c8luauqkuSSJC+tqke7+127X6yqTu/YPNvdZwfOD4A1IjyWV1VtJdma\naOz2gsSRD3Q8yeeTvDjJI0k+luTG7j63z/i3Jnl3d79zj33d3VZFANiT8Fgt4z7XB618dPdjVXVz\ntv9jOJbkju4+V1U3jfbfNuT1ASARHutm0MrHNFn5AGAvwmM1jftcd4dTAJaW8FhP4gOApSQ81pf4\nAGDpCI/1Jj4AWCrCY/2JDwCWhvDYDOIDgKUgPDaH+ABg4YTHZhEfACyU8Ng84gOAhREem0l8ALAQ\nwmNziQ8A5k54bDbxAcBcCQ/EBwBzIzxIxAcAcyI8eJz4AGDmhAc7iQ8AZkp4sJv4AGBmhAd7ER8A\nzITwYD/iA4CpEx6MIz4AmCrhwUHEBwBTIzyYhPgAYCqEB5MSHwAMJjw4DPEBwCDCg8MSHwAcmfDg\nKMQHAEciPDgq8QHAoQkPhhAfAByK8GAo8QHAxIQH0yA+AJiI8GBaxAcABxIeTJP4AGAs4cG0iQ8A\n9iU8mIXB8VFVp6rqfFXdV1W37LH/R6vqU1X16ar6SFU9b+gxAZg94cGsVHcf/clVx5J8PslLkjyc\n5ONJbuzuczvG/IUkn+vur1fVqSSnu/uaPV6ru7uOPBkApkZ4MNS4z/WhKx9XJ7m/ux/o7keT3JXk\nhp0Duvs3u/vro82PJrli4DEBmCHhwawNjY/Lkzy4Y/uh0c/28+NJ3jPwmADMiPBgHo4PfP7E39lU\n1V9K8neS/MDAYwIwA8KDeRkaHw8nOblj+2S2Vz++y+gi09uTnOru39vvxarq9I7Ns919duD8AJiA\n8GCoqtpKsjXR2IEXnB7P9gWnL07ySJKP5YkXnD4zyW8k+Zvd/VtjXssFpwALIDyYhXGf64NWPrr7\nsaq6Odv/0R5Lckd3n6uqm0b7b0vyT5L80SRvqaokebS7rx5yXACmQ3iwCINWPqbJygfAfAkPZmmW\nv9UWgBUkPFgk8QGwYYQHiyY+ADaI8GAZiA+ADSE8WBbiA2ADCA+WifgAWHPCg2UjPgDWmPBgGYkP\ngDUlPFhW4gNgDQkPlpn4AFgzwoNlJz4A1ojwYBWID4A1ITxYFeIDYA0ID1aJ+ABYccKDVSM+AFaY\n8GAViQ+AFSU8WFXiA2AFCQ9WmfgAWDHCg1UnPgBWiPBgHYgPgBUhPFgX4gNgBQgP1on4AFhywoN1\nIz4AlpjwYB2JD4AlJTxYV+IDYAkJD9aZ+ABYMsKDdSc+AJaI8GATiA+AJSE82BTiA2AJCA82ifgA\nWDDhwaYRHwALJDzYROIDYEGEB5tKfAAsgPBgkw2Oj6o6VVXnq+q+qrplnzFvGu3/VFU9f+gxAVaZ\n8GDTDYqPqjqW5M1JTiV5bpIbq+o5u8Zcl+R7u/uqJD+Z5C1DjgmwyoQHDF/5uDrJ/d39QHc/muSu\nJDfsGnN9kjuTpLs/muSiqjox8LgAK0d4wLbjA59/eZIHd2w/lOT7JxhzRZILA48NK6kqleTSJLXo\nuTBXVyZ5R4QHDI6PnnDc7l9k93xeVZ3esXm2u88eYU6wtKpyLMlbk7w8ybcWPB3m67EID9ZYVW0l\n2Zpk7ND4eDjJyR3bJ7O9sjFuzBWjnz1Bd58eOB9YWjvC47Ikl3fnmwueEsDUjBYMzj6+XVW/sN/Y\nodd83J3kqqq6sqqenORVSd61a8y7krx6NJFrknytu33lwkbZFR7XCw9gkw1a+ejux6rq5mxfQHUs\nyR3dfa6qbhrtv62731NV11XV/Um+keQ1g2cNK0R4AHy36p70so3ZqqrubhfgsVaEB7Cpxn2uu8Mp\nzIjwANib+IAZEB4A+xMfMGXCA2A88QFTJDwADiY+YEqEB8BkxAdMgfAAmJz4gIGEB8DhiA8YQHgA\nHJ74gCMSHgBHIz7gCIQHwNGJDzgk4QEwjPiAQxAeAMOJD5iQ8ACYDvEBExAeANMjPuAAwgNgusQH\njCE8AKZPfMA+hAfAbIgP2IPwAJgd8QG7CA+A2RIfsIPwAJg98QEjwgNgPsQHRHgAzJP4YOMJD4D5\nEh9sNOEBMH/ig40lPAAWQ3ywkYQHwOKIDzaO8ABYLPHBRhEeAIsnPtgYwgNgOYgPNoLwAFge4oO1\nJzwAlov4YK0JD4DlMyg+quriqjpTVfdW1Qeq6qI9xpysqg9V1W9X1Wer6u8POSZMSngALKehKx+3\nJjnT3c9K8sHR9m6PJvnZ7v7TSa5J8lNV9ZyBx4WxhAfA8hoaH9cnuXP0+M4kr9g9oLu/2N2fHD3+\nP0nOZfsDAWZCeAAst6HxcaK7L4weX0hyYtzgqroyyfOTfHTgcWFPwgNg+R0/aEBVnUly6R67Xr9z\no7u7qnrM6/zhJG9P8tOjFRCYKuEBsBoOjI/uvna/fVV1oaou7e4vVtUzknxpn3FPSvKOJP++u391\nzOud3rF5trvPHjQ/SIQHwKJV1VaSrYnGdu+7WDHJgd6Y5Kvd/YaqujXJRd19664xle3rQb7a3T87\n5rW6u+vIk2FjCQ+A5TPuc31ofFyc5G1JnpnkgSSv7O6vVdVlSW7v7pdV1V9M8uEkn07y+MFe193v\nm3SSsB/hAbCcZhYf0yQ+OCzhAbC8xn2uu8MpK0l4AKwu8cHKER4Aq018sFKEB8DqEx+sDOEBsB7E\nBytBeACsD/HB0hMeAOtFfLDUhAfA+hEfLC3hAbCexAdLSXgArC/xwdIRHgDrTXywVIQHwPoTHywN\n4QGwGcQHS0F4AGwO8cHCCQ+AzSI+WCjhAbB5xAcLIzwANpP4YCGEB8DmEh/MnfAA2Gzig7kSHgCI\nD+ZGeACQiA/mRHgA8DjxwcwJDwB2Eh/MlPAAYDfxwcwIDwD2Ij6YCeEBwH7EB1MnPAAYR3wwVcID\ngIOID6ZGeAAwCfHBVAgPACYlPhhMeABwGOKDQYQHAIclPjgy4QHAUYgPjkR4AHBUR46Pqrq4qs5U\n1b1V9YGqumjM2GNVdU9Vvfuox2N5CA8Ahhiy8nFrkjPd/awkHxxt7+enk3wuSQ84HktAeAAw1JD4\nuD7JnaPHdyZ5xV6DquqKJNcl+eUkNeB4LJjwAGAahsTHie6+MHp8IcmJfcb9yyT/MMm3BxyLBRMe\nAEzL8XE7q+pMkkv32PX6nRvd3VX1hK9UquqvJPlSd99TVVtDJsriCA8ApmlsfHT3tfvtq6oLVXVp\nd3+xqp6R5Et7DHthkuur6rokfzDJ06vq33X3q/d5zdM7Ns9299mD/gGYLeEBwCRGiwxbE43tPto1\noFX1xiRf7e43VNWtSS7q7n0vOq2qFyX5B9398n32d3e7JmSJCA8Ajmrc5/qQaz5+Mcm1VXVvkr88\n2k5VXVZVv7bPc/xulxUhPACYlSOvfEyblY/lITwAGGpWKx+sIeEBwKyJD75DeAAwD+KDJMIDgPkR\nHwgPAOZKfGw44QHAvImPDSY8AFgE8bGhhAcAiyI+NpDwAGCRxMeGER4ALJr42CDCA4BlID42hPAA\nYFmIjw0gPABYJuJjzQkPAJaN+FhjwgOAZSQ+1pTwAGBZiY81JDwAWGbiY80IDwCWnfhYI8IDgFUg\nPtaE8ABgVYiPNSA8AFgl4mPFCQ8AVo34WGHCA4BVJD5WlPAAYFWJjxUkPABYZeJjxQgPAFad+Fgh\nwgOAdSA+VoTwAGBdiI8VIDwAWCfiY8kJDwDWjfhYYsIDgHUkPpaU8ABgXYmPJSQ8AFhn4mPJCA8A\n1t2R46OqLq6qM1V1b1V9oKou2mfcRVX19qo6V1Wfq6prjj7d9SY8ANgEQ1Y+bk1yprufleSDo+29\n/Osk7+nu5yR5XpJzA465toQHAJuiuvtoT6w6n+RF3X2hqi5Ncra7n71rzB9Jck93/4kJXq+7u440\nmRUnPABYN+M+14esfJzo7gujxxeSnNhjzB9P8uWqemtVfaKqbq+qpw445toRHgBsmrHxMbqm4zN7\n/HX9znG9vXyy1xLK8SQvSPJL3f2CJN/I/l/PbBzhAcAmOj5uZ3dfu9++qrpQVZd29xer6hlJvrTH\nsIeSPNTdHx9tvz1j4qOqTu/YPNvdZ8fNb5UJDwDWSVVtJdmaaOyAaz7emOSr3f2Gqro1yUXd/YSw\nqKoPJ3ltd987iovv6e5b9hi3Mdd8CA8A1t24z/Uh8XFxkrcleWaSB5K8sru/VlWXJbm9u182Gvd9\nSX45yZOT/E6S13T31w8zyXUiPADYBDOJj2nbhPgQHgBsiln9bhcOQXgAwDbxMQfCAwB+n/iYMeEB\nAN9NfMyQ8ACAJxIfMyI8AGBv4mMGhAcA7E98TJnwAIDxxMcUCQ8AOJj4mBLhAQCTER9TIDwAYHLi\nYyDhAQCHIz4GEB4AcHji44iEBwAcjfg4AuEBAEcnPg5JeADAMOLjEIQHAAwnPiZUlUpyR4QHAAxy\nfNETWBXd6aq8I8kHhQcAHF1196LnkCSpqu7uWvQ8AIDhxn2u+9oFAJgr8QEAzJX4AADmSnwAAHMl\nPgCAuRIfAMBciQ8AYK7EBwAwV+IDAJgr8QEAzJX4AADmSnwAAHMlPgCAuRIfAMBcHTk+quriqjpT\nVfdW1Qeq6qJ9xr2uqn67qj5TVf+hqp5y9OkCAKtuyMrHrUnOdPezknxwtP1dqurKJD+R5AXd/WeS\nHEvyIwOOuXBVtbXoObA/78/y8x4tP+/RcluH92dIfFyf5M7R4zuTvGKPMf8ryaNJnlpVx5M8NcnD\nA465DLYWPQHG2lr0BDjQ1qInwIG2Fj0Bxtpa9ASGGhIfJ7r7wujxhSQndg/o7t9N8i+S/M8kjyT5\nWnf/+oBjAgAr7vi4nVV1Jsmle+x6/c6N7u6q6j2e/yeT/EySK5N8Pcl/rqof7e5fOfKMAYCVVt1P\naIbJnlh1PslWd3+xqp6R5EPd/exdY16V5Nrufu1o+8eSXNPdP7XH6x1tIgDAUuru2uvnY1c+DvCu\nJH8ryRtGf//VPcacT/KPq+p7kvzfJC9J8rHDTBAAWC9DVj4uTvK2JM9M8kCSV3b316rqsiS3d/fL\nRuN+Lttx8u0kn0jy2u5+dApzBwBW0JHjAwDgKNzhdIyq+uujG6T9v6p6wZhxp6rqfFXdV1W3zHOO\nm+4QN7t7oKo+XVX3VNWeX/0xXZOcF1X1ptH+T1XV8+c9x0120PtTVVtV9fXROXNPVf2jRcxzU1XV\nv62qC1X1mTFjVvb8ER/jfSbJDyf58H4DqupYkjcnOZXkuUlurKrnzGd6ZIKb3Y10ti+Qfn53Xz23\n2W2oSc6Lqrouyfd291VJfjLJW+Y+0Q11iF+3/uvonHl+d//TuU6St2b7/dnTqp8/4mOM7j7f3fce\nMOzqJPd39wOja1nuSnLD7GfHyCQ3u3uci5rnZ5Lz4jvvXXd/NMlFVfWE+wUxE5P+uuWcWZDu/m9J\nfm/MkJU+f8THcJcneXDH9kOjnzEfB97sbqST/HpV3V1VPzGfqW20Sc6LvcZcMeN5sW2S96eTvHC0\npP+eqnru3GbHJFb6/BnyW23Xwpgbqf18d797gpdwxe6MDb3Z3cgPdPcXquqPJTlTVedH/2fBbEx6\nXuz+P2vn03xM8u/5E0lOdvc3q+ql2b6dwrNmOy0OaWXPn42Pj+6+duBLPJzk5I7tk9kuUKZk3Hs0\nuiDr0h03u/vSPq/xhdHfv1xV/yXby87iY3YmOS92j7kiq/9nP62KA9+f7v7fOx6/t6p+qaouHv2x\nGSzeSp8/vnaZ3H7ffd6d5KqqurKqnpzkVdm+ARvz8fjN7pJ9bnZXVU+tqqeNHv+hJD+Y7YuJmZ1J\nzot3JXl1klTVNdn+s58uhHk48P2pqhNVVaPHV2f71gzCY3ms9Pmz8Ssf41TVDyd5U5JLkvxaVd3T\n3S/deSO17n6sqm5O8v4kx5Lc0d3nFjjtTfOLSd5WVT+e0c3ukmTXze4uTfLO0a+jx5P8Snd/YDHT\n3Qz7nRdVddNo/23d/Z6quq6q7k/yjSSvWeCUN8ok70+Sv5bk71bVY0m+meRHFjbhDVRV/zHJi5Jc\nUlUPJvmFJE9K1uP8cZMxAGCufO0CAMyV+AAA5kp8AABzJT4AgLkSHwDAXIkPAGCuxAcAMFfiAwCY\nq/8PAFIY33RARSoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f18913d98d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(9,6)\n",
    "t = arange(-1,1,0.001)\n",
    "plot(t, prox_gamma_g(t,0.3))\n",
    "axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The proximity operator of $\\gamma$ times the indicator function of $\\Omega$ is projection onto $\\Omega$ \n",
    " and does not depends on $\\gamma$.\n",
    "$$ \\mathrm{prox}_{\\gamma f}(x)=\\mathrm{prox}_{\\iota_\\Omega}(x)=P_\\Omega(x)  = x + A^* (A A^*)^{-1} (y-Ax). $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "pA = pinv(A) # pseudo-inverse. Equivalent to pA = A.T.dot(inv(A.dot(A.T)))\n",
    "def prox_f (x, y) :\n",
    "    return x + pA.dot(y-A.dot(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the values of $\\gamma$ and $\\rho$.\n",
    "Try different values to speed up the convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = 0.1 # try 1, 10, 0.1\n",
    "rho = 1     # try 1, 1.5, 1.9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbiter = 700"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>\n",
    " <font color='magenta'>\n",
    "__Exercise: Implement nbiter iterations of the Douglas-Rachford algorithm.\n",
    "Keep track of the evolution of the $\\ell^1$ norm.__</font></h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = zeros(N)\n",
    "En_array = zeros(nbiter)\n",
    "skm1=s\n",
    "for iter in range(nbiter):  # iter goes from 0 to nbiter-1\n",
    "    # put your code here\n",
    "    x=prox_f( skm1, y )\n",
    "    skm1=s\n",
    "    s = skm1+gamma*(prox_gamma_g(2*x-skm1,gamma)-x)\n",
    "    \n",
    "    En_array[iter] = norm(x, ord=1)  \n",
    "x_restored = x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We display the original and the recovered signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f188f264450>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6YAAAG0CAYAAAAhAzpSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH6hJREFUeJzt3Xuw7WdZH/DvQwLFaBxk6ERJjgLKKGgVRKP1mpQqJ2mH\naKdVo7VWRTJto87oKJXplDDVtjrjZSiWUhsZLx3ijLZOHBEoxTM6KlBsEkATTFCmuWC8oyKMiXn7\nx14nZ2VnX9bZ+7f281vrfD4ze85Zt/f3vOt9z37O97fWXrvGGAEAAIAuT+guAAAAgAubYAoAAEAr\nwRQAAIBWgikAAACtBFMAAABaCaYAAAC0EkwBAABodXF3AcDjVdWpJJXk8jHGb3TXAwAA6+QVU5in\nJyR5epIv6y4EAADWTTCFefqyJN+b5EPdhQAAwLrVGKO7BmCXqnpmkkeSXDHG+LXuegAAYJ28Ygrz\n9Eh3AQAAcFIEU5ivJyb5jO4iAABg3XwqLzSrqhcn+ZskX5zk3UlOJ7k7yVcm+d+NpQHAxtunz37f\nGOOu1sKAx/AzptCoqj4xyZPGGPdU1W8meWGSL0ryx0nuS/L0McbbO2sEgE11QJ996xjjr06ohucm\nOTXGeNNJHA82lWAKM1BVlyX5mTHGVd21AMC2OW6frapPHmO874iPvSjJT48xrj/K4+FC4WdMoVFV\nfVpVfVaSa5P8yuK6a3urAoDtMEWfXXxS/nVHrWGM8TdJPnLUx8OFws+YQq8vT3Jpkg8keXJVfUV2\n3sILABzfnn22qr45Oz9z+sEk35nkk5O8IMnHZeeT8e9I8vljjB9KcmWSz10E3A8l+bvZ+ZGbv1o8\n7uw4P5rkcxZjvDXJQ4vj35bkkpOYLGwywRQajTFe1V0DAGyr/fpsVf1hkueOMb5zcfnlSf5zdkLs\nZ2fnFdY3LO7+tiSfNsa4o6pel+QlY4y/qaqfSfLdZ8epqp/MTji9NMmpJF+b5F+MMT5cVS9d3yxh\nOwimAABcaCrJny9dHkneN8b4o6r6cJI/SfIdSX59cVuq6tOTfHSSJyX58OLvWRrnkaUxLkvyddn5\nNGBgBT78CACAC0pV/bMk/yTJS8cYH6iqT07yj7PzttvPSfJrST5+jPEzVfXEJD+Q5OeSPJjk6iS/\nn5236v7ts+Nk5+26Z8d4b5InJ3lRdn5Fzfcn+aoxxvtPao6waQRTAAAAWvlUXgAAAFqd2M+YVpWX\nZgGY1BijumvYZHozAFM7am8+0Q8/2pb/QFTVTWOMm7rrOK5tmUdiLnO1LXPZlnkkWzcXoWoCevO8\nbMs8EnOZq22Zy7bMI9m6uRy5N3srLwAAAK0EUwAAAFoJpkdzpruAiZzpLmBCZ7oLmNCZ7gImdKa7\ngImc6S5gQme6C4A1OdNdwETOdBcwoTPdBUzoTHcBEzrTXcBEznQXMKEz3QXMwYn9upiqGtvycywA\n9NNXjs9zCMCUjtNXvGIKAABAK8EUAACAVoIpAAAArQRTAAAAWgmmAAAAtBJMAQAAaCWYAgAA0Eow\nBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGglmAIAANBKMAUAAKCVYAoA\nAEArwRQAAIBWgikAAACtBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAKAABAq0ODaVX9\neFU9WFXvPuA+r6qqu6vqjqp6/rQlAgDL9GYAts0qr5i+Lsnp/W6sqmuTfMoY49lJXprkNRPVBgDs\nTW8GYKtcfNgdxhi/WlXPOOAuL07yE4v7vr2qnlJVl40xHpymxHmr+oxrk1Pfllz65OQvPpLc+6ox\n3vOGdY8xxXGndtI1ret461iPOazXttfQsR+m3isdazSHfbGKversrqmT3nwwvbmvJr35/Gx7DXrz\n9PObk8l78xjj0K8kz0jy7n1u+4UkX7B0+S1JXrDH/cYqx9qkr+TTr01ecncyxrmvl9ydfPq16xxj\niuPO8bmYw/HWsR5zWK9tr6FjP0y9VzrWaA774nh1ZnTX1vu86M3nt1/0Zr1Zbz7JGvTmeT1vJ1dn\nxpHHXO3Ahza/L1y6/JYkn73H/Y5c5Fy/ktNvfOxinP265pfWOcYUx53jczGH461jPeawXtteQ8d+\nmHqvdKzRHPbF8erM6K6t93nRm89vv+jNevO5x8xhvba9Br15Xs/bydWZcdQxD30r7wruT3Jq6fIV\ni+sep6puWrp4ZoxxZoLjN7r0yXtf/zEftd4xpjju1E66pnUdbx3rMYf12vYa5rQfpv63u841msO+\nWMXZOs8svliB3vw4evM5evM81mvba5jTftCbpzd9b54imN6a5MYkt1TV5yf5s7HPz7CMMW6a4Hgz\n8hcf2fv6v/zweseY4rhTO+ma1nW8dazHHNZr22vo2A+jzv+YU493XHPYF6s4W+dVi6+zXnnypWwO\nvflx9OZz9OZ5rNe216A3H80c9sUqpu/Nq/y6mNcn+fUkn1pV91bVN1XVDVV1Q5KMMd6Q5Her6p4k\nr03yL49czca591XJt9zz2Ote8r7k//2n9Y4xxXGndtI1ret461iPOazXttfQsR+m3isdazSHfbGK\n/eq8cOnNB9Gbz9Gb93/MHNZr22vQm49mDvtiFWvozSf3PuSjv994zl87P/h7zS+dey/6+f9g8lHG\nmOK4c3wu5nC8dazHHNZr22vo2A9T75WONZrDvjhqndvaV072ed3O51Bv7qtJb57H8zWXGvTmeT1v\nJ1HncfpK7Qy6flU1xtjv5fDNV7XzfJ70GFMcd2onXdO6jreO9ZjDem17DR37Yeq90rFGc9gXq1iu\nc9v7yknY9udQbz5Hb9abO2vQm49mDvtiFVP15kPfygsAAADrJJgCAADQSjAFAACglWAKAABAK8EU\nAACAVoIpAAAArQRTAAAAWgmmAAAAtBJMAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAA\nAK0EUwAAAFoJpgAAALQSTAEAAGglmAIAANBKMAUAAKCVYAoAAEArwRQAAIBWgikAAACtBFMAAABa\nCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAKAABAK8EUAACAVoIpAAAArQRTAAAAWgmmAAAAtBJM\nAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGglmAIA\nANBKMAUAAKCVYAoAAEArwRQAAIBWgikAAACtBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAFAACg\nlWAKAABAK8EUAACAVoIpAAAArQRTAAAAWgmmAAAAtBJMAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvB\nFAAAgFaHBtOqOl1Vd1XV3VX1sj1uf1pVvbGqbq+q91TVP19LpQBAEr0ZgO1zYDCtqouSvDrJ6STP\nTXJ9VT1n191uTHLbGON5Sa5K8oNVdfEaagWAC57eDMA2OuwV0yuT3DPGeP8Y46EktyS5btd9PpDk\nYxd//9gkfzzGeHjaMgGABb0ZgK1z2NnTy5Pcu3T5viSft+s+P5bkrVX1QJJLk3zVdOUBALvozQBs\nncNeMR0rjPHyJLePMZ6e5HlJfrSqLj12ZQDAXvRmALbOYa+Y3p/k1NLlU9k5M7vsC5J8X5KMMd5X\nVb+X5FOTvHP3YFV109LFM2OMM+dZLwAXqKq6KnlFql55U3ctzfRmAGZhyt5cY+x/4nXxQQnvTfLC\nJA8keUeS68cYdy7d54eSfHCM8cqquizJbyb5zDHGn+waa4wx6rgFz1VVxhg51vyOMsYUx53aSde0\nruOtYz3msF7bXkPHfph6r3Ss0Rz2xSqW69z2vrIfvXl1evM5erPe3FmD3nw0c9gXq5iqNx/4iukY\n4+GqujHJm5JclOTmMcadVXXD4vbXJvn3SV5XVXdk563B37278QEA09CbAdhGB75iOumBnJVdyxhz\nPJPirKyzsp01OCt7NHPYF6vwium0tv051JvP0Zv15s4a9OajmcO+WMVUvfmwDz8CAACAtRJMAQAA\naCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGglmAIAANBK\nMAUAAKCVYAoAAEArwRQAAIBWgikAAACtBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAK\nAABAK8EUAACAVoIpAAAArQRTAAAAWgmmAAAAtBJMAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAA\ngFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGglmAIAANBKMAUAAKCVYAoAAEArwRQAAIBWgikAAACt\nBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAKAABAK8EUAACAVoIpAAAArQRTAAAAWgmm\nAAAAtBJMAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEA\nAGglmAIAANBKMAUAAKCVYAoAAEArwRQAAIBWhwbTqjpdVXdV1d1V9bJ97nNVVd1WVe+pqjOTVwkA\nPEpvBmDbXHzQjVV1UZJXJ/n7Se5P8n+q6tYxxp1L93lKkh9N8qIxxn1V9bR1FgwAFzK9GYBtdNgr\nplcmuWeM8f4xxkNJbkly3a77fG2Snxtj3JckY4w/mr5MAGBBbwZg6xwWTC9Pcu/S5fsW1y17dpKn\nVtUvV9U7q+rrpywQAHgMvRmArXPgW3mTjBXGeGKSz07ywiSXJPmNqnrbGOPu4xYHADyO3gzA1jks\nmN6f5NTS5VPZOTO77N4kfzTG+HCSD1fVryT5rCSPa35VddPSxTNjjDPnWzAAF6aquip5RapeeVN3\nLc30ZgBmYcreXGPsf+K1qi5O8t7snHF9IMk7kly/6wMWPi07H8LwoiR/K8nbk3z1GOO3d401xhh1\n3ILnqipjjBxrfkcZY4rjTu2ka1rX8daxHnNYr22voWM/TL1XOtZoDvtiFct1bntf2Y/evDq9+Ry9\nWW/urEFvPpo57ItVTNWbD3zFdIzxcFXdmORNSS5KcvMY486qumFx+2vHGHdV1RuTvCvJI0l+bHfj\nAwCmoTcDsI0OfMV00gM5K7uWMeZ4JsVZWWdlO2twVvZo5rAvVuEV02lt+3OoN5+jN+vNnTXozUcz\nh32xiql682GfygsAAABrJZgCAADQSjAFAACglWAKAABAK8EUAACAVoIpAAAArQRTAAAAWgmmAAAA\ntBJMAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGgl\nmAIAANBKMAUAAKCVYAoAAEArwRQAAIBWgikAAACtBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAF\nAACglWAKAABAK8EUAACAVoIpAAAArQRTAAAAWgmmAAAAtBJMAQAAaCWYAgAA0EowBQAAoJVgCgAA\nQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGglmAIAANBKMAUAAKCVYAoAAEArwRQAAIBW\ngikAAACtBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAKAABAK8EUAACAVoIpAAAArQRT\nAAAAWgmmAAAAtBJMAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0ODaZVdbqq7qqq\nu6vqZQfc73Or6uGq+kfTlggALNObAdg2BwbTqrooyauTnE7y3CTXV9Vz9rnf9yd5Y5JaQ50AQPRm\nALbTYa+YXpnknjHG+8cYDyW5Jcl1e9zvW5P8bJI/nLg+AOCx9GYAts5hwfTyJPcuXb5vcd2jqury\n7DTE1yyuGpNVBwDspjcDsHUOC6arNLIfSfKvxxgjO28V8nYhAFgfvRmArXPxIbffn+TU0uVT2Tkz\nu+wFSW6pqiR5WpJrquqhMcatuwerqpuWLp4ZY5w534IBuDBV1VXJK1L1ypu6a2mmNwMwC1P25to5\nmbrvgS5O8t4kL0zyQJJ3JLl+jHHnPvd/XZJfGGP8jz1uG2OMrT1jW5UxxvHOSB9ljCmOO7WTrmld\nx1vHesxhvba9ho79MPVe6VijOeyLVSzXue19ZT968+r05nP0Zr25swa9+WjmsC9WMVVvPvAV0zHG\nw1V1Y5I3Jbkoyc1jjDur6obF7a89ykEBgKPRmwHYRge+YjrpgZyVXcsYczyT4qyss7KdNTgrezRz\n2Ber8IrptLb9OdSbz9Gb9ebOGvTmo5nDvljFVL35sA8/AgAAgLUSTAEAAGglmAIAANBKMAUAAKCV\nYAoAAEArwRQAAIBWgikAAACtBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAKAABAK8EU\nAACAVoIpAAAArQRTAAAAWgmmAAAAtBJMAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAA\nAK0EUwAAAFoJpgAAALQSTAEAAGglmAIAANBKMAUAAKCVYAoAAEArwRQAAIBWgikAAACtBFMAAABa\nCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAKAABAK8EUAACAVoIpAAAArQRTAAAAWgmmAAAAtBJM\nAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGglmAIA\nANBKMAUAAKCVYAoAAEArwRQAAIBWgikAAACtBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAFAACg\nlWAKAABAK8EUAACAVisF06o6XVV3VdXdVfWyPW7/uqq6o6reVVW/VlWfOX2pAMBZejMA2+TQYFpV\nFyV5dZLTSZ6b5Pqqes6uu/1uki8ZY3xmkn+X5L9OXSgAsENvBmDbrPKK6ZVJ7hljvH+M8VCSW5Jc\nt3yHMcZvjDE+uLj49iRXTFsmALBEbwZgq6wSTC9Pcu/S5fsW1+3nm5O84ThFAQAH0psB2CoXr3Cf\nsepgVXV1km9K8oVHrggAOIzeDMBWWSWY3p/k1NLlU9k5M/sYiw9V+LEkp8cYf7rXQFV109LFM2OM\nMytXCsAFraquSl6Rqlfe1F3LDOjNALSbsjfXGAefdK2qi5O8N8kLkzyQ5B1Jrh9j3Ll0n09M8tYk\n/3SM8bZ9xhljjDpuwXNVlTFGjjW/o4wxxXGndtI1ret461iPOazXttfQsR+m3isdazSHfbGK5Tq3\nva8cRG9ejd58jt6sN3fWoDcfzRz2xSqm6s2HvmI6xni4qm5M8qYkFyW5eYxxZ1XdsLj9tUn+bZKP\nS/KaqkqSh8YYVx6lIADgYHozANvm0FdMJzuQs7JrGWOOZ1KclXVWtrMGZ2WPZg77YhVeMZ3Wtj+H\nevM5erPe3FmD3nw0c9gXq5iqN6/yqbwAAACwNoIpAAAArQRTAAAAWgmmAAAAtBJMAQAAaCWYAgAA\n0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGglmAIAANBKMAUAAKCV\nYAoAAEArwRQAAIBWgikAAACtBFMAAABaCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAKAABAK8EU\nAACAVoIpAAAArQRTAAAAWgmmAAAAtBJMAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAA\nAK0EUwAAAFoJpgAAALQSTAEAAGglmAIAANBKMAUAAKCVYAoAAEArwRQAAIBWgikAAACtBFMAAABa\nCaYAAAC0EkwBAABoJZgCAADQSjAFAACglWAKAABAK8EUAACAVoIpAAAArQRTAAAAWgmmAAAAtBJM\nAQAAaCWYAgAA0EowBQAAoJVgCgAAQCvBFAAAgFaCKQAAAK0EUwAAAFoJpgAAALQSTAEAAGglmAIA\nANDq0GBaVaer6q6quruqXrbPfV61uP2Oqnr+9GUCAGfpzQBsmwODaVVdlOTVSU4neW6S66vqObvu\nc22STxljPDvJS5O8Zk21zkZVXdVdwxS2ZR6JuczVtsxlW+aRbNdcLlR68962ZW9vyzwSc5mrbZnL\ntswj2a65HMfFh9x+ZZJ7xhjvT5KquiXJdUnuXLrPi5P8RJKMMd5eVU+pqsvGGA/uHqzq6oeS1BSF\nH98lSR766+QjjyQXVZInZeXaPukJVVc/8thx3pyqqz90fuPsruV8xpjiuMvz2D3uUZ6XKWo6H+eO\nl3zSm6uurmmOt9c8LqmDn5fD5n4+z81e6zKFdazPYftl91zWuUfWNfYlSZ45kuwx7tR75ST+DZ3k\nmiw7zveW3XX+gz9Pnvzw9DVujC3uzWcdZb9sc29eHvt8/w1tc29OHfy8XKi9eXncvZ6XC7E3775t\n9/Nyofbm5WMd9//+x+/NhwXTy5Pcu3T5viSft8J9rkjyuOaXPOuw452Qj0/yoiQ/8VFHe/yfJXnW\nRY8f51mXHL+WVcaY4rjJuXnsN+75mKqmox7v456YPGsN4z7rksOfl8Pmfr7Pze51mcI61meV/bI8\nl3XukXWNfXbcb8/jx516r5zUv6GTWpNlx/2eu/vxv3jpzp8zy1InZ0t781lH3S/b2pv3GntV29yb\n97p+lcesevtum9Kb9xp3twutN+9321Ees029ea9jHeexx+/NhzWjseI4uyvY53GfsOJw6/a9Sf5N\njl7PpYvH7h7nKOMdZYwpjpucm8d+456PqWo66vF2z2Wqcfda51Uecz637zbVXI5Tw1HG3MvyXNa5\nR9Y19tlxL91j3Kn3ykn9GzqpNVl23O+5x3381tnS3nzWUdd7W3vzXmOvapt7817Xr/KYVW/fbVN6\n817j7nah9eb9bjvKY7apN+91rJN67N5qjP37W1V9fpKbxhinF5e/J8kjY4zvX7rPf0lyZoxxy+Ly\nXUm+dPfbhapq1UYKACsZY1xwL5vqzQDM2VF782GvmL4zybOr6hlJHkjy1Umu33WfW5PcmOSWRbP8\ns71+huVC/M8DAKyB3gzA1jkwmI4xHq6qG5O8KclFSW4eY9xZVTcsbn/tGOMNVXVtVd2T5ENJvnHt\nVQPABUpvBmAbHfhWXgAAAFi3A3+P6RRW+SXgc1ZV76+qd1XVbVX1jsV1T62q/1VVv1NVb66qp3TX\nuZeq+vGqerCq3r103b61V9X3LNbprqr68p6q97bPXG6qqvsWa3NbVV2zdNss51JVp6rql6vqt6rq\nPVX1bYvrN25dDpjLRq1LVT25qt5eVbdX1W9X1X9YXL+Ja7LfXDZqTZZV1UWLmn9hcXnj1mWO9OY+\n29Kbt6UvJ3rzHOeiN89zLmetrTePMdb2lZ23GN2T5BlJnpjk9iTPWecx1zCH30vy1F3X/UCS7178\n/WVJ/mN3nfvU/sVJnp/k3YfVnp1f0n77Yp2esVi3J3TP4ZC5vCLJd+xx39nOJTufrf28xd8/Jsl7\nkzxnE9flgLls4rpcsvjz4iRvS/JFm7gmB8xl49ZkqcbvSPLfk9y6uLyR6zKnr+jN3bVvRW/eZx4b\n+b3mgH62ieuiN89sHgfMZePWZKnGtfTmdb9i+ugvAR9jPJTk7C8B3zS7Pxzi0V9cvvjzK062nNWM\nMX41yZ/uunq/2q9L8voxxkNj55e235Od9ZuFfeaS7P3LkmY7lzHG748xbl/8/S+T3Jmd3ze4cety\nwFySzVuXv1r89UnZ+U/7n2YD1yTZdy7Jhq1JklTVFUmuTfLfcq7+jVyXmdGbG21Lb96WvpzozZnv\nXPTmGc5lnb153cF0r1/wffk+952rkeQtVfXOqvqWxXWXjXOfbvhgkst6SjuS/Wp/enbW56xNWatv\nrao7qurmpbcNbMRcaucTNZ+f5O3Z8HVZmsvbFldt1LpU1ROq6vbsPPe/PMb4rWzomuwzl2TD1mTh\nh5N8V5JHlq7byHWZGb15frZpX2/i95pH6c3zmYve/KhZzSVr7M3rDqbb8MlKXzjGeH6Sa5L8q6r6\n4uUbx87r1Bs5zxVqn/u8XpPkmUmel+QDSX7wgPvOai5V9TFJfi7Jt48x/mL5tk1bl8VcfjY7c/nL\nbOC6jDEeGWM8L8kVSb6kqq7edfvGrMkec7kqG7gmVfUPk/zBGOO27H1GeaPWZWa24XnRm+dp477X\nLNObHzWLuejNjx1i7UWuYN29ed3B9P4kp5Yun8pjU/PsjTE+sPjzD5P8z+y8/PxgVX18klTVJyT5\ng74Kz9t+te9eqysW183WGOMPxkJ23k5w9q0Bs55LVT0xO43vp8YYP7+4eiPXZWkuP312Lpu6Lkky\nxvhgkl9M8oJs6JqctTSXz9nQNfmCJC+uqt9L8vokf6+qfiobvi4zoTfPz1bs6w39XpNEb156+Kzm\nkujNmddc1tqb1x1MH/0l4FX1pOz8EvBb13zMyVTVJVV16eLvH53ky5O8Oztz+IbF3b4hyc/vPcIs\n7Vf7rUm+pqqeVFXPTPLsJO9oqG9li41/1ldmZ22SGc+lqirJzUl+e4zxI0s3bdy67DeXTVuXqnra\n2bfPVNVHJfmyJLdlM9dkz7mcbRYLs1+TJBljvHyMcWqM8cwkX5PkrWOMr88GrssM6c3zsxX7etO+\n/5+lN89vLnrzPOey9t481v+pTddk5xPB7knyPes+3sS1PzM7nyR1e5L3nK0/yVOTvCXJ7yR5c5Kn\ndNe6T/2vT/JAkr/Ozs8TfeNBtSd5+WKd7kryou76D5nLNyX5ySTvSnLH4h/AZXOfS3Y+he2RxZ66\nbfF1ehPXZZ+5XLNp65Lk7yT5v4t5vCvJdy2u38Q12W8uG7Ume8zrS3Puk/82bl3m+BW9ubP+rejN\ne8xjI/vyoja9eWZzOaCfbeKa6M0rzqUWDwAAAIAW634rLwAAABxIMAUAAKCVYAoAAEArwRQAAIBW\ngikAAACtBFMAAABaCaYAAAC0EkwBAABo9f8BPOablArDCpgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f18912992d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (subfig1,subfig2) = subplots(1,2,figsize=(16,7)) # one figure with two horizontal subfigures\n",
    "subfig1.stem(xsharp)\n",
    "subfig1.set_ylim(0,1.1)\n",
    "subfig2.stem(x_restored)\n",
    "subfig2.set_ylim(0,1.1)\n",
    "subfig1.set_title('$x^\\sharp$')\n",
    "subfig2.set_title('$x_\\mathrm{restored}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the original signal is highly sparse, it is perfectly recovered.\n",
    "\n",
    "We display the convergence speed of the $\\ell^1$ norm on the first half iterations, in log\n",
    "scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/gorthaur/lib/python2.7/site-packages/IPython/kernel/__main__.py:1: RuntimeWarning: divide by zero encountered in log10\n",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f188e1551d0>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh4AAAFwCAYAAAD31XL9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm4nWV57/HvTQIKiCJKGQNJERBQxjI4gEFBARkdoUUU\n0FY9VwG1hYI9lFNFj1rRepRiDTNIkTkMAmGIVqERCQEEIoPMMwQQAQ3Dff5435CduLOzs9c7rOH7\nua51Ze01PO/NQ3D/fKY3MhNJkqQmLNV2AZIkaXAYPCRJUmMMHpIkqTEGD0mS1BiDhyRJaozBQ5Ik\nNabj4BEREyLi6oi4JSJ+ExEHVVGYJEnqP9HpOR4RsSqwambOiojXAdcDe2bmbVUUKEmS+kfHIx6Z\n+Uhmziqf/wG4DVi903YlSVL/qXSNR0RMBDYDZlTZriRJ6g+VBY9ymuVs4OBy5EOSJGkB46toJCKW\nBs4BTsvM8xd6z5vBSJLURzIzxvrdKhaXBnAy8GRmfmGY97OTAvtdRByVmUe1XUe3sn8Wzb4Zmf0z\nMvtnZPbPonX6e72KqZZ3AfsC20fEDeVjpwralSRJfabjqZbM/AUeRCZJkkbBwNC+6W0X0OWmt11A\nF5vedgFdbnrbBXS56W0X0OWmt11Av+p4jcdiL+AaD0mS+kY3rPGQJEkaFYOHJElqjMFDkiQ1xuAh\nSZIaY/CQJEmNMXhIkqTGGDwkSVJjDB6SJKkxBg9JktQYg4ckSWqMwUOSJDXG4CFJkhpj8JAkSY0x\neEiSpMYYPCRJUmMMHpIkqTEGD0mS1BiDhyRJaozBQ5IkNcbgIUmSGmPwkCRJjTF4SJKkxhg8JElS\nYwwekiSpMQYPSZLUGIOHJElqjMFDkiQ1xuAhSZIaY/CQJEmNMXhIkqTGGDwkSVJjDB6SJKkxBg9J\nktQYg4ckSWqMwUOSJDXG4CFJkhpj8JAkSY0xeEiSpMYYPCRJUmMMHpIkqTEGD0mS1JhGgkcEhzdx\nHUmS1N0iM+u9QERCPgFsm8nsWi8mSZJqFRGZmTHW7zc11fKvwA8jnNqRJGmQNRUEjgWWBfZv6HqS\nJKkLdRw8IuKEiHg0Im5e1GcyeRn4W+DrEazS6TUlSVJvqmLE40Rgp8V9KJNZwPHAlAjGPDckSZJ6\nV8fBIzP/G3hqlB//F2A14LOdXleSJPWeRhd7ZjIX+BvgXyPYsMlrS5Kk9jW+yyST3wJHAD+O4LVN\nX1+SJLVnfBMXiYijhvw4HXIKsAPwHxEckEm9h4lIkqQxiYjJwOTK2qviALGImAhcmJlvH+a9YQ8a\niWB54JfAyZl8p+MiJElS7Vo/QCwizgCuAdaLiPsjYlRndWTyHLAH8I8R7NZpHZIkqfs1cmT6SMko\ngi2Bi4FPZHJZrcVIkqSOtD7i0alMrgP2BE6NWPx5IJIkqXe1HjwAMrkG2As4MYKDPWBMkqT+1PpU\ny4KfZRJwFvAk8JlM7qu1OEmStER6fqplqEzuBt4B/ByYGcERESzXclmSJKkiXRU8ADJ5MZOjKQLI\npsDtERwQwbiWS5MkSR3qqqmW4b/PNsA3gVWAbwOnZvJCVfVJkqTR6/j3ercHj6INAngP8CWKkZBz\ngB8D/53JK51XKUmSRmMggseC7bEWsDewD7AyxWLUc4FrMnm5qutIkqQ/N3DBY8G22RD4CPAhYDXg\nfIrRkKszebGOa0qSNMgGOngseB3eQnEWyEeAtYCTgeMzuaPua0uSNCgMHsNekw2AA4H9gFuBKcA5\nLkqVJKkzBo8Rr80ywG7Ap4GtgNOBH2ZySxv1SJLU6/rqALGqZTI3k3My2RnYHHgauDyCX0awXwTL\ntlyiJEkDpa9HPIYTwXhgV+BvWXAU5NZWC5MkqQc41dKBCCZSrAU5EPgdcCpwViZz2qxLkqRuZfCo\nQARLAzsDfwPsBPwMOBO4JJOn2qxNkqRuYvCoWAQrMH9b7mTgOuACYGom97RXmSRJ7TN41Ki8M+6O\nwB7ALsAfKe6cO+9xRyb1dqAkSV3E4NGQ8n4x6wHblY/3AMsA1wDXAzOBmZk82lqRkiTVzODRogjW\nBt4JbEaxXXdz4AXKEALcANwI3OvN7CRJ/cDg0UXKUZG1mR9CNgc2BlYEbqM4RfVW4BbgTopA4mmq\nkqSeYfDoARG8AdgA2AjYsPzzLynuKfM0cE/5uHvI83sogskfGy5XkqRFMnj0sAiWAlYFJpaPSUOe\nTwQmAHMopmtmMn8tyb0uapUktcHg0cciGAesCWxCMW2zRfnnssCvgV+UjxmZPNdWnZKkwWHwGEAR\nrApsCby7fGxKsW7kCuAy4NpM5rZXoSSpXxk8RHmzu60ozhz5AMW2359RhJDLgTudmpEkVcHgoT8T\nwZuBHYD3UwSRPwGXABcB012wKkkaK4OHRlRu8d2I4uTVXSnWi0ynCCEXZ/JQe9VJknqNwUNLJII3\nUYyC7EpxQ7y7KULIRcD1HnQmSRqJwUNjFsF4ipNXdy0fKwEXU4SQKzJ5tsXyJEldyOChykSwDvBB\nihDyDuBaytGQTH7XZm2SpO5g8FAtIliBYnHqB8vHk8yfkrkmk5daLE+S1BKDh2pXnrD6V8yfklmb\nYqvuRcClmcxpsTxJUoMMHmpcBGswf0pmMjALuBC4IJPbWyxNklQzg4daVR5eNhnYA9id4qZ35wMX\nANe5S0aS+ovBQ12jnJLZkiKE7AmsCEylCCJXZ/KnFsuTJFXA4KGuFcF6zA8hG1GsCzkf+GkmT7dZ\nmyRpbAwe6gkRrALsRhFE3kOxVfd04NxM/tBmbZKk0TN4qOdE8DqKxan7AttSHFp2GjDNbbqS1N0M\nHuppEawMfJwihEwE/gs4FZjpHXUlqfsYPNQ3IliXIoDsC8wFTgFOyeTBVguTJL3K4KG+U95R9x3A\nJ4GPAjOAE4GpmfyxzdokadAZPNTXIlgO2AvYH9gUOJMihFzvVIwkNc/goYERwdrAfsCngOeBH1FM\nxbg1V5IaYvDQwCkPKtsO+DtgJ+As4NhMZrVamCQNAIOHBloEqwIHUoSQB4BjgbNdCyJJ9TB4SEAE\n4ynOBvkcsDnFOpDjMrm71cIkqc90+nt9qQoK2CkiZkfEHRFxWKftSWORyUuZXJDJTsC7gHHAdRFc\nHMEHIxjXcomSJDoc8YiIccBvgR2AB4HrgH0y87Yhn3HEQ60o75z7ceDzwMrAccAJmTzeamGS1MPa\nHvHYCrgzM+/JzBcpTp3co8M2pUpk8kImJ2WyFfAxYH3g9ghOi2DzlsuTpIHUafBYA7h/yM8PlK9J\nXSWT6zI5AFgHuAG4MIKpEWzRcmmSNFDGd/h9D3BST8lkDvDtCH4AfBq4IIKZwFGZzGy3unpF8Bpg\nR+C9wCbA6hRrYZ4HHgVuojgl9rJMnm2rTkn9rdPg8SAwYcjPEyhGPRYQEUcN+XF6Zk7v8LpSR8rt\ntt+PYApFALkogunAl/ttJ0wEKwJfoFjrchvwU+CbFP+tvggsRzFSuQnF1uQTyr74AXC5J8RKgy0i\nJgOTK2uvw8Wl4ykWl74PeAj4FS4uVQ+KYHngi8AhFDen+2omT7ZbVeci2Av4PnA58LVM7hjFd14P\nfJiiL8YD/wycbwCRBF1wjkdE7Ax8l2LI9vjM/HqVBUpNimAV4EiK3TDfBL6bydx2q1py5fbhr1Es\nqt03k1+OoY0APkDRD88An8vkN5UWKqnntB48FnsBg4d6UATrUQTqtYDPZvKLlksatfJI+eMpFtJ+\nKJMnOmxvHMV01FeBfwP+LZOXOy5UUk8yeEg1Kf8f/4cpAsilwKHl4tSuVi6c3Qj4YCbPVdjuROAk\nYC6wdy/0haTqtX2Oh9S3MslMzgY2pNj5cXMEu7dc1ogi+CzFIrDdqwwdAJncQ3FY4E3AryLYqMr2\nJQ0GRzykUYpgO+AE4Frg4G77f/zloWiXAu/M5M6ar7UfxdqPPTKZUee1JHUXRzykhmTyc4otp3Mo\nRj92bbmkV0XwWuBU4JC6QwdAJqcwfxvy9nVfT1L/cMRDGoNy9ONU4Gzg8LZ3vkTwVWAD4CNNbnuN\nYDJwFsXUzrVNXVdSe1xcKrUkgjcBJwKrAB8v10C0Ucck4NfAxpk82ML1d6Hoh/e53Vbqf061SC0p\nDxjbAzgTmBHR2g0SvwV8p43QAZDJJRSHjV0SwWpt1CCpdzjiIVUggq2Bn1Ccn/GVpqY7ypvcTQXe\nkskLTVxzhFqOBHYCts/kT23WIqk+TrVIXSKCVYELgLuAA8r7wdR9zQuAKzP5Xt3XGkUtS1GEr98D\nB3rEutSfnGqRukQmj1CcobEUML0MIrWJYDNgS+BHdV5ntDJ5BfgksA2wb8vlSOpSBg+pQuV0xz4U\n52nMqPmQrSOBb7Q9xTJUeWjZPsAxEazTdj2Suo9TLVJNItgX+DawZ9VbTSPYFLgEWKebgsc8ERxM\nEUDenclLbdcjqTpOtUhdKpPTgP2BqRF8oOLmjwS+1Y2ho/Q94FmK3S6S9CpHPKSaRfAu4FzgoEzO\nrKC9jYHLKEY7nu+0vbqUUy0zgC0zubvteiRVwxEPqctl8ktgR4p1D5+poMkjKW5N37WhAyCTuyim\nmo4t7/QrSY54SE2JYF1gGvC9TI4ZYxsbA5dTjHZUevfZOkSwNDAT+JdMzm27Hkmd8xwPqYdEMAG4\nAjgD+D9LetZFBBcB0zL59zrqq0MEOwLHAhu1fU8bSZ1zqkXqIZncD2wH7Al8uzx0a1Qi2BZ4G3Bc\nTeXVIpNpFIeq/V3btUhqnyMeUgsieCNwMcUv5AMXNxJQBpRrgR+Ut6TvKRG8nWKkZ/1Mnm67Hklj\n54iH1IMyeQrYAXg9cFEEKyzmK/sDrwCn1V1bHTK5mSJofantWiS1yxEPqUURjAd+QHH0+S7lsesL\nf2Y14Iby/ZkNl1iZCP4S+BXFDe0c9ZB6lCMeUg8rT/X8LHAecE0Ebxv6fgTLAxcB/6+XQwdAJr+j\nGPX4+7ZrkdQeRzykLhHBfhTnXkwBzgZWBv4FuI0+udtrBOsDv6DYDvz7tuuRtOTcTiv1kQjWAL4M\nvBt4CjgdmFLe+bUvRHAGMCuTb7Rdi6QlZ/CQ1FMi2Ay4EJiUyYtt1yNpybjGQ1JPyeQG4HfAh9qu\nRVLzDB6S2vDvwMFtFyGpeQYPSW24AFgjgi3bLkRSswwekhpXbiP+PnBQ27VIapaLSyW1IoI3URwZ\nP6k8yVVSD3BxqaSelMmTwGXAX7ddi6TmGDwktWkKcGDbRUhqjsFDUpuuBFYqz/aQNAAMHpJaU57I\neiKOekgDw8WlkloVwVoUd99dI5M/tl2PpJG5uFRST8vkPmAWsGvbtUiqn8FDUjc4Bdiv7SIk1c+p\nFkmti+B1wAPAepk81nY9khbNqRZJPS+TPwBTgX3arkVSvQwekrrFycAn2y5CUr0MHpK6xXRg5Qje\n3nYhkupj8JDUFTJ5GTgVF5lKfc3FpZK6RgRvBa4C1irvYCupy7i4VFLfyGQ2cD/wvrZrkVQPg4ek\nbnMKLjKV+pZTLZK6SgRvAu6imG75fdv1SFqQUy2S+komT1Ks8/ho27VIqt6Yg0dEfDQibomIlyNi\n8yqLkjTwTsbdLVJf6mTE42ZgL+DnFdUiSfP8FNgggkltFyKpWmMOHpk5OzNvr7IYSQLIZC7wXzjq\nIfUd13hI6lY/BD4bwWvbLkRSdcaP9GZETANWHeatIzLzwtFeJCKOGvLj9MycPtrvShpMmdwSwQ3A\nvsCUOq4RwdLA+sAE4E3ASsCKzP/fxjnAo8DDwG+BRzKpdyug1GUiYjIwubL2Ot1OGxFXA1/KzJmL\neN/ttJLGJILtgf8ANiqPVO+0veWA9wK7AO8C1qU4sOxuipDxJPA08BIQFEFkFWB14K0UgeRm4Grg\nCmBGOS0kDYxOf6+POOKxJHVU1I4kDTUdeAw4EPjPsTQQwQRgN2BX4N3A9cDFwAnALZm8sARtrQxs\nShFejgHeEsHFwFnAZZn8cSw1SoNkzCMeEbEX8D3gzcAzwA2ZufMwn3PEQ9KYRbApcBmwQSZzRvmd\nFYG/BfYG1gIuAS4ELs/kmQprWw34EMWZI5uU1zkLuNQQon7V6e91Ty6V1PUiOBZYJpNPL+ZzawCH\nAAdQhIAfAdc0ccO5CFZlfgjZGDgDOCGTYaehpV5l8JDU9yJ4PTALODqT44d5fwPgHyjOFjoF+E4m\n9zZb5QL1TKS438z+wPPA2eXjZhenqtcZPCQNhAjWA6ZRTJmcCNxLsWZjf2Ab4PvAseWR610hgqWA\nrYEPAx8B5gLnUISQmYYQ9SKDh6SBEcGbgSOAHSi2wN4InA6cnsnzbda2OBEEsAXzQ8h4ihNary8f\nt2TyYnsVSqNj8JCkHlOGkI0pzkb4K4pAMhG4BbgduAP4HcUZIo+Vfz5uMFE3MHhIUh+IYAXg7RRn\ni6xLEUT+onysQrGD8Fnmh5HHFvP8WadyVAeDhyQNgHK9yBspQshfsGAoGe75MiwYRB4B7qNYG3Mv\ncA/wgKMoWlIGD0nSn4lgWRYMKKtTnGkyEVi7fKwKPEAxtXNn+eddFKe5Pgg86aiJFmbwkCSNSQTL\nUASRedM7bwHWAdYE1gCWBR6iCCF3UWxpngXMqvIgNvUWg4ckqRblvW3WoAgi61GczroZxVqUR4Cr\nKLY4X9VN25hVL4OHJKlREYwDNgTeB+wIbEuxI+dc4NxM7mqxPNXM4CFJalU5ZTOZ4oySPYGHKU6Q\nPS2Tx1osTTUweEiSukY5GrIt8CmKEHIVxZ2AL23injmqn8FDktSVynvsfIzipn0TgZOBH2TyQJt1\nqTOd/l5fqspiJEmaJ5PfZzIlk3cC76XYJXNTBKdFsFnL5aklBg9JUu0ymZ3JIcBfAjcBF0ZwVQST\n261MTXOqRZLUuHJB6j7AkRT3pfnnTGa0W5VGw6kWSVLPyWRuJicDbwV+ApwTwYkRrNxyaaqZwUOS\n1JpMXszkR8AGwFPALREcWN7BV33IqRZJUteIYFNgCsXN7Q7M5OGWS9JCnGqRJPWNTGYB7wB+BcyK\n4GMtl6SKOeIhSepKEWwJnAr8D/D5TJ5vuSThiIckqU9lch2wBTAOuDaCt7Rckipg8JAkda1MngP2\nA44Drolgj5ZLUoecapEk9YQItgbOAr6byTFt1zOovFeLJGlgRDABuBT4KXBoJq+0XNLAcY2HJGlg\nZHI/xd1vtwZOiWB8yyVpCRk8JEk9JZM5wPuBNwMnRTCu5ZK0BAwekqSek8kLwF7AasCUCH+f9Qr/\nRUmSelIZPnYH1gO+0nI5GiWDhySpZ5XbbfcE9o7gUy2Xo1FwUY4kqadl8ngEHwR+FsE9mUxvuyYt\nmiMekqSel8lsYB/gzAjWbrseLZrBQ5LUFzK5CvgW8JMIlmm7Hg3PA8QkSX0jggAuAO7K5Att19OP\nPEBMkqRSJgl8CvhQBDu3XI6GYfCQJPWV8oCxA4D/jOCNbdejBTnVIknqSxH8AFg+0222VXKqRZKk\n4R0GbBfBbm0XovkMHpKkvpTJH4D9geMieEPb9ajgVIskqa9FMAV41l0u1ej097rBQ5LU1yJYGbgV\neG8mN7ddT69zjYckSSPI5HHgKOD75TkfapHBQ5I0CI4DVgT2aLuQQWfwkCT1vUxeBv4Z+EoE49qu\nZ5AZPCRJg+Ii4A/A3m0XMshcXCpJGhgRbA/8CNggkxfbrqcXtbq4NCK+FRG3RcSNEXFuRLhPWpLU\ntTK5GrgX+ETbtQyqTqdaLgc2ysxNgNuBwzsvSZKkWh0NHOpaj3Z0FDwyc1pmvlL+OANYs/OSJEmq\n1dXAM8CebRcyiKpcXHoAcEmF7UmSVLlMEvg68E+e69G8xQaPiJgWETcP89htyGe+DMzNzB/XWq0k\nSdWYCiwPvK/tQgbN+MV9IDN3HOn9iPgUsAsj/MuLiKOG/Dg9M6ePrjxJkqqXySsRHAMcDFzRdj3d\nLCImA5Mra6+T7bQRsRPwbeA9mfnEIj7jdlpJUteJYDngPmDLTO5uu55e0epN4iLiDmAZYE750rWZ\n+fkqC5QkqS4R/BvwSiaHtl1Lr/DutJIkjVEE6wD/A6ydyfNt17OwCF4DbApMotg5+lrgOeBu4Dbg\n9nKxbIM1GTwkSRqzCC4Czs3khLZrAYjgDRTHun8c2Aq4o3w8ALwArEARRDYGxgGXAScBv2gihBg8\nJEnqQAS7Akdk8s6W61gFOAL4JDANOBWYnsnvF/H5AN4C7A58BngS+IdMrq23ToOHJEljFsHSwP3A\ndpnc3lINewPfA84Avp7JI0v4/XHAvsBXKQ5IOyiTpysvlJbv1SJJUq8rbxZ3OsVIQ6MiWDaC04Cj\ngF0yOXhJQwdAJi9ncjLwVoo78N4YwXbVVlsNg4ckSXAy8Ikm799STq1cDQSweSa/7rTNTJ7L5PPA\n54CzIvh0p21WzeAhSRp4mdwEPAFs38T1IlgbuAa4FNi36h01mVwCbEtxLPyXqmy7UwYPSZIKJ9HA\ndEsEqwNXAv+eyVF17UQp16tMBj4XwRfquMZYuLhUkiQggtWAW4BVM5lb0zWWA34OnJ/JV+u4xjDX\nXAuYDhydyfGdt+fiUkmSOpbJw8BsKrwvyVDl9tfjy2scXcc1hpPJfcDOwNcieG9T110Ug4ckSfOd\nD+xVU9uHAusAn2n6tNFMfgvsA5xeTvW0xqkWSZJKEawL/AxYM5NXKmx3a2AqsEUmD1TV7hjqOJJi\nRGfHTF4eWxtOtUiSVIlM7gCeArasqs0IVqA4J+TzbYaO0tHA0sBn2yrAEQ9JkoaI4GhgXCb/VFF7\nxwOZ2R1nakSwIcWoziaZPLTk33fEQ5KkKp0H7FlFQ+Xpoe+H7tnOmsmtwHEUR7Q3zuAhSdKCrgdW\niGD9Thop7wFzLPDFTJ6tpLLqfA3YJIKdmr6wwUOSpCHKHSdTKe762okDgUeBszsuqmKZvAAcBnwj\notksYPCQJOnPXQDsMdYvlweF/W/gsKa3zi6B84AXgL9u8qIuLpUkaSERvIZitGLdTB4fw/cPo9g6\n+7HKi6tQBNsDPwQ2GO32WheXSpJUsUz+BEwDdl3S75ajHV+kuNV9t5sOPE3n00qjZvCQJGl4FzC2\nX8ifBn5Z7h7pauU00DeBf2jqmk61SJI0jAhWAu6muGncC6P8zjLAncCHM7muzvqqEsE44A7g46Op\n2akWSZJqkMkcYCbwviX42t8Av+2V0AFQru2YAs0ccGbwkCRp0aYyyt0t5bkdh1OckdFrTgY+GsHy\ndV/I4CFJ0qJNBXYb5VkXBwL3ZnJ1zTVVLpMHgWuAj9Z9LYOHJEmLkMldwBPAViN9rhwpOBKqub9L\nS04G9qn7IgYPSZJGdg6LP2TrEODnmVzfQD11+SmwTQRvrPMi7mqRJGkEEawJ3ARMyuSZYd5/MzAb\n2CaTO5uur0oRnAecl8kpi/6Mu1okSapNJg8AlwEHLOIjRwJn9nroKJ0LfLjOCzjiIUnSYkSwNXAm\nsH55qum817cBzgfePpaj1btNBCsC9wOrZPL88J9xxEOSpFplMgP4DXDQvNfKBaXHAwf1Q+gAyORp\n4EbgXXVdY3xdDUuS1GcOAa6JYDbwCPAtYAZwVqtVVe9KikPTptXRuMFDkqRRyOTOCPYCTgJeohjt\n+G4X3/Z+rK4Ejqmrcdd4SJKkV5X3m3kCmFgeG7/Q+67xkCRJFclkLvBLYHId7Rs8JEnSwn5OTQtM\nDR6SJGlh1wFb1tGwazwkSdICymPT7wNWzOTlBd9zjYckSapQJk9RbBneoOq2DR6SJGk4tUy3GDwk\nSdJwDB6SJKkxtQQPF5dKkqQ/E8EKFOs8Vsjklfmvu7hUkiRVLJNngTnA2lW2a/CQJEmLchsV72wx\neEiSpEUxeEiSpMbMBt5aZYMGD0mStCi/AyZV2eCYg0dEfCUiboyIWRFxZURMqLIwSZLUurupOHiM\neTttRKyQmc+Wz/8e2CQzPz3M59xOK0lSD4rgNcDvgeXm3bOlte2080JH6XXAE2NtS5IkdZ9M/gQ8\nDqxZVZvjO/lyRBwNfAJ4HtimkookSVI3mTfdcm8VjY044hER0yLi5mEeuwFk5pczcy3gJOA7VRQk\nSZK6yn1AZes4RxzxyMwdR9nOj4FLFvVmRBw15MfpmTl9lO1KkqRWnTgOpu4Xcf46VbTWyeLSdTPz\njvL53wNbZeYnhvmci0slSepREXwBWDuTQ4qfO/u93skaj69HxPrAy8BdwOc6aEuSJHWnh4B3VNXY\nmINHZn6kqiIkSVLXehhYvarGPLlUkiSN5CEMHpIkqSEPA6tFUMl6TYOHJElapEyeA14CXl9FewYP\nSZK0OE8CK1XRkMFDkiQtzhwMHpIkqSGOeEiSpMbMAd5URUMGD0mStDhOtUiSpMY8iSMekiSpIY54\nSJKkxri4VJIkNcbFpZIkqTFPA2+ooiGDhyRJWpzngOWraMjgIUmSFsfgIUmSGmPwkCRJjXkeg4ck\nSWqIIx6SJKkxc4GlIli604YMHpIkaUSZJBWNehg8JEnSaBg8JElSYwwekiSpMc8By3XaiMFDkiSN\nhiMekiSpMQYPSZLUGIOHJElqjMFDkiQ1xuAhSZIaU8n9WsZXUIgkSep/5wEvdNpIZGYFtYxwgYjM\nzKj1IpIkqRGd/l53qkWSJDXG4CFJkhpj8JAkSY0xeEiSpMYYPCRJUmMMHpIkqTEGD0mS1BiDhyRJ\naozBQ5IkNcbgIUmSGmPwkCRJjTF4SJKkxhg8JElSYwwekiSpMQYPSZLUGIOHJElqTMfBIyK+FBGv\nRMRKVRQkSZL6V0fBIyImADsC91ZTzuCJiMlt19DN7J9Fs29GZv+MzP4Zmf1Tn05HPI4BDq2ikAE2\nue0CutzktgvoYpPbLqDLTW67gC43ue0CutzktgvoV2MOHhGxB/BAZt5UYT2SJKmPjR/pzYiYBqw6\nzFtfBg4H3j/04xXWJUmS+lBk5pJ/KeJtwJXA8+VLawIPAltl5mMLfXbJLyBJkrpWZo55sGFMwePP\nGom4G9giM+d03JgkSepbVZ3j4aiGJElarEpGPCRJkkaj1pNLI2KniJgdEXdExGF1XqsbRcQJEfFo\nRNw85LVgG/Z9AAAEDUlEQVSVImJaRNweEZdHxIpD3ju87KvZEfH+4VvtHxExISKujohbIuI3EXFQ\n+bp9BETEayNiRkTMiohbI+Lr5ev2TykixkXEDRFxYfmzfVOKiHsi4qayf35Vvmb/lCJixYg4OyJu\nK//72tr+KUTE+uXfm3mPZyLioMr6JzNreQDjgDuBicDSwCxgg7qu140PYFtgM+DmIa99Ezi0fH4Y\n8H/L5xuWfbR02Wd3Aku1/c9Qc/+sCmxaPn8d8FtgA/togT5arvxzPPA/wLvtnwX654vA6cDU8mf7\nZn7f3A2stNBr9s/8vjgZOKB8Ph54g/0zbD8tBTwMTKiqf+oc8dgKuDMz78nMF4H/Avao8XpdJzP/\nG3hqoZd3p/gLT/nnnuXzPYAzMvPFzLyH4l/cVk3U2ZbMfCQzZ5XP/wDcBqyBffSqzJy3c2wZijD/\nFPYPABGxJrALMIX52/ntmwUtvPPA/gEi4g3Atpl5AkBmvpSZz2D/DGcHit/l91NR/9QZPNYA7h/y\n8wPla4Nulcx8tHz+KLBK+Xx1ij6aZ6D6KyImUowOzcA+elVELBURsyj64erMvAX7Z57vAP8IvDLk\nNftmvgSuiIhfR8Rnytfsn8Ik4PGIODEiZkbEjyJieeyf4ewNnFE+r6R/6gwerlpdjCzGqEbqp4Ho\nw4h4HXAOcHBmPjv0vUHvo8x8JTM3pTgrZ7uI2H6h9weyfyJiV+CxzLyBRRxeOKh9M8S7MnMzYGfg\nf0XEtkPfHPD+GQ9sDhybmZsDzwH/NPQDA94/AETEMsBuwFkLv9dJ/9QZPB6kmBOaZwILJqJB9WhE\nrAoQEasB8w5cW7i/5h3K1tciYmmK0HFqZp5fvmwfLaQcBr4Y2AL7B+CdwO7lGUJnAO+NiFOxb16V\nmQ+Xfz4OnEcx9G3/FB6guOXHdeXPZ1MEkUfsnwXsDFxf/h2Civ7+1Bk8fg2sGxETy9T0cWBqjdfr\nFVOBT5bPPwmcP+T1vSNimYiYBKwL/KqF+hoTEQEcD9yamd8d8pZ9BETEm+etGo+IZSnuBH0D9g+Z\neURmTsjMSRRDwVdl5iewbwCIiOUiYoXy+fIUt7e4GfsHKNaXAfdHxHrlSzsAtwAXYv8MtQ/zp1mg\nqr8/Na+G3Zlip8KdwOFtrsxt41H+C3sImEux3mV/YCXgCuB24HJgxSGfP6Lsq9nAB9quv4H+eTfF\n/Pwsil+oNwA72Uev/rO+HZhZ9s9NwD+Wr9s/C/bTe5i/q8W+Kf5ZJ5V/b2YBv5n3v7/2zwJ9tAlw\nHXAjcC7Frhb7Z/4/7/LAE8AKQ16rpH88QEySJDWm1gPEJEmShjJ4SJKkxhg8JElSYwwekiSpMQYP\nSZLUGIOHJElqjMFDkiQ1xuAhSZIa8/8BsTIYYIpnsyUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f18912f2090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(log10(En_array-En_array.min()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The convergence is linear in practice. Convergence up to machine precision (1e-14 here) is achieved after a few hundreds of iterations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>\n",
    " <font color='magenta'>\n",
    "__Exercise: test the recovery of a less sparse signal.\n",
    "What do you observe?__</font></h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "S = 31\n",
    "random.seed(0)\n",
    "sel = random.permutation(N)\n",
    "sel = sel[0:S]   # indices of the nonzero elements of xsharp\n",
    "xsharp = zeros(N)\n",
    "xsharp[sel] = 1\n",
    "\n",
    "y = A.dot(xsharp)\n",
    "\n",
    "# put your code here\n",
    "s = zeros(N)\n",
    "En_array = zeros(nbiter)\n",
    "skm1=s\n",
    "for iter in range(nbiter):  # iter goes from 0 to nbiter-1\n",
    "    # put your code here\n",
    "    x=prox_f( skm1, y )\n",
    "    skm1=s\n",
    "    s = skm1+gamma*(prox_gamma_g(2*x-skm1,gamma)-x)\n",
    "    \n",
    "    En_array[iter] = norm(x, ord=1)  \n",
    "x_restored = x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f188cf5b6d0>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6YAAAG0CAYAAAAhAzpSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xu0HeV55/nfowtIQkowTcL1yAQJY8NxMvEFI9wI0bGR\ndHxrexJsPNOxMZLpOEj0mqwJdsBBTuPpdlbPTCR8Y2Fw4+7pAGvsZNnD0cVu+1iejHBMQjCywUbY\nNIdrg/FFMhJI6Jk/apd37Tq171X7rar9/ax11tm7du23nnrr3Zdnv2+9Ze4uAAAAAABCmRc6AAAA\nAADAeCMxBQAAAAAERWIKAAAAAAiKxBQAAAAAEBSJKQAAAAAgKBJTAAAAAEBQJKYAAAAAgKAWhA4A\nwFxmNiHJJJ3m7ntCxwMAAAAUiR5ToJzmSTpV0ptDBwIAAAAUjcQUKKc3S7pB0i9DBwIAAAAUzdw9\ndAwAUszstyQdlXS6u/9d6HgAAACAItFjCpTT0dABAAAAAKNCYgqU10JJk6GDAAAAAIrGrLxAYGb2\ndkkvSbpQ0v2S1kl6SNI7Jf3XgKEBAFB5bT5nP+7uDwYNDEALzjEFAjKz5ZKOcfd9ZvYPkn5P0j+X\n9BNJj0k61d2/HTJGAACqqsPn7Nfd/fkRxXCOpAl33zmK7QFVRWIKlICZnSTpDndfEzoWAADqZtjP\nWTNb4e4PD/jc+ZL+s7tfNsjzgXHBOaZAQGb2SjP7HUlTknY3lk2FjQoAgHrI43O2MVP+OwaNwd1f\nknRo0OcD44JzTIGwLpG0TNKTkhaZ2b9UNIQXAAAML/Nz1syuUHTO6c8l/YmkFZJeK+llimbGv0/S\n+e7+f0g6T9LrGwnuLyWtUnTKzfON58XlfErS6xplfF3S4cb275W0ZBQ7C1QZiSkQkLtvCx0DAAB1\n1e5z1syekXSOu/9J4/6fSfq0oiT2NYp6WKcbq98t6ZXufp+ZfV7SBnd/yczukPSncTlm9gVFyeky\nSROS3ivpj9z9oJl9sLi9BOqBobwAAFSMmd1qZk+b2f1tHl9jZj83s3sbf9eNOkag5EzSLxL3XdLD\n7v41SX8n6f+SdHniMZnZuZKOk3RMY/lxjf9xOUcTZdzX2MZLRe0AUDdMfgQAQMWY2YWSDkj6gru/\nOuPxNZL+F3d/+6hjA6rAzP5Q0h9I+qC7P2lmKyT9vqJht69TlJye7O53mNlCSX8p6YuSnpZ0saSn\nFA3V/Y24HEXDdeMyfiBpkaS1ii5R8wlJl7r7I6PaR6BqSEwBAKggMztD0lc6JKZ/4u5vG3FYAAAM\nhKG8AADUj0u6wMzuM7PpxnUUAQAorZFNfmRmdM0CAHLl7hY6hpL6R0kT7v68ma2X9LeSXpFeic9m\nAEDeBv1sHumsvHX5AmFmW9x9S+g4hlWX/ZDYl7Kqy77UZT+k2u0LSVUb7r4/cXu7mX3azE5w9+cy\n1uWzuUTqsh8S+1JWddmXuuyHVLt9GfizmaG8AADUjJmdZGbWuH2eojkl5iSlAACUBdcxBQCgYszs\nryVdJOlEM5uVdL2khZLk7jcpmhn0j8zsiKTnJb0nVKwAAPSCxHQwM6EDyMlM6AByNBM6gBzNhA4g\nRzOhA8jJTOgAcjQTOgAMz90v6/L4pyR9akThlMVM6AByMhM6gBzNhA4gRzOhA8jRTOgAcjITOoAc\nzYQOoAxGdrkYM/O6nMcCAAiPz5XhUYcAgDwN87nCOaYAAAAAgKBITAEAAAAAQZGYAgAAAACCIjEF\nAAAAAARFYgoAAAAACIrEFAAAAAAQFIkpAAAAACAoElMAAAAAQFAkpgAAAACAoEhMAQAAAABBkZgC\nAAAAAIIiMQUAAAAABEViCgAAAAAIisQUAAAAABAUiSkAAAAAICgSUwAAAABAUCSmAAAAAICgSEwB\nAAAAAEGRmAIAAAAAgiIxBQAAAAAERWIKAAAAAAiKxBQAAAAAEBSJKQAAAAAgKBJTAAAAAEBQJKYA\nAAAAgKC6JqZmdquZPW1m93dYZ5uZPWRm95nZ7+YbIgAAAACgznrpMf28pHXtHjSzKUkr3f0sSR+U\n9JmcYgMAAAAAjIEF3VZw92+Z2RkdVnm7pNsa637bzI43s5Pc/elBAjKbnJImNkvLFklPLZUWmXT8\nfmn/IWl2m/ve6UHKbb+N/Mptln3CX0i/dqZ0rKT9P5ae+Ogg5RcZ5yDl97N+t3Wbj790qjT/FGnx\nk9LBJ+L1hqnHYeot67nRI+1j7aXc7ttq39Z7r8ti2kn/+zJ8nfe673nEN+w2OreZ0RyTLm1pjzSx\nKmT7AAAA6Mrdu/5JOkPS/W0e+4qkCxL3vybptRnrefftnDslbXhIcpe+6dKfeXQ7/tvwkHTuVC8x\n97aN/Mptlv3Ox+fGffkT/ZZfZJyDlN/P+t3WbT7e9hhfP2g9DlNv2c995+PRdvNtj7229d7rsph2\nUnRb7fTcvPatyG10bjOjOSad29I3XdrwYsj2UWzbk4eOoep/1CF//PHHH395/g3zudLrBrolpm9M\n3P+apNcMEqS0bkfzy9O1qS/q8d/67cNVVnIb+ZXbLDufuIuMc5Dy+1m/27rNx9vV1bpnBq3HYeot\n+7nXeudYBzsevbb13uuymHZSdFvt9Ny89q3IbXRuM6M5Jp3b0mhjGfUfSRV1yB9//PHHX7n+hvlc\n6TqUtwePS5pI3D+9sWwOM9uSuDvj7jOtayxb1Lx9Q5vNLV3cd4Rtt5FnuXHZ7aq03/KLjHOQ8vtZ\nv9u68eNt62rh4PU4TL1lPXdB6v8g5Xbb1iBtPV2X/Tw3b3nXebfnjuK11Os2OrWZQcvsV6e2NOpY\nimVmayStCRwGAAAoQB6J6ZclXSXpdjM7X9LPvM35pe6+pXNR+w9139yBg/0G2Ns2hi03LvtIm8f6\nLb/IOAcpv5/1u60bP962rg4PXo/D1FvWc4+k/g9Sbq/byirbrfN2i24nvci7zuPndtv3XhW5jU5t\nZtAy+9WpLY06lmI1fsycie+b2fXBggEAALnq5XIxfy3p/5N0tpnNmtkHzOxKM7tSktx9WtKPzGyf\npJskfWjwcGa3SRv3Ne9fm3p8w8PSozcOXn7WNvIqNy77+0/MjfvyJ/svv8g4Bym/n/W7rRs/fknG\ndjY8LM1+cvB6HKbesp77vSelDzzZPtZBj0evbb3XuswrrkHkXee97nse8Q27jU5tZti4e9WpLV0i\naePh0cUCAAAwoLKNN44m8li/PToP6sLvSG++p3k+WJ4T/8TbyK/cZtkXfqd5Lteb7xlugpxi4hyk\n/H7W77Zu6+PrnkmvN0w9DlNvWc/tFms+9d++rfdXl/m3kxB1nve+FbmN7m2m+GPSpS1dH7p9FLff\n8tAxVP2POuSPP/744y/Pv2E+V6xRQOHMzN3bDZvLWj+KL30735iKKTcuW5LyKL/IOAcpv5/1u60b\nP95uvWHqcZh6y3put1gH1Wtb77Uu84prEHnXeR7ljmobndrMoGUOE0N622VoH3nr93MFc1GHAIA8\nDfO50nUoLwAAAAAARSIxBQAAAAAERWIKAAAAAAiKxBQAAAAAEBSJKQAAAAAgKBJTAAAAAEBQJKYA\nAAAAgKBITAEAAAAAQZGYAgAAAACCIjEFAAAAAARFYgoAAAAACIrEFAAAAAAQFIkpAAAAACAoElMA\nAAAAQFAkpgAAAACAoEhMAQAAAABBkZgCAAAAAIIiMQUAAAAABEViCgAAAAAIisQUAAAAABAUiSkA\nAAAAICgSUwAAAABAUAtCBwAAAFAks8kpaWKztGyRtP+QNLvNfe906LiqiLoEUBQSUwAAUFtRIrVq\nq3TzyubSjSvMJkVC1R/qEkCRGMoLAABqbGJzayIlRfeXbwoTT5VRlwCKQ2IKAABqbNmi7OVLF482\njjqgLgEUh8QUAADU2P5D2csPHBxtHHVAXQIoDokpAACosdlt0sZ9rcs2PCw9emOYeKqMugRQHBJT\nAABQW9GkPHuulqZ2REumdkh3b2aynv5RlwCKZO4+mg2Zubtb7+vL3WXp2/nGVEy5cdmSlEf5RcY5\nSPn9rN9t3fjxdusNU4/D1FvWc7vFOqhe23qvdZlXXIPIu87zKHdU2+jUZgYtc5gY0tsuQ/vIW7+f\nK5hr3Oqwjq+DUKhLAFmG+VyhxxQAAAAAEBSJKQAAAAAgKBJTAAAAAEBQJKYAAFSMmd1qZk+b2f0d\n1tlmZg+Z2X1m9rujjA8AgH6RmAIAUD2fl7Su3YNmNiVppbufJemDkj4zqsAAABgEiSkAABXj7t+S\n9NMOq7xd0m2Ndb8t6XgzO2kUsQEAMAgSUwAA6uc0SbOJ+49JOj1QLAAAdLUgdAAAAKAQ6evIZV64\n3My2JO7OuPtMUQEBAOrFzNZIWpNHWSSmAADUz+OSJhL3T28sm8Pdt4wiIABA/TR+zJyJ75vZ9YOW\nxVBeAADq58uS/lCSzOx8ST9z96fDhgQAQHv0mAIAUDFm9teSLpJ0opnNSrpe0kJJcveb3H3azKbM\nbJ+kX0q6PFy0AAB0R2IKAEDFuPtlPaxz1ShiAQAgDwzlBQAAAAAERWIKAAAAAAiKxBQAAAAAEBSJ\nKQAAAAAgKBJTAAAAAEBQJKYAAAAAgKBITAEAAAAAQZGYAgAAAACCIjEFAAAAAARFYgoAAAAACIrE\nFAAAAAAQFIkpAAAAACAoElMAAAAAQFAkpgAAAACAoEhMAQAAAABBkZgCAAAAAIIiMQUAAAAABEVi\nCgAAAAAIisQUAAAAABBU18TUzNaZ2YNm9pCZXZPx+IlmtsPM/snM9prZ+wuJFAAAAABQSx0TUzOb\nL+mTktZJOkfSZWb2qtRqV0m6193/B0lrJP3vZraggFgBAAAAADXULYE8T9I+d39EkszsdknvkPRA\nYp0nJf124/avSfqJux/JOU4AAFBjZpNT0sRmadkiaf8haXab+97p0HEBAEajW2J6mqTZxP3HJL0h\ntc7Nkr5uZk9IWibp0vzCAwAAoRWdNEblr9oq3byyuXTjCrNJkZwCwHjolph6D2X8maR/cvc1ZrZC\n0lfN7Hfcff/w4QEAgJBGkzRObG4tX4ruT22SRGKKkaHnHginW2L6uKSJxP0JRb2mSRdI+rgkufvD\nZvZjSWdLuiddmJltSdydcfeZPuMFAIwpM1ujaC4DjNQoksZli7KXL12cT/lAd/TcA2F1S0zvkXSW\nmZ0h6QlJ75Z0WWqdByW9SdLfmdlJipLSH2UV5u5bhogVADDGGj9mzsT3zez6YMGMlVEkjfsPZS8/\ncDC/bQDd0HMPhNRxVt7GJEZXSdop6fuS7nD3B8zsSjO7srHa/ybpdWZ2n6SvSfpTd3+uyKABAMCo\njCJpnN0mbdzXumzDw9KjN+a3DaAbeu6BkLpe1sXdt0vanlp2U+L2s5Leln9oAAAgvNlt0sYVrT1J\n+SaN7nunzSYV9UxNr5OmdkiP3sjwSYwWPfdASObey/xGOWzIzN3del9f7i5L3843pmLKjcuWpDzK\nLzLOQcrvZ/1u68aPt1tvmHocpt6yntst1kH12tZ7rcu84hpE3nWeR7mj2kanNjNomcPEkN52GdpH\n3vr9XMFcvdZhdO7d8pEkjUV/NtftdRBKHesy+xzTDQ9Ld2/mRxKgN8N8NpOYkpiSmPb4XBLT7khM\nSUxHicR0eMN8NhcXE4lpFdS1Lkf5IwxQRySmA8dEYjpI+SSmJKZFxEBiOjwSU/SLxBSDqntd1n3/\ngKIM89nccfIjAAAAAACK1nXyIwAAAFRfNEx1YnM0++z+Q9LsNoapAigLElMAAICay57YZ+MKs0mR\nnNYTP0SgakhMAQAAam9ic2tSKkX3pzZJIlmpGX6IQBVxjikAAEDtLVuUvXzp4tHGgdFo90PE8k1h\n4gG6IzEFAACovf2HspcfODjaODAa/BCB6iExBQAAqL3ZbdLGfa3LNjwsPXpjmHhQLH6IQPWQmAIA\nANRcdF7hnqulqR3Rkqkd0t2bOd+wrvghAtVj7j6aDQ1xEe+iLnJc9EW8JSmP8ou+yHO/5fezfrd1\n48fbrTdMPQ5Tb1nP7RbroHpt673WZV5xDSLvOs+j3FFto1ObGbTMYWJIb7sM7SNvw1zEG5FhPpuL\ni6nYz+a6vQ4GkUc91L0u67B/0QRIyzdJ0+uiHyIevZEfIlC0YT6bSUxJTElMe3wuiWl3JKYkpqNE\nYjo8EtPxRGLaXZ32r077gvIb5rOZobwAAAAAgKBITAEAAAAAQZGYAgAAAACCIjEFAAAAAARFYgoA\nAAAACIrEFAAAAAAQFIkpAAAAACAoElMAAAAAQFALQgcAAAAAADGzySlpYrO0bJG0/5A0u81973To\nuFAsElMAAAAApRAlpau2SjevbC7duMJsUiSn9cZQXgAAAAAlMbG5NSmVovvLN4WJB6NCYgoAAACg\nJJYtyl6+dPFo48CokZgCAAAAKIn9h7KXHzg42jgwaiSmAAAAKDWzySmz9TvMLp2J/k9OhY4JRZnd\nJm3c17psw8PSozeGiQejwuRHAAAAKC0mwxkv7nunzSYlTW2SptdJUzukR2/kWNefuftoNmTm7m69\nry93l6Vv5xtTMeXGZUtSHuUXGecg5fezfrd148fbrTdMPQ5Tb1nP7RbroHpt673WZV5xDSLvOs+j\n3FFto1ObGbTMYWJIb7sM7SNv/X6uYK5hPpuLi6nYz+a6vQ4GkUc9jLIuzdbvkLavnfvI1A736fXF\nbLM+baXK+1Ll2MfVMJ/NDOUFAABAiTEZDjAOSEwBAABQYkyGA4wDElMAAAD8SvkmGmIyHGAcMPkR\nAAAAJJVzoiEmwwHGA5MfMfkRkx/1+FwmP+qOyY+Y/GiUmPxoeEx+NJ46vxf2NtFQqLoc1Xbr1Faq\nvC9Vjn1cMfkRAAAAcsBEQwDCIDEFAABAAxMNAQiDxBQAAAANTDQEIAwmPwIAAIAkJhrCXNGEWBOb\no2He+w9Js9toDygCiSkAAAB+pZF0TEcTzzQnPML4KeMszagvhvICAAAAyDCxuTUplaL7yzeFiQd1\nRmIKAAAAIAOzNGN0SEwBAKgYM1tnZg+a2UNmdk3G42vM7Odmdm/j77oQcQKoOmZpxuhwjikAABVi\nZvMlfVLSmyQ9Luk7ZvZld38gteo33f3t+W+fiVCA8TG7Tdq4onU4L7M0oxgkpgAAVMt5kva5+yOS\nZGa3S3qHpHRiasVsnolQgHHBLM0YJYbyAgBQLadJmk3cf6yxLMklXWBm95nZtJmdk9/mmQgFGCfu\ne6fj2Zndp9eTlKIo9JgCAFAt3sM6/yhpwt2fN7P1kv5W0iuyVjSzLYm7M+4+039ITIQCAOPIzNZI\nWpNHWSSmAABUy+OSJhL3JxT1mv6Ku+9P3N5uZp82sxPc/bl0Ye6+ZfiQmAgFAMZR48fMmfi+mV0/\naFkM5QUAoFrukXSWmZ1hZsdIerekLydXMLOTzMwat8+TZFlJ6WA27mu9z0QoAIDhkZgCAFAh7n5E\n0lWSdkr6vqQ73P0BM7vSzK5srPb7ku43s3+S9FeS3pNfBHuujiZAkaL/d2/mnDMAwLDMvZdTVXLY\nkJm7e88zBJrJ3WXp2/nGVEy5cdmSlEf5RcY5SPn9rN9t3fjxdusNU4/D1FvWc7vFOqhe23qvdZlX\nXIPIu87zKHdU2+jUZgYtc5gY0tsuQ/vIW7+fK5hrHD+b6/Y6GESv9TCK9+V+jWq7dWoro/gMLUqd\njsO4GOazmR5TAAAAAEBQJKYAAADoymxyymz9DrNLZ+L7gUMCUCPMygsAAICOoiR01dbW69iu2mo2\nKc4xRtlE7XVis7RskbT/kDS7jXZafiSmAAAA6GJic2tSKkX3pzZJ4gs/SiP7R5SNK/gRpfwYygsA\nAIAuli3KXr508WjjALpp9yPK8k1h4kGvSEwBAADQxf5D2csPHBxtHEA3/IhSVSSmAAAA6GJ2m7Rx\nX+uyDQ9Lj94YJh6gHX5EqSoSUwAAUHnJGWOj/8wYm6fo3Lw9V0tTO5pL797MOXsoH35EqSpz99Fs\naAwv4i1JVbigcb/l97N+t3Xjx9utN0w9DlNvWc/tFuugem3rvdZlXnENIu86z6PcUW2jU5sZtMxh\nYkhvuwztI2/DXMQbkbp8NreZ7GSftOfqZOJUhdfBKGYT7bUeivhsHtaojmEe2ynLzLCj+AwtyqDb\njOp++SZpel30Y8qjN/IjymgM89nMrLwAAKDi6jFjLLOJ9qcsiV8WjmVYjTqejhLb6fWh40FvSEwB\nAEDF1WWyk3ok2HlrTUDvVHOYdpkTP45l3ZT5h5C6IDEFAAAVV5fJTuqSYOcnu+dx1Vbpv/2i3Ikf\nx7JO6AEfja6TH5nZOjN70MweMrNr2qyzxszuNbO9ZjaTe5QAAABt1WWyk7ok2Hlq1/O47Izs9cuS\n+HEs64Vro45Cx8TUzOZL+qSkdZLOkXSZmb0qtc7xkj4l6W3uPinp9wuKFQAAYI65M8ZO7ajmjLF1\nSbDz1K7n8YU2k6uUJfGr7rFkhuss9ICPQrehvOdJ2ufuj0iSmd0u6R2SHkis815JX3T3xyTJ3Z8t\nIE4AAIC26jDZifveabNJRcNRmU000q7n8RcPSxuPb+3FKk/iV9Vj2WnIqrQ3WFzh0QM+Ct2G8p4m\naTZx/7HGsqSzJJ1gZt8ws3vM7F/lGSAAAMC4cN87HSfW7tPry57IFK9dz+Nz15e9l7yax5Ihq9mq\n2wNeJd16THu5yOlCSa+R9HuSlkjaY2Z3u/tDwwYHAACA8TW351FKJaCV7iVPCz/zK0NWs1S1B7xq\nuiWmj0uaSNyfUNRrmjQr6Vl3PyjpoJntlvQ7kuYkpma2JXF3xt1n+g0YADCezGyNpDWBwwAwYslh\n2on7tVOOmV8ZstpOHU4XKDtzb98pamYLJP1AUW/oE5L+XtJl7v5AYp1XKpogaa2kYyV9W9K73f37\nqbLc3ducqJ61bbm7LH07T0WVG5ctSXmUX2Scg5Tfz/rd1o0fb7feMPU4TL1lPbdbrIPqta33Wpd5\nxTWIvOs8j3JHtY1ObWbQMoeJIb3tMrSPvPX7uYK56vjZPIr3klEo+jtKL2UX8dk8jHbbLaKu8iyz\n9/pev0PavnbuI1M7hkmE+vvelpUcb3hYunuztPeuEMe8qPeDMpVZJ8N8NnfsMXX3I2Z2laSdkuZL\nusXdHzCzKxuP3+TuD5rZDknflXRU0s3ppBQAAABAJ+GH0XYasmqkYihYt6G8cvftkranlt2Uuv8f\nJP2HfEMDAAAAxkU5htEyZBWhdJuVFwAAAMAA4muA9nZNUGZ+LQLXZa2Orj2mAAAAAPrTPF9Tku68\nKPrffjIjZn4tSugJpdCrjpMf5bqhGk6w0K1sicmPmPyov20x+RGTHw2DyY/Qrzp+NjP5UX5lM/nR\ncGUOM5lRiMmXuj0vxOsnj/eD7EcGn1CqSu8jIRQ2+REAAACAQYSfzAjthD0G4a9XW04kpgAAADU2\n90vw9u5PQg7KMZkRsoQ7BuW4Xm05kZgCAADUVPaX4Gj5uH8J7mb4Xq3ZbdLGFXOvCcpkRv0a7lhs\n3BfiGLSPeWJz+vUY3Z/aJGmsX5MkpgAAoPaiL4l7ZXbpzHgNncv6EixJy8f+S3B3w/VqNScz0l3R\nEiYzGsTwPYx7rh71hFKdYpbOYYh3GySmAACg1vqdHbVeOM9xcMP3akXJaXyba4IOZrgexjDXZe0U\nM0O82+E6pgAAoObafUlcvilMPL3J59qLfAnOFwn96FXxx5VOMaevV3udpPccko785rhfY5UeUwAA\nUHPV+mIbDztuvdTIoD28Wec5SpznOCgS+tGr4o8r7WNuDvFe/TFp9+ukGyRpkaTXSBu3jsdIjmz0\nmAIAgEowm5warBexal9sJzbPXTZYD2/0BXfP1dG5dVL8f1y/+PYn2aslMXFRKOkeRqn8x6JzzNHr\n77ifzH1e+UdyFInEFAAAlF7zPNHta6PzRLevlVZt7S05rdoX23x7eN33Tsfn1nGeYz/SCf3dm0no\nRy/7x5VyH4veYq7WSI5RYCgvAACogMEnQJk7O+rqe6RF86Rz/tRs/ebyzdBbtR7eegozaQ6yVPFY\ndI+Z13kaiSkAAKiA4XoXkrOjSmcfX+6L289uk7S2dVmZe3jHz/DXOAV4nacxlBcAAFRAnr0L5Z6h\nt5ngVGfo4jgZblg5EOF1PheJKQAAqICizxMt33ldyfNCx/nLavlU8/JDKCde500M5QUAAKXXPE90\napM0vS7qXXj0xvy+yFXzvC6GlIbApDWj1NrG75TZ5BRtvJ5ITAEAQCXkNwHKxn2tPV7VPK+rOaS0\nzOfL1hGT1oxWuo2vGutrfdYZQ3kBAMCYqdalJ9pfu5UhpWFU7fJD+Rr8esKDoo2PC3pMAQDAWKnK\npSeiL/x7FU2uE0v2iDKkNITih5VHzNbvKNsQ7fL00tPG64jEFAAAoJQmNs9dlrx2K0NKQynyx43u\nP0iENPj1hPNFG89bGc5XJzEFAAC1UoYvWPno1iM6u03auKIO58siqdsPEiGF6KWvxznhZdZvT3hR\n77EkpgAAoGayv2BFvVBV0rlHdFRDSjFqZR6iHaKXfs/VzTYuFXFOeH1+zBpU7z3hRQ7nJjEFAAA1\n0/YLVsXMbpO0tnVZa29RVc6XRT/KPER79L30yTaeuJ+b8pw3G1I/P4YUN5ybWXkBAMAYKENvU6tu\ns5s2vxRXYwbh0c/WWlez2+YuK8fw1ajtVWtW6+6Y3bq/H0OK69GnxxQAAIyBMvQ2NfXTS+M+vb7s\nPaL0OuUnGqItRUlf+YZo16+XvsxDp0eln57w4nr06TEFAAA1U4VrTNatl6Zu+xNenPS5T68vS1Ja\nJ3EPv/T8b2evUa4fs4rUX094cdfxpccUAADUTHKylGZvU9QLVRZ166Wp2/6gzlp7+HdLulbSxxNr\nlPHHrEhREzX12hNe5KRrJKYAAKBWqjHUsMwT3AyibvuDekv28K9uLPuopH+rsg2dTirLkPmi3mMZ\nygsAADByxQ2HC6Nu+1NPTFAVS/fwr1aUlJZ96HS9h8zTYwoAADBidbsG6aj3p3U4450ym5yqat2N\nSll628q/5dH1AAAgAElEQVShXQ9/2dV7yDyJKQAA6Ev0BXevzC6dGc+L0eejGkOOezeq/clOsFZt\nHc8Eqx/FXX+yetrNQvu5FeFi6kW9h8wzlBcAAPSsmRRI0p0XSdvXNpKCMR0SiNGr93DGtPyG39a7\nt60f7WahDRtVL+o9ZJ4eUwAA0Ad6XRDauCVYeQ2/rXdvW7+yevjLNXP3XHU7BSCNxBQAAPRh3JKC\nVkVdqgH9GLcEK68fgtoNX61Hb9u4qNspAEkkpgAAoA/jlhQ0MXlMWZBgDfJDUN1721B9JKYAAKAP\ngycF1e9tZBhzGcxNsCTp7s3VakvDGuyHoDr3tqH6SEwBAEDPmkmB7oqW9NbrUo/exvEexlwmyQQr\ncb/0BvtxZuO+8e4dxrggMQUAAH2JktP4dq+9LnXobRzfYcwY3uA/zuy5muG3GAckpgAAYATq0NtY\n3LmNc3vStg9bJEpnsB9nGH6LcUFiCgAARqD6vY3FTh6T7kmLk9W9wxeNkqjDjzNAcUhMAQDACNRj\nJtXieq/SPWmStHxTfuUjvOr/OAMUicQUAAAUjktVDGK0PWnjNpx49LNE1+PHGaAoJKYAAGAkOFeu\nX8X2pLUmZk8tlc47Vbr1lPQ6dfzxIMQs0fw4A3RGYgoAABBc+pIgUqMnbV0RW5ubmF0n6YaMNZdX\naNbkfoSZJZofZ5qqf11j5G1e6AAAAACqwGxyymz9DrNLZ+L7+ZW+5+qoB02K/xf7JT2dmLXrq6jr\nxDxMRBRS84eR7WulOy+K/q/amu9rClVDjykAAEAX2UM/V23Na+hnuifNTD5smZ2lE7Mjbdar68Q8\nTEQUVh2ua4y8kZgCAICx131YYd2+SKcTs0skXSvp46n16joxT38TERU57DSr7LJdJii+dFE0WiCP\nGOmxrguz9Tta2+7gSEwBAKgYM1sn6a8kzZf0OXf/RMY62yStl/S8pPe7+72jjbI6epsIp25fpNOJ\n2WpJ/+eT0iVPSLteGw0nnl5X13P++p+IKLt9DJtAtmt7QxWas2aMUjTsVho+xvL2WKcTrbq+BobV\nvM7y9rXNpRtXSN8buEwSUwAAKsTM5kv6pKQ3SXpc0nfM7Mvu/kBinSlJK939LDN7g6TPSDq/t/IH\n6xnq9rzWx++U2eT10sSq5oywi0w6fn/jsZ5mgu3U09Sphyda9tRSabfMLvmuNHG2dPMxUam7Je2S\ndNpKaeI2s8lPRnH6a7KjeOq4dl9koxhO+Atp/tmSFkkyaYmkYw9KzzwlLf5Fcp+j8uKY4zpqud/3\nDLkZ9T4l7U0kZqs/Ju1+nfSWw9Kxy6IYpUGHEzfr/a3PSYcXSnZUOvKctHCptPhJ6eATIb7szz0W\nS0x6ReM4SNK8N0hnv8Hskh9LT3w03Y5aE8fdkk5aKZ39nxrrJNpysteoWe9zY4mP68S50s0nNsvd\nJWneymi9N/1Y+prM3vXdqN629/g6O+EvpF87UzpW0n//ifStRjt/cbl03Dxp0RHpizK75B+i9pd8\n/cW3tVDaJbO3/EJ6+SLp5oXRFq5TlD6ctLK5vex4mo+9dKo0/5TE8d8jnXCidMVB6ZbEDzvvfFI6\n8ptRrPH6X0rEeWdjH5Ovaal1f/f/WHrio/2/TtbviOLcpXSilTVcP2O/90gnvC1qX0uOlY59QTrw\nw7gttba/5OO7usTV7Xi3vPfsSbfDHuLOaD971XoMnn2x2SaSbej0+XMjvnml9Ll+qr51f90LPoUh\n3pCZu7v1vr7cPXpzTN7ON6Ziyo3LlqQ8yi8yzkHK72f9buvGj7dbb5h6HKbesp7bLdZB9drWe63L\nvOIaRN51nke5o9pGpzYzaJnDxJDedhnaR976/VypCzNbJel6d1/XuP9hSXL3f59Y57OSvuHudzTu\nPyjpInd/OlWWS2sO6FdJ0wuSzp4vXa7oC9Mzkn74krT4Remll6SdS6W37G98sX1ZVMrFB6QXFjef\n9wVJT0o6KMmPSvJmuZ9X9OV7taJ1L5e0VdLxkm5JPPa2o9ILL0iHjkrLjkbb299IGCRp713SBQel\nsxdHZcbe9YT0pVOlS+6XTl8h3bpE2iDpgKSfHZV2NCZ93CBpoaJ8/V83dv+vJH1a0ocacXxB0o8k\nvVzS6yXNSDpLrUNd//lh6RXzpVvnNZ/z46PS0cY+nzRfOiGx/smS1qb2OfamQ9LyedKtxzTr4QqX\nbrHm/akXpANHpPkm6ZgouYqT3G+tjI5FnHTFy886LroETLKM6WOli5+Pyjh+vvQ31kyKnpF0k6Tf\ne1r6ryc110sm1V9aGsV8yf3Rl9fDB6LE8/ABadnJ0t8siuryu5J+u/H/vW3ahhSVe/jF6HjH+5Z8\nv2q33cVPSl96dbRs/c5mAnT4gHTkZZItbJZ1/DzpBGs9Dl+QtF/SHYljfkDST13aadG+/9ai6PhK\nieOi1vZ6hUvvs+bznzwqnWnS51PvT3Fd/tZ86VaL2uFzkr6k1ra3U1E7i4dWX5c4Lh84Ir1/QfvX\n2UnzpXMaz0uWdama7Te5rfTrb2ejbm5Ts31ukfQvGvuqVDyXu3S5zY3nBY9e869PHf/4NZV8r7hL\n0gVHpbPnNdf/bCLOZN3Ht3dJ+n9flF62QDpnXvN4xjEsPiodfkHatTh6zzpizfevIyYdtOi9Zf8C\naXcjOf5QY5u7E2W9IOlYlw4fkuYdlpb8t6jNXf6CdPmxURw/UHSO9glqbVvJtvTOo9IJ85qPZ73W\n3vRIa7t9wZrvm/F+/8Uh6djDzbjT71VZ78FLPHr9PP4L6Zx/1ow7/f5+0KSVC6P3oLjtxa/j9zba\nRDL+hyX9deLYvEfRjwNf0KCfzczKCwBAtZwmaTZx/7HGsm7rnJ5d3JlLpTMXSGfOlyYayeVORecc\nHpF0wXzpw4ul0xtJwV3LpKtf1vr8+Hm3Kco1Tpb0W5LOnNcsN05Kb2s8L15/gZpfiuPH/td50sRi\n6cLjou198WXSrtdIK/9T9CdFjyeT0t2SXnVqdPu8V0vvXxLdPknS7ZJel/jOc7Kkyxq336voS9xu\nSd9sLIv34zcacX5TUW/NWkkfTWzztIXNpDR+zhmJfT5X0imJv/gLf7zPcdySdP6i6AuhGvUvRUlp\n8v4px0aJ5plLpAsWSHfNl/7NUukVK5vHIl7+paXSRSubSWmyDEm6YIn0xgXSuYltXCKp0XmnM05q\nlLkkah8XzG9uL/b6V0sfOVF6wxnS9hOl88+Qzm0Mef5mo27j/+3axgXzpWvmR8cz3re4PcaPZ213\n+4nNpFSSPrK2Gc/Lz5BW/Hqznt44P9rP5HHY2YjlrMbzb1PUVj4k6fRGnbxxifT+RFIat89kEidF\nSWny+SfPayal8fFN1vn7G4+5oh9Ikm1vl5o/fsQdd8nj8v4FnV9n5yaevytRxivUTFaT7TzZFuNt\n72rUU+yIWnv24vcGKUpKs+KJ3xPSxz9+Te2U9IeN50jSxLzm+p9NxRnHFovb6knHSOfOa76u4hhW\nSbqm8R4iSdcsa75/xbcvPC66fUpjnd2KEsm4Xj7XKOsNkq4x6fWLpR2/Jv2bRpu7/NhmHPPUfK3H\nsaTb0rnzWtte8pie33itJdvtmQua75txfLdJesOi1rjT9Zo+FhfMa75+zji1Ne7k+/vrlkqrj2u+\nByXbymfVbBPJ+I8m6kuK3mfj24MhMQUAoFp6HeqU/sW6zfOSidOr1PxyurqxLOuL6q7U81+VWOeW\nVJnx4/Hzkl+CT1HUW9XusfgL9m5FvTSLTpAmG12QcZnJmOL1FyRiTC6L3ZB4fLWkP278j8uM9yPe\nr1cl1v23iXLS+5Xe5wWNvxsaf6sb692RKCPuBUrGl56AKL6fLD+ZfJySeDz53AWJdZJlxGUuSKwT\nH/P0elnbS8aVrve4vDvUrNO4LWW1jWT7yvprt920dFKVfP4CtR6LeD9vSMR7i+a29WQ7Wp1Yntye\n1Npek+0lfl4y7mSZtyhqe59KPCfZDlYn/qe31e511q6dx9eoXa1m+0u+/pLrLFDrNW1vkPRoKq64\nLtrFk3ztJI9/1nuMNPe1tqtNbFJrPcf1eUsihvRxzLqdPla7FPU0KrFucr10O0zuQ/K1Hr/Gs9pS\n/Hj6tZb1+k6+b0pz34/bvVelj0Uy7k7v73HbTEoeg7hNJJ//x6n6Gh5DeRnKy1DeHp/LUN7uGMrL\nUN5RGuOhvOdL2pIYyvsRSUeTEyA1hvLOuPvtjftth/Jen7i/pvEHAEAvZhp/sY9JAw/lZfIjAACq\n5R5JZ5nZGZKekPRuNcelxr4s6SpJtzcS2Z+lk9LYx3Rt4vYlinpv7lDUQylFv5InbytxPxY/7yxl\nu0TR8C/LWD/9WHp716m1tyTefvy8rN6z+Jyn5HN3J/Ytvv3Hat3f9H5cklgvua1kfNe2if0SzR3W\ndiC1z/H5g3GMcVk3qDlUNL6frvPksUjXWfLcwqRkGfF6tzXiiusga1vttpc+LhsUDYF8t5p1G/9v\n1zayjnc/281anl4W7+ctis4Pjfc3uf9nZcSSbDPJ45o8Lsl14m1fotbzMZP7Eq9/VmJZsj2vTfxP\nbrtdPElxmztZzaGY8TGIl5uyX3/xeZHxxDVXJMqJ/UuX/tZa66LdsU2+drJeU+l9SL/W4jjjdpx+\nPUjN1/m7UzG0a1NZt9PH80NqPS88Xvbpxv2s4xvfj9tWui0l3wuyXmvtYmz3PtHpvSp9LNKvCSXu\nJ7eXfN+I22Bc76vV+hpKlpf+bJDmDtbpHUN5AQCoEHc/oijp3Cnp+5LucPcHzOxKM7uysc60pB+Z\n2T5Fs2t8qG2BejLxt1PSRYomBLqksUyJ29em7qef98PU8uTjayU9Jemh1Prpx65NbS/rN/RkmR9V\nNDlLUvyl8XvejHl1Y5sf8Obt/5La3/R+7Eysl9zWtw9Je/dH5SZjT9bnTknvU/Ql7SlJP5Y0ezS6\nlEIy9mcaMSXLuqIRY/J+us6Tx0Kp5cmhh7E/UmsZ1zbq5H2KzjW7XM1rmaa3ld7etYn/UvSFXGom\nNUcUnRuXrON2bSO5PbV5vN124+UbD7cuT5aV3E9JeuSoNPtStL/J/f+e5ra9uJ0kXw/J9vpHiXWU\neH587ml8fL/nc8v8oVr3Py73q5KelfSJl6THfipd/lJz35Lbavc6e5+kpxX9VrXVpR//RPr8C83l\n9ypKMJKvv9Maxyk5m+qTL0jf+mU0WdD/+Fy07KH/In3g+da6aHdsk6+drNdUsk6zXmtPKZpYJ36d\nfF3S3x+SZp6IZvNVop6PKHorzGpT3W6nj+d7UvtxhaIJkJLtMN7uv07dNzWOd2pbyfeCZzT3tdbu\nNZB+n4h1eq9KH4v0+0Qy7uT7e9y+4jZ474vRY/+x8f+rar4XX+5z406+JgfHUF6G8jKUt8fnMpS3\nO4byMpR3lMZ1KG+eGrPy7pdssRTXZTxD6rwXpSM/lb52RvTl9LCkQ/Mas+Qeji5DIUkLTo6ev9ha\nZ1a1Y5plxuUee1B6pvFt6Nf/mfTiwuwyk48dd6y0PZWdbpD0oyPRbJPHvCAd+EF07U0puh7l0sXR\n9RYnt0gnvDV1KYmvSMvPj9Z56rjGZRCOmTu7bDw75nELpRePRPWx5FHp0OPSozdG29p7l/TW5xKX\n5ThLujhRn/E+H/iB9MSfNy7RMiWd8LG5Me3dIl36zeg6jo/e3Ywxef/IadHlGmxhs36XKJrcZPUv\nm/X+sgXSF1NHe7ekv3xOevHxZhnHLWzWX1wvR06bO7NtclvxMVz8i6jedr1aWndQ2rF47vb+3bPS\nS0+2Ly/ZNg6/GM2gevHzndtOcrvvur9xPBL1s+vV0eymXzujWdZia87qO/mW5nFYvqm1juO2ktUu\n4+3+6nIZieXHH4ja24UPNdvurmWt8Z3w1uhSIcctiuo8fh0sODlaNreNtcYZH5d4JuILH2q+9rLq\nKt3m0vuabE+Lj5v7GoteS+7T66Xkd6BkPHFdxO8N6dd9+r0k6zXVbvmSR6VnX2hcxuZAI+5Eney9\nK5qdOW4Hz74g+a9nvx8tO9r6/pW+Pe9o6/bi42wLpWMWR7Pqxj2HWxTNJfc5RduP1/vGEuldB5qv\npXZtKfkelHytxe87yXabfg+J3+Pe+lxU7ryj2fWXfg/uNJt28v09rof4PW7vXdLkW6L/8XvT9Lpo\nWbJNxq/t+HXxpRMG/mx295H8RZvqZ333rNv5xlRMuXHZeZVfZJyDlN/P+t3WjR9vt94w9ThMvWU9\nt1useWyrU9m91mXIv7zrPO99K3IbndpMiPpPb7sM7SP//VXt9qmMdRi67UjnTkkbHoo/D6K/K/ZJ\n506l4xzFd4dOdVPUe2Dv9bLhoWj5uh2ty+O/9duLqAPpD2ayt/cHM0XVwdy6T+539H/YOg7VFvJu\nK93L6378Qr8PZNVB6/FufT8ovm7cO70G+99Ouv3OLWPQ78RR2Rfe04zxzf8Ql9/L96J+31eH+Wxm\nKC8AACit6OLve66O7l36zej/3ZvTF44fPxObo4vZJ928MurJmN0mbdzX+tiGh5s9vXnbfyh7+YGD\nxWyvVdSDtmqrtL1xbZTof7Qc3YU9fv1qHtfk8V61tZjj3a5upM6vwd5lt9/m/phNTpmt39G6fj9l\nv+Jm6cLXNpfueo103ufK+PogMQUAAKUWJ6Hud65J3h9vyxZlL1+6ePTJ/Oy2ucuKTITTshIEqd8E\nYXyN+oeMYU1snrus/4SwN+3qRur0GuxvG+0T3LlJq9RfEj6xWTrn1LnnnN96ShlfH8zKCwAAMITo\nS+LEZmm7op6N2W3FJ8+de7mic/GiZD46P7C4eOJtNc/vbZ4PWNQ2W+WVIIynxnmbkqY2Nc/PLvOo\nhNEd77l1E7VtSXfl19PcaX/aJa1TmyT1cHyWLWqf7pXv9UFiCgAAMKBmj0b85XH7WmnjCrPJgnt2\nZ7dJG1e0fmkN28sVT5QzetUailpGjbY6LcUTHZU1KZVGfbyTdROLfojJ6zXYaX+WtUkee00q9x+K\nZi1uV365dB3Ka2brzOxBM3vIzK7psN7rzeyImb0r3xABAADKKp/zzPrVHK47tSMarju1o9y9XEXK\nGm4plXcoKoYTeuh4JL/XYKeh1MMm4bPbpO8/MfdyLpc/WcbXR8ceUzObL+mTkt4k6XFJ3zGzL7v7\nAxnrfULSDg1zVVUAAIAhjH5YbbhhpFk9OeMoe7jl9LrxTNLrL/zQ8dZYNORrsN1w4ebywXtlG2Vs\nlJ79mPTW5CWq/rxZj+XRbSjveZL2ufsjkmRmt0t6h6QHUuttkvR/S3p93gECAAD0btTDahlGWgbp\nBCG+DjrqK9zQ8fy1S3DzOP+3Sj9gdRvKe5qk2cT9xxrLfsXMTlOUrH6msYg3AgAAEMioh9VWbUZT\nAFXivnfafXr9qGYlT16eJnmZmlHo1mPaS5L5V5I+7O5uZiaG8gIAgFIpblhtp2F4RW0T/QkzazJQ\nPdmTuUXLR/Ga6ZaYPi5pInF/QlGvadJrJd0e5aQ6UdJ6Mzvs7l9OF2ZmWxJ3Z9x9pt+AAQDjyczW\nSFoTOAyUTDLpaK/YYbVVGio3bsLNmgxUUcdrAgdPTO+RdJaZnSHpCUnvlnRZcgV3PzO+bWafl/SV\nrKS0se6WIWIFAIyxxo+ZM/F9M7s+WDAohblJhyRtPCzdvLB5n2G1423Y60AC9dN+FEHYawJ3TEzd\n/YiZXSVpp6T5km5x9wfM7MrG4zeNIEYAADCm5n6BSvaMZiUd/2qhtP5Zadn3GFaL0F+0gTzkORy9\n0ygCaSLoZG5dr2Pq7tvd/Wx3X+nu/66x7KaspNTdL3f3LxURKAAAGC/NL1DReU7J852i+1lJx2pJ\ny77nfueaaMIQktLxxqzJw0hPhNN87WFUst8HV20d/Fh0uvZy2GsCd01MAQAAwuh4vpNIOuqh2OSH\nWZMHlX9ChMF0SiQH0X4UQfRD3p6ro2vEXvrN6H/xMwHHup1jCgAAEEi3YZiz24a5+DzCK3pyImZN\n7t3c4aKn/gbn55ZB3sPRO/+gF/KawCSmAACgpLp/gRr24vMIrfjJiZg1ubvsHwje02bkAefnjlbe\nI0PK+4MeiSkAACip5Beo3ZJ2SbpB0pHfjK+rl0w6zOQkpVXD5ETlkPUDwco2x4Ch8qOVbyJZ5lEE\nJKYAAKCUml+gVn9MOutc6ZbGF+Vdr5E2buValHXAecLlkPUDwSWSrjjYfN1JZelZGydFJJJlHUVA\nYgoAAEor+lK2fnPrl2OJc93qorzDCsdL1g8EqyVd9z1p6tmy9ayNm7ImknkjMQUAYIwNe0280WC4\nZ12VeVjheGn3A8Fz17vv5lhgJEhMAQAYa/nOglqM+g73nDsTatl/JMjfuPQGlRk/EKAMSEwBABh7\nZR8WW8/hnkVfKgXoBz8QIDQSUwAAoDIPi61vb07xl0oBgKogMQUAACr7sNh69uZw7iyA8hvVKQfz\n8i4QAABUTfWHxVZTfc+dBerIbHIqSsyiieOihG0crNoanWogRf9XbS1i30lMAQAYa1M7pLs3V39Y\nbBXNbpM27mtdxo8EVTa+iUv9Nc8JLz5BK5+sUw6Wb8p7KwzlBQBgjLlPrw8dw7iq77mz44nJrOqO\nc8Jb5X/KAYkpAABAIPU8d3ZckbjUG+eEt8r/lAOG8gIAAGBgDF2NkbjU2zifEz6aUw7oMQUAAEDf\noiR0r1rPuRvnoavjnLiMg3peT7k3e64exSkHJKYAAAAYwMTmucvGeejqOCcu9TfO54SP6pQDElMA\nAAAMgKGrSeOcuIwLzgkvFokpAAAABsDQ1TQSF2BwTH4EAACAAcxum7uMoasABkOPKQAAAPoWDV2V\npKkdDF0FMCwSUwAAAAzMfXp96BgAVB9DeQEAAGrEbHIqeW3R0PEAQC9ITAEAAGoiurboqq2t1xaN\nlwNAeZGYAgAA1MbE5tbraF4XL7+N5BRAmZGYAgAA1EZ8bdHdjfs3NP5vP1FatZXkFEBZkZgCAADU\nRnxt0V0Zj928Ulq+aaThAECPSEwBAABqY3abtHFf+wsvLF080nAAoEdcLgYAAKAmomuLTkqauE3S\niXPXOHBw1DEBQC/oMQUAAKgR973T0uz7op7TpA0PS4/eGCYqAOiMHlMAAICaafacTm2Khu8eOCg9\nemOUtAJA+ZCYAgAA1FAjCS0kEY1m953YLG2X2fod0uw2kt5iUeeoO4byAgAAoGdRgrRqq7R9bbRk\n+1ouRVMs6hzjgMQUAAAAfZjYHF16JolL0RSLOkf9kZgCAIBKM5ucioY2Smbrd9CLVLRli7KXcyma\n4lDnqD8SUwAAUFkMcQxh/6Hs5VyKpjjUOeqPxBQAAFQYQxxHb3Ybl6IZNeoc9cesvAAAoMIY4jhq\nZb4UTV1nri1znQN5ITEFAAAVxhDHEIq8FM2gmsO64x707WuljSvMJlWHBK6MdQ7kiaG8AACgwhji\niBjDuoEqo8cUAABUFkMc0cSwbqDKSEwBAEClMcQREYZ1A1XGUF4AAADUAMO6gSqjxxQAAKAPdZ35\nteoY1g1UG4kpAABAj+o+82vVMawbqC6G8gIAAPSMmV8BoAgkpgAAAD1j5lcAKAKJKQAAQM+Y+RUA\nikBiCgAA0DNmfgVQf2aTU9HkbpLZ+h3R+fXFYvIjAAAwBzPPZmPmVwB1F2qSN3P3ospu3ZCZu7v1\nvr7cXZa+nW9MxZQbly1JeZRfZJyDlN/P+t3WjR9vt94w9ThMvWU9t1usg+q1rfdal3nFNYi86zyP\ncke1jU5tZtAyh4khve0ytI+89fu5grna1eHcLyVS1Eu452oSMJRB2d7TyhZP3VHfxYp+jNy+du4j\nUzvcp9d3fu7gn80M5QUAACnMPAsA4yvMJG8kpgAAIIWZZwGUT4jzHsdTmEneSEwBAEAKM88CKJfm\nKQbxENPta6VVW0lOixBmkjcmPwIAoELM7ARJd0h6uaRHJF3q7j/LWO8RSb+Q9JKkw+5+Xu9bmd0m\nbVzROpyXmWcBhNTuFIOpTZI49z1HoSZ5IzEFAKBaPizpq+7+l2Z2TeP+hzPWc0lr3P25fjfAzLMA\nyodTDEap8X4/0vd8ElMAAKrl7ZIuaty+TdKMshNTSYPPWhniSwkAtMcpBnXHOaYAAFTLSe7+dOP2\n05JOarOeS/qamd1jZhtHExoAFCXMeY8YHXpMAQAoGTP7qqSTMx66NnnH3d3M2l2Q/I3u/qSZ/Yak\nr5rZg+7+rbxjBYBR4BSD+iMxBQCgZNz9ze0eM7Onzexkd3/KzE6R9N/blPFk4/8zZvY3ks6TNCcx\nNbMtibsz7j4zTOwAUBROMSgfM1sjaU0uZbm3+6E1X2bm7t7zuS5mcvfo3Jjk7XxjKqbcuGxJyqP8\nIuMcpPx+1u+2bvx4u/WGqcdh6i3rud1iHVSvbb3XuswrrkHkXed5lDuqbXRqM4OWOUwM6W2XoX3k\nrd/Plbows7+U9BN3/4SZfVjS8e7+4dQ6SyTNd/f9ZnacpF2SPubuu1LrjWUdovrK9p5WtniAUIb5\nXOnpHFMzW2dmD5rZQ40ZANOP/09mdp+ZfdfM/s7MfnuQYAAAQFf/XtKbzeyHkv5F477M7FQzu6ux\nzsmSvmVm/yTp25L+n3RSCgBAmXTtMTWz+ZJ+IOlNkh6X9B1Jl7n7A4l1Vkn6vrv/3MzWSdri7uen\nyqHHdIiy6DGlx5QeU3pM+4mBHlP0gjpEVZXtPa1s8QChFN1jep6kfe7+iLsflnS7pHckV3D3Pe7+\n88bdb0s6fZBgAAAAAADjp5fE9DRJs4n7jzWWtXOFOCkZAAAAANCjXmbl7Xl2JDO7WNIHJL1x4IgA\nAAAAAGOll8T0cUkTifsTinpNWzQmPLpZ0jp3/2lWQUxJDwAYVJ5T0gMAgHLpZfKjBYomP/o9SU9I\n+nvNnfxouaSvS/qf3f3uNuUw+dEQZTH5EZMfMfkRkx/1EwOTH6EX1CGqqmzvaWWLBwhlmM+Vrj2m\n7u/WHikAAAxoSURBVH7EzK6StFPSfEm3uPsDZnZl4/GbJP25pJdJ+oyZSdJhdz9vkIAAAAAAAOOl\na49pbhuix3SosugxpceUHlN6TPuJgR5T9II6RFWV7T2tbPEAoRR9uRgAAAAAAApDYgoAAAAACIrE\nFAAAAAAQFIkpAAAAACAoElMAAAAAQFAkpgAAAMAAzCanzNbviG6v32E2ORU6JqCqSEwBAACAPkVJ\n6Kqt0va10ZLta6VVW0lOgcGQmAIAAAB9m9gs3byyddnNK6Xlm8LEA1QbiSkAAADQt2WLspcvXTza\nOIB6IDEFAAAA+rb/UPbyAwdHGwdQDySmAAAAQN9mt0kb97Uu2/Cw9OiNYeIBqm1B6AAAAACAqnHf\nO202KWlqUzR898BB6dEb3fdOh44NqCISUwAAAGAAjSSURBTIAUN5AQAAAABBkZgCAAAAAIIiMQUA\nAAAABEViCgAAAAAIisQUAAAAABAUiSkAAAAAICgSUwAAAABAUCSmAAAAAICgSEwBAAAAAEGRmAIA\nAAAAgiIxBQAAAAAERWIKAAAAAAiKxBQAAAAAEBSJKQAAAAAgKBJTAAAAAEBQJKYAAAAAgKBITAEA\nAAAAQZGYAgAAAACCIjEFAAAAAARFYgoAAAAACIrEFAAAAAAQFIkpAAAAACAoElMAAAAAQFAkpgAA\nAACAoEhMAQAAAABBkZgCAAAAAIIiMQUAAAAABEViCgAAAAAIisQUAAAAABAUiSkAAAAAICgSUwAA\nAABAUCSmAAAAAICgSEwBAAAAAEGRmAIAAAAAgiIxBQAAAAAERWIKAAAAAAiKxBQAAAAAEBSJKQAA\nAAAgKBJTAAAAAEBQJKYAAAAAgKBITAEAAAAAQZGYAgAAAACCIjEFAAAAAARFYgoAAAAACIrEFAAA\nAAAQFIkpAAAAACAoElMAAAAAQFAkpgAAAACAoEhMAQAAAABBdU1MzWydmT1oZg+Z2TVt1tnWePw+\nM/vd/MMEAABm9gdm9j0ze8nMXtNhva6f3UAVmU1Oma3fEd1ev8Nscip0TADy0TExNbP5kj4paZ2k\ncyRdZmavSq0zJWmlu58l6YOSPlNQrKVhZmtCx5CHuuyHxL6UVV32pS77IdVrX8bU/ZLeKWl3uxV6\n+eyuo7q07brsh5T/vkRJ6Kqt0va10ZLta6VVW0eRnHJcyqcu+yHVa1+GsaDL4+dJ2ufuj0iSmd0u\n6R2SHkis83ZJt0mSu3/bzI43s5Pc/el0YWYXH5ZknTe5RNLhF6VdMnvLL6QjJu2U2cW/lHRM9+f3\nYomkYw9KX5LZhQ9JC06WtKj3sl8+z+zio93KliSzi/f3V3a6rLgu8tz/2Mp5Zpccai1/iUXbPHRU\nmm+t2+wnnm7rJh9/04+lr6XW67cek8dkmHrLeu4Sa7aXyeulvUMe16xtNdu69PJdZhebeq73ottJ\nL/sw77B05Lm5x7HTayVdRlb8ee3bsNvo9ppPPn/JQunYF6RnnpK+lVNb6abT++aShdLhI834Xr7I\n7GIvLpZuMWa9twxS1rEH84utOtz9QUky61h1vXx219EaSTOBY8jDGtVjP6Tc92Vis3TzytZlN6+U\npjZJms5vO5nWiONSNmtUj/2Q6rUvA+uWmJ4maTZx/zFJb+hhndMlzUlMpTO7bO9kSWsl3bY4un/N\nskbOK+nMJV1i7dHJkj4uaffS6P4rVnZaO9vPJJ05v3vZknTm0rnr9SJdF3ntf7L8JyTNS5Sf3uag\n8XRbN358Z+Pxl5/Rut4g9Rgfk2HqLeu5cSxqbP+VH+0tnn63lWzrL1sondkhpm4xj0ryOK5dLN32\na3NjaPdaSZeRFX9e+5bHNjq95pPPv2BJo+0ulG5rvLcM21a66dSWLljSeGxhI5YljX0pNqSuMQ5b\nVvyaHNHvL9XTy2c3UEHLFmUvX5rDewuA0Lolpr3+qp7+dtDmead0KeYGSdcl1tuVuN3tub26IVH2\n6gHLXdbmeemyNWD5cVnJushr/5Plr1ZrnOltDhpPt3Xjxz+eevyUxONSf/UYH5Nh6i3ruTek1nnl\n/P7L7WVbybaebF+91mVR7aST5HFsF0O710q6jKzn5rVveWyj02s++fy4TRfx3tVOp7aUdWy6HZMi\ndHpvGaSsejOzryrKwNP+zN2/0kMRAXrEgVHYfyh7+YGxHEEB1I25t//8MrPzJW1x93WN+x+RdNTd\nP5FY57OSZtz99sb9ByVdlB7Ka2Z8UAIAcuXuY9ltambfkPQn7v6PGY91/exOrMtnMwAgV4N+Nnfr\nMb1H0llmdoaicZ/vlnRZap0vS7pK0u2ND8OfZZ1fOq5fHgAAKEi7z9VePrsl8dkMACiPjrPyuvsR\nRUnnTknfl3SHuz9gZlea2ZWNdaYl/cjM9km6SdKHCo4ZAICxZGbvNLNZSedLusvMtjeWn2pmd0nt\nP7tDxQwAQC86DuUFAAAAAKBoHXtM81D1i3yb2SNm9l0zu9fM/r6x7AQz+6qZ/dDMdpnZ8aHjzGJm\nt5rZ02Z2f2JZ29jN7CON4/SgmV0SJupsbfZli5k91jg295rZ+sRjpdwXM5sws2/8/+3dT4hVdRjG\n8e+jKWUGIYGVCs6ihYtgpFaV/YP+GGG1cxNS0CoqECxy07IIonZtclEWborEaBNSizYakeNfKgRX\nZdpCImmR4Nvi/O7McbjnOoWn83svzwcuc++5d2Z+zzx37uthrudIOinphKSXy/Z0vUzIkqoXSddL\nOixpTtIpSW+W7Rk76cqSqpM2ScvLmr8ot9P1UiPP5uFMy2yelrkMns01ZvFsrjPLSG+zOSJ6uwDL\ngdPARmAFMAds6vN79pDhDLBm0ba3gVfL9deAt4ZeZ8fatwCbgeNXWzvNSdjnSk8bS2/Lhs5wlSxv\nADvHPLbaLDRH2pwt11cDPwGbMvYyIUvGXlaVj9cBh4D7MnYyIUu6Tlpr3Al8Ahwot1P2UtMFz+ah\n1z4Vs7kjR8rXmgnzLGMvns2V5ZiQJV0nrTX2Mpv7/ovp/Em+I+ISMDrJdzaLDw6xjYUTBX4IPP3/\nLmdpIuJb4MKizV1rfwrYFxGXojkp+2ma/qrQkQXGH/yj2iwR8VtEzJXrF2lOeL+OhL1MyAL5evmr\nXF1J84/2CyTsBDqzQLJOACStB54APmBh/Sl7qYxn84CmZTZPy1wGz2bqzeLZXGGWPmdz3zum407y\nva7jsbUK4KCk7yW9ULatjYUjD58D1g6ztP+ka+230/QzkqWrlyQdlbSn9baBFFnUHDFzM3CY5L20\nshwqm1L1ImmZpDman/03EXGSpJ10ZIFknRTvAruAy61tKXupjGdzfabpeZ3xtWaeZ3M9WTyb51WV\nhR5nc987ptNwZKV7I2IzsBV4UdKW9p3R/J06Zc4lrL32XO8DM8AscBZ4Z8Jjq8oiaTXwGfBKRPzZ\nvi9bLyXLpzRZLpKwl4i4HBGzwHrgfkkPLbo/TSdjsjxIwk4kPQmcj4gjdJwWJVMvlZmGn4tnc53S\nvda0eTbPqyKLZ/OVX6L3RS5B37O57x3TX4ANrdsbuHKvuXoRcbZ8/B34nObPz+ck3Qog6Tbg/HAr\n/Ne61r64q/VlW7Ui4nwUNG8nGL01oOosklbQDL69EbG/bE7ZSyvLx6MsWXsBiIg/gC+Bu0jayUgr\ny91JO7kH2CbpDLAPeFjSXpL3UgnP5vpMxfM66WsN4Nnc+vSqsoBnM3Vl6XU2971jOn+Sb0kraU7y\nfaDn73nNSFol6aZy/UbgUeA4TYYd5WE7gP3jv0KVutZ+ANguaaWkGeAO4LsB1rdk5Yk/8gxNN1Bx\nFkkC9gCnIuK91l3peunKkq0XSbeM3j4j6QbgEeAIOTsZm2U0LIrqOwGIiN0RsSEiZoDtwNcR8SwJ\ne6mQZ3N9puJ5ne31f8Szub4sns11Zul9Nkf/R23aSnNEsNPA631/v2u89hmaI0nNASdG6wfWAAeB\nn4GvgJuHXmvH+vcBvwJ/0/x/oucmrR3YXXr6EXhs6PVfJcvzwEfAMeBo+QVYW3sWmqOwXS7PqSPl\n8njGXjqybM3WC3An8EPJcQzYVbZn7KQrS6pOxuR6gIUj/6XrpcYLns1Drn8qZvOYHCnnclmbZ3Nl\nWSbMs4ydeDYvMYvKJ5iZmZmZmZkNou+38pqZmZmZmZlN5B1TMzMzMzMzG5R3TM3MzMzMzGxQ3jE1\nMzMzMzOzQXnH1MzMzMzMzAblHVMzMzMzMzMblHdMzczMzMzMbFDeMTUzMzMzM7NB/QMF69JzBnKg\nowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f188ced3d50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (subfig1,subfig2) = subplots(1,2,figsize=(16,7)) # one figure with two horizontal subfigures\n",
    "subfig1.stem(xsharp)\n",
    "subfig1.set_ylim(0,1.1)\n",
    "subfig2.stem(x_restored)\n",
    "subfig1.set_title('$x^\\sharp$')\n",
    "subfig2.set_title('$x_\\mathrm{restored}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/gorthaur/lib/python2.7/site-packages/IPython/kernel/__main__.py:1: RuntimeWarning: divide by zero encountered in log10\n",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f188ced3c10>]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh4AAAFwCAYAAAD31XL9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm4XXV97/H3lwRQJgFBwhBkSlBECIMgdSCMBqsgtrRV\nK9Za67V1uK23DuhteXrrUFq9aFvaXgHrSKuiKFiRoBy0VAZJAmFORDTM8ySgDN/7x28dzknYOUnO\nXmuvPbxfz7Oevc/e+6z1zc/g/uQ3rchMJEmSemG9tguQJEmjw+AhSZJ6xuAhSZJ6xuAhSZJ6xuAh\nSZJ6xuAhSZJ6puvgERGzI+KCiLg6Iq6KiPfUUZgkSRo+0e0+HhExC5iVmUsiYhPgcuB1mXltHQVK\nkqTh0XWPR2benplLqucPA9cC23V7XkmSNHxqneMRETsB+wCX1HleSZI0HGoLHtUwy9eB91Y9H5Ik\nSSuZWcdJImJ94EzgS5l51irveTMYSZKGSGbGdH+3jsmlAXweuCcz/6zD+9lNgcMuIk7MzBPbrqNf\n2T6rZ9tMzfaZmu0zNdtn9br9Xq9jqOVlwO8Dh0TE4upYUMN5JUnSkOl6qCUz/ws3IpMkSWvBwNC+\nsbYL6HNjbRfQx8baLqDPjbVdQJ8ba7uAPjfWdgHDqus5Hmu8gHM8JEkaGv0wx0OSJGmtGDwkSVLP\nGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwk\nSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLP\nGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwkSVLPGDwk\nSVLPdB08IuL0iLgjIpbWUZAkSRpedfR4fA5YMNUHIjiuhutIkqQB13XwyMwfAfet4WOnRLBjt9eS\nJEmDrVdzPE4GTo0genQ9SZLUh3oVPP4W2BJ4e4+uJ0mS+tDM3lwmPgK7XQK/c3LEZQ9mnvfvvbmu\nJEnqRkTMB+bXdr7M7P4kETsBZ2fmizu8l5kZ5TkfAF4FHJ7JU11fWJIk9dTk7/XpqGM57RnAfwNz\nI2JFRLx1io9/EtgAeH+315UkSYOnlh6PKS+wSjKKYDbwE+D1mVzU6MUlSVKtWu/xWFeZrADeBpwR\nwda9vr4kSWpPK1umZ3IO8EXgGxFs2EYNkiSp93o+1DLxOusBXwUeBY7PpNlCJElS1wZuqGVctarl\neOAFwMfcXEySpOHX6t1pM3kEOAp4LfCRNmuRJEnN69EGYquXyd0RHA5cGMFM4ESHXSRJGk6t9niM\ny+R24GDK5mJfjWDLlkuSJEkN6IvgAU+Hj/nArcBVEfxWuxVJkqS6tbaqZerf4WXAacBVwLuqUCJJ\nklo2sKtaplLtaDoPuIHS+/GBCJ7VclmSJKlLfRk8ADJ5LJMTgJcCBwHXRPDbLruVJGlw9eVQS+fz\ncCjwKcqGYx/I5IddFydJktbJUA61dJLJD4B9gX8CPh/BdyLYq+WyJEnSOhiY4AFlt9NMvkTZ7fQ8\nYGEEX4hgp3YrkyRJa2Oggse4TH6VyaeBOcBNwOURnOzdbiVJ6m8DGTzGZfJgJn8J7EHZhfXaCP46\nglktlyZJkjoY6OAxLpM7MnkXZQXMNpQA8oUI9mu5NEmSNMlQBI9xmSzP5B3ArpTNx74RwUUR/E4E\n67dcniRJI29gltNO79rMBI4B3gvsDJwCfDaTu9uoR5KkQTcyy2mnI5MnMjkzk1cCRwNzgWURnBrB\ny9yMTJKk3hrqHo9OInge8IfA8cCGwJeAL2WyrNXCJEkaAN1+r49c8BhX9XbsC7wZeAPwM+CLwH84\nFCNJUmcGjxpUc0GOBH4feDVwISWEnJPJY23WJklSPzF41CyCzYDXU3pC5gFfAf41k6taLUySpD5g\n8GhQBDsCfwS8jTIUcwrwtUweb7UwSZJaYvDogWoo5jXAe4DdgE8Dp2byQKuFSZLUYy6n7YFqWe5Z\nmRwKHAvsB9wYwd9HMLvl8iRJGhgGj3WUyeWZvBHYBwjgimp79j1aLk2SpL5n8JimTH6RyfuAXYDr\ngAsi+GYEL2m5NEmS+pbBo0uZ3J/Jxyhbsv8AODOChREc6s6okiStzMmlNYtgA+BNwAeB+4CPUfYD\nearVwiRJqoGrWvpUBDMoE1FPADYAPk7ZFfWJVguTJKkLBo8+Vw23HEkJILOBk4B/c0dUSdIgan05\nbUQsiIjrImJZRHyg2/MNm0wyk+9lcjBlN9TXALdFcEYEx0ewZwTrt1ymJEk90VWPR0TMAK4HDgdu\nAS4D3pCZ1076zEj3eHQSwSzgtZR2m0fpCVlO2R31JuDnwD3AQ8CDkx7Hnz/snBFJUhtaHWqJiIOA\nv8rMBdXPHwTIzE/UVeAoiGATYC7wfGCn6nELYDNg0+px8vONgEeYCCJ3UYLLcmDZ+PNMHuzln0OS\nNPy6/V6f2eX1twdWTPr5ZuDALs85cjJ5GFhUHWsUwXrAJkyEkW0oW7nvBvzu+PMIHmYiiFwNXFkd\nt2fS7OQeSZI66DZ4+OXVgmqYZXzoBeBaYGzyZ6pJrdsCc6rjRcACYC8gIp4OIePHNZk80ov6JUmj\nq9vgcQusdK+S2ZRej5VExImTfhzLzLEur6s1qHo0bq2OC8dfrwLJNpQAshcwH3g3sHsENwI/Gj8y\nV+rNkiSNoIiYT/muqOd8Xc7xmEmZXHoY5QvuUpxcOpCqlTV7Aa8AXlk9PkwJIT+sHm9wiEaSRlvr\n+3hExFHAycAM4LTM/HidBaodVc/I7kyEkFcAz2ZSjwhwRSZPtlakJKnnWg8ea7yAwWNoRLAjK/eI\nbAf8mBJCLgQuzeTx9iqUJDXN4KHWRLA18HJKCDmEcqO8HwLnAwuB6xyakaThYvBQ36iCyKHAEdUx\nE/j++JH5zInHkqTBYvBQX6rmiOxGmXh8GKVH5F5KCPkBcEEmd7dXoSRpOgweGgjVpmcvpoSQQynD\nM8spQzILgYu8cZ4k9T+DhwZStXz3pUwMy+wJ/DclhJwPXOn9aCSp/xg8NBQi2JwyHDMeRDajDMss\nBBY6P0SS+oPBQ0Mpgp2YCCGHAXcyMSwzlslD7VUnSaPL4KGhF8EMYB8mgsgBwGImgshlmTzRXoWS\nNDoMHho5EWxE2cRsPIjMptwkbyFl2/5lmTzQWoGSNMQMHhp5EcwCDq+OvSl3430U+MUqx8+BZZR7\nzjzaTrWSNNgMHtIqqj1EtgaeD+xYHc+vjrnALpSbGl676pHJ/W3ULEmDwuAhraNqKe8uwAs7HA8A\nVwPXVI9XA0udzCpJhcFDqkm1ydls4EXVsQdlf5E9gFuARZRJrYuAxe68KmkUGTykhkUwE3gBZWXN\nPsC+1eOtlJvi/ZCyBfytrRUpST1i8JBaUC3x3YuyuuZgYD5wM3BudVyUya9aK1CSGmLwkPpA1Svy\nEuBVwALKfJELKSHke5n8tMXyJKk2Bg+pD0XwXMoeIwuq4yEmekPGMvlli+VJ0rQZPKQ+V01a3YuJ\nELIfcDFVbwhwdSbN/ocoSTUxeEgDJoLNKDfEWwAcBWzEyitmFgE3GkYk9SODhzTAqs3OtqOslBlf\nLbMv5e68S4DLKb0jP/YOvZL6gcFDGkIRbEUJIfsDB1XHY5R70fykOi7P5N7WipQ0kgwe0gioekZ2\nowSR/SnzRPYF7qb0ioyHkUWZ3NdWnZKGn8FDGlHVpNW5lBAyHkb2Ae6g6hFhIox4t15JtTB4SHpa\ntbHZXFbuGZlH2WX1J5QJrDdUx41uciZpXRk8JE2pCiMvoASReZRgMody197bgGWTjhuqx5syebyV\ngiX1NYOHpGmp7tL7fEoImcNEIJlDWWmzghJCbqQElDuq4/bx5/aYSKPH4CGpdhFsCOxMCSM7A9tU\nx6xJj88DfsmkIDLp+fjjrcCyTB7u8R9BUkMMHpJaUU1u3YLOoWT8cQdgV+Be4PoOxy8yebLnxUua\nNoOHpL5WBZQdgd07HFsAVwFXUDZMWwIstYdE6l8GD0kDq9o+fi/KpNe9q8cXATczEUbGH29xG3mp\nfQYPSUMlgpmUuSWTw8g8YAbwY+DC6licyRNt1SmNKoOHpJEQwXbAy4GDq2MH4PvAd4HvZbKixfKk\nkWHwkDSSIpgFHEm5y++RlJU03wXOBf7Lpb5SM1oNHhFxHHAiZXOil2TmoroLlKQ1qTZJ258SQhYA\ne1CGY84Fzs3kxhbLk4ZK28HjBcBTwL8C7zN4SOoHETwXOIKJIPIAVQgBxjJ5tMXypIHWF0MtEXEB\nBg9Jfahazrs3JYAcRbmR3tVM3EhvEXB9Jo+1VqQ0QAwekrQOItiUEj72Y+ImertStoi/ujquqY67\ngIeBh93oTCoaDx4RsZCyA+GqTsjMs6vPGDwkDazqvjVzKHuIjB8vBLYENgU2Bh6lhJCHJj2u7vn4\n432UPUlWAPe4D4mGQbff6zPX9IHMPGK6Jx8XESdO+nEsM8e6Pack1aW6E+94L8fXVn2/Gq7ZiBJC\nNqkeV/f8ucBOk55vD8wGnhXBzUwEkV8A1wHXAtdl8svG/oBSFyJiPjC/tvPVONTyvzLz8g7v2eMh\naeRFsAkTIWQHSjh5AWUFzhzgTkrwuWLScYObpKnftL2q5VjgM8BWlFnjizPzqDoLlKRhVy0H3pky\nxLMXZTLs3sB2lB6RKyj3tLkWuAFYUfXSSD3XF5NLp7yAwUOSpqXqJXkxJYTsUR27AdsCtwA3TjqW\nU+5p81PnkqhJBg9JGjERbEC54+8uk465lNU6z6EsEV4EXEbZt+SOlkrVEDJ4SJKeFsHWwL7VcSBl\nUuDPgYXAfwI/dN6IumHwkCStVnW33wMoO7m+hjKX5NvAl4ELMnmqxfI0gAwekqS1FsHzgWOBt1KW\nCH8W+LdM7my1MA0Mg4ckaZ1FEJShmHdQgsh3gE9ksrTVwtT3uv1eX6/OYiRJgyGTzOTiTN5K2VPk\nCuC8CL4Vwf7tVqdhZvCQpBGXyf2ZnERZHXM+8K0Ivh7BC1suTUPI4CFJAiCTRzP5B8pOqpcCF0Zw\negQ7tlyahojBQ5K0kkweqXpA5gK3AosjODmC57VcmoaAwUOS1FE1BPMRylbu6wHXRPDeaot3aVpc\n1SJJWisR7A78K7Ax8MeZLG65JLXAVS2SpJ7I5HrgEOCfgXMj+LsINm65LA0Yg4ckaa1Vy3BPp9y8\nblvgqggWtFyWBohDLZKkaYvgVZQekIuBP/OGdMPPoRZJUmsy+R6wJ7ACWBrB26pdUaWO7PGQJNUi\ngnnA/wMeAX4/k5tbLkkNsMdDktQXMlkCHAScB1wcwQEtl6Q+ZI+HJKl2ERwNnEpZdntW2/WoPt1+\nr8+ssxhJkgAy+Xa12uU7EWycyZfbrkn9weAhSWpEJosiOJxy19sZmXyh7ZrUPoOHJKkxmVwdwWHA\nWAT3ZnJO2zWpXU4ulSQ1KpPrgGOAz0XwG23Xo3YZPCRJjcvkEuB44MwIdmq5HLXI4CFJ6olMvgt8\nHDg7gk3brkftcDmtJKlnql1N/wXYHnhdJk+0XJLWkRuISZIGRiYJvAvYAPiM26uPHns8JEk9F8Fm\nwEXAP2XyL23XM+wi2AR4PrARJfT9ErgXuCWTJ9ftXN19rxs8JEmtiGAuJXwcmcnitusZBhFsCOxH\n2bp+b2A3YFdgU+DnwMPA48AmwJbAFsC1wH9Ttrofy+Thqa9h8JAkDagIfhf4KLBfJg+0Xc8gqoar\nDqesGnoNcCPwY2ARsLw6bquGuVb93U2BFwOvBI4EXgL8BDgH+GImdz7zdwwekqQBFsEpwPOA4zp9\nOaqzCNaj7I/yYeDZlEm7Z2Zyaxfn3ASYD/wWcCywkHLPnYWZPFU+Y/CQJA2wCJ5FGXL5fCafabue\nQRDBkcCngF9ReozOGg8GNV7jOcAbgXdUL30E+A7EUwYPSdJAi2BX4GLg8EyuaLueflUNjfwT8HLg\nvcA5TfcSVUM5xwB/AzwIcZDLaSVJAy2TnwLvA74SwbPbrqcfRTAPuBx4DNgzk7N7MTSVSWZyFmWy\n6j90e76uejwi4u8oE1l+DfwUeGtmPrDKZ+zxkCStUfUv668A92Tyrrbr6ScRHAx8HXhvJl9pt5Z2\nNxA7D3hRZu4N3AB8qMvzSZJGVPWv93cCr4ngtW3X0y8iOIgSOn6v7dBRh66CR2YuzMzxySyXADt0\nX5IkaVRlcj/wJuC0CF7Udj1tq+a+fAN4Sybfb7ueOtQ5x+MPgf+s8XySpBGUyUXAnwPnRLBN2/W0\nJYItKd+rf505PN+vM9f0gYhYCMzq8NYJmXl29ZkPA7/OzI5dQBFx4qQfxzJzbN1LlSSNiky+FMEc\n4PwI3pLJorZr6qVqvsupwLmZ/HO7tcR8yt4e9Zyv2+W0EfEHwNuBwzLzsQ7vO7lUkrTOqi/f44GT\ngNOBv8zk8Xar6o0I3gy8H9g/k1+1Xc9krW4gFhELgE8CB2fm3av5jMFDkjRt1XDL5yg3ODsuk7ta\nLqlREewPfJc+vYdN26ta/oFyo5mFEbE4Ik7p8nySJK0kkzuA11J2N70sgn1bLqkx1Y3zzgb+qB9D\nRx3cuVSSNDAiOA44BfifmXy57XrqFMG2lHD1sUxObbue1en2e32Nk0slSeoXmXwtguuAb0WwfSYn\ntV1THar7opwLnNbPoaMOBg9J0kDJZGkErwDGIvh1Jie3XVM3qpvkfQv4IfCxlstpnEMtkqSBFMFs\nyo3l3pnJt9uuZzoimAF8FXgCeGMmT7Zc0hq1uqplrS5g8JAkNSSCA4BzgIOqG80NjGq58CnAXODV\n/bZsdnXaXtUiSVJrMrkU+Fvgs9UX+UCoav0EcCBw7KCEjjoYPCRJg+5kYGPKZpZ9rwodJwFHAkdk\n8mDLJfWUQy2SpIEXwZ7ABcC8TG5pu57VqULH31O2ID8ik3vbrWjdOdQiSRp5mVwFfBb4aNu1rMH/\nBg4BDh/E0FEHezwkSUMhgs2A64Hf7MebykXwJuBvKBNhb2+7nulyVYskSZUI3gH8ViZHtl3LZBFs\nB1wBHJbJlW3X0w2DhyRJlQjWB5YBb8jkx23XMy6CfwEezOT9bdfSLYOHJEmTRPA/gKMzeXXbtcDT\nvR1XAXMyuafterrl5FJJklb2OWDvCF7cdiGVtwP/MQyhow4GD0nSUKk24zod+KO2a6mGfv6YskOp\nMHhIkobTacCbqhuwtekY4KeZLG25jr5h8JAkDZ1MbgIuB17fcilvpoQgVQwekqRhdSotbqMewRaU\nHUq/2VYN/cjgIUkaVt8C9ohgTkvXPxY4f9TuxbImBg9J0lDK5NfAF2hvkukbgDNaunbfch8PSdLQ\nimB34EJg+0ye7OF1ZwNLgB0yebRX1+0F9/GQJGk1MrkeuBM4oMeXfjtwxrCFjjrMbLsASZIa9p/A\nq6E3W6hHsAvwJ/Q+7AwEezwkScNuPHj0yl8Bn8nkxh5ec2A4x0OSNNSq3UPvBPbI5LaGr7UbcDGw\nayYPNHmttjjHQ5KkKWTyOLAQWNCDy30I+MdhDR11MHhIkkbBQuCwJi8Qwc7A64BPN3mdQWfwkCSN\ngh8Ah0bQ5ND/nwKnZnJfg9cYeAYPSdIouBF4DJjXxMmrQHMc8MUmzj9MDB6SpKGXSVJ2MX1rQ5c4\nAHgEuLqh8w8Ng4ckaVR8DnhjBM9q4Ny/A3ytCjiagsFDkjQSMvk5sIhy87baRLAeZZjla3Wed1gZ\nPCRJo+Q04G01n/NQ4J5MltZ83qE07eAREf8nIq6IiCUR8f2ImF1nYZIkNeAsYO9q6Wtd/pQSaLQW\npr1zaURsmpkPVc/fDeydmc+49bA7l0qS+kkEpwGLM/nHGs51AHAmMHdUbgjX2s6l46Gjsglw93TP\nJUlSD10AHFLTud4H/N2ohI46dHWvloj4KPBmyhKil2bm/R0+Y4+HJKlvRLADsATYuptVKBFsDSwD\ndszkwbrq63fdfq9PGTwiYiEwq8NbJ2Tm2ZM+90Fg98x8xvpog4ckqd9EcCtwULXSZbrn+APgNzM5\nrrbCBkC33+szp3ozM49Yy/N8hXLb4Y4i4sRJP45l5thanleSpCZcCewF0w8ewNHAN+spp39FxHxg\nfm3n62Jy6ZzMXFY9fzdwQGa+ucPn7PGQJPWVCE4CHsjko9P8/fUocxv3yOT2Wovrc432eKzBxyNi\nd+BJ4KfAO7s4lyRJvXQl8Noufv+FwL2jFjrqMO3gkZm/XWchkiT10HWUFSnT9XLgv2qqZaS4c6kk\naRQtB3ar7io7HQcBP66xnpFh8JAkjZxM7gceA7aZ5ikOBC6pr6LRYfCQJI2q5cBu6/pLEWwOzAau\nqr2iEWDwkCSNqmkFD2B/YFEmT9Rcz0gweEiSRtUyYM40fs9hli4YPCRJo2q6PR4HYPCYNoOHJGlU\nrXPwqFbBHAhc2khFI8DgIUkaVctY9yW1uwG/BlY0U9LwM3hIkkZSJvcBTwBbr8OvvRL4UTd3tR11\nBg9J0ihb1+GWQ4AfNVTLSDB4SJJG2VoHjwg2A36TEbgjbZMMHpKkUbYuS2pfD4xlckeD9Qw9g4ck\naZSty1DL0cA3GqxlJBg8JEmjbK2CRwQbAIcB5zZe0ZAzeEiSRtlyYM5aLKndHbg1k7t6UNNQM3hI\nkkbZPUACz13D5+YC1zdfzvAzeEiSRla1H8faDLfMAW5ovqLhZ/CQJI26tQkeczF41MLgIUkadctZ\n85LauZSlt+qSwUOSNOrOBd4ewbZTfGYWcEuP6hlqBg9J0kjL5CLgK8BfTvGxrSgTUdWlyGz2PjcR\nkZm5Lnf+kySppyLYCvgZsE0mj6zy3vrAI8CGmTzVRn39pNvvdXs8JEkjL5O7gV/Qea7HlsB9ho56\nGDwkSSqup2wUtqqtgLt7XMvQMnhIklRch8GjcQYPSZKK64EXdHj9uTixtDYGD0mSCodaesDgIUlS\ncT0wt8MN41xKWyODhyRJQCb3AY/BMzYS2wK4r/cVDSeDhyRJEzpNMN0MeKCFWoaSwUOSpAmdJphu\nBjzYQi1DyeAhSdKE24DnrfKawaNGBg9JkibcDzxnldcMHjXqOnhExPsi4qmI2LKOgiRJatH9wOar\nvGbwqFFXwSMiZgNHAD+vpxxJklrVKXhsCjzUQi1Dqdsej08B76+jEEmS+sAD2OPRqGkHj4g4Brg5\nM6+ssR5Jktq0Uo9HtZmYwaNGM6d6MyIWArM6vPVh4EPAkZM/XmNdkiS1YdWhlg2BzORXLdUzdKYM\nHpl5RKfXI2JPYGfgiogA2AG4PCIOyMw7O3z+xEk/jmXm2HQLliSpQasGj5Hv7YiI+cD82s6Xmd2f\nJOJnwH6ZeW+H9zIz7Q2RJPW9CGZStk1fP5OMYDfge5ns2nJpfaPb7/W69vHoPr1IktSyTJ6gBI9N\nqpc2ZcR7POo25VDL2srMXeo4jyRJfWB8uOUhylCLS2lr5M6lkiStbAUTPR4jP8ejbrX0eEiSNCwy\nOWjSjwaPmtnjIUnS6hk8ambwkCRp9QweNTN4SJK0eq5qqZnBQ5Kk1XNVS80MHpIkrZ5DLTUzeEiS\ntHoGj5oZPCRJWj2DR80MHpIkrZ7Bo2YGD0mSVs9VLTUzeEiStHr2eNTM4CFJ0uptBjzcdhHDxOAh\nSdLqbQD8qu0ihonBQ5Kk1ZsBPNl2EcPE4CFJUgcRBAaP2hk8JEnqbD0gM3mq7UKGicFDkqTOZgJP\ntF3EsDF4SJLUmcMsDTB4SJLUmT0eDTB4SJLUmT0eDTB4SJLUmT0eDTB4SJLUmT0eDTB4SJLU2Qzs\n8aidwUOSpM5mYo9H7QwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9H\nAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HA6YdPCLixIi4OSIWV8eCOguTJKll9ng0YGYXv5vA\npzLzU3UVI0lSH7HHowHdDrVELVVIktR/7PFoQLfB490RcUVEnBYRm9dSkSRJ/cEejwZMGTwiYmFE\nLO1wHA38M7AzMA+4DfhkD+qVJKlX7PFowJRzPDLziLU5SUScCpw9xfsnTvpxLDPH1ua8kiS1yB4P\nICLmA/PrOt+0J5dGxLaZeVv147HA0tV9NjNPnO51JElqiT0eQNVZMDb+c0T8VTfn62ZVy99GxDzK\n6pafAe/ophBJkvqMPR4NmHbwyMzj6yxEkqQ+Y49HA9y5VJKkzuzxaIDBQ5KkzuzxaIDBQ5Kkzuzx\naIDBQ5KkzuzxaIDBQ5KkzuzxaIDBQ5KkzuzxaIDBQ5KkzuzxaIDBQ5KkzmZg8KidwUOSpM5m4lBL\n7QwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwe\nkiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1Zo9HAwwekiR1NgN7PGpn8JAk\nqbOZ2ONRO4OHJEmd2ePRAIOHJEmdbQg83nYRw8bgIUlSZ88F7m67iGFj8JAkqbPnAXe1XcSwMXhI\nktTZ1hg8amfwkCRpFRFsAGwC3Nd2LcPG4CFJ0jNtBdyTyVNtFzJsDB6SJD2TwywNMXhIkvRMWwN3\ntl3EMDJ4SJL0TPZ4NMTgIUnSM90MnNN2EcMoMnP6vxzxbuBPKHvZfyczP9DhM5mZMf0SJUlSv+j2\ne33aPR4RcQhwNLBXZu4J/P10zzXKImJ+2zX0M9tn9Wybqdk+U7N9pmb7NKeboZZ3Ah/PzMcBMtOx\nsOmZ33YBfW5+2wX0sfltF9Dn5rddQJ+b33YBfW5+2wUMq26CxxzglRFxcUSMRcT+dRUlSZKG08yp\n3oyIhcCsDm99uPrdLTLzpRHxEuCrwC71lyhJkobFtCeXRsR3gU9k5oXVz8uBAzPznlU+N/3Zq5Ik\nqe90M7l0yh6PNTgLOBS4MCLmAhusGjq6LU6SJA2XboLH6cDpEbEU+DVwfD0lSZKkYdXVPh6SJEnr\notGdSyP8QJlbAAAEJElEQVRiQURcFxHLIuIZm4sNu4g4PSLuqHqFxl/bMiIWRsQNEXFeRGw+6b0P\nVW11XUQc2U7VvRMRsyPigoi4OiKuioj3VK/bRkBEPCsiLomIJRFxTUR8vHrd9qlExIyIWBwRZ1c/\n2zaViLgpIq6s2ufS6jXbpxIRm0fE1yPi2uq/rwNtnyIidq/+3owfD0TEe2prn8xs5ABmAMuBnYD1\ngSXAC5u6Xj8ewCuAfYClk147CXh/9fwDlAm6AHtUbbR+1WbLgfXa/jM03D6zgHnV802A64EX2kYr\ntdFG1eNM4GLg5bbPSu3z58CXgW9XP9s2E23zM2DLVV6zfSba4vPAH1bPZwLPsX06ttN6wG3A7Lra\np8kejwOA5Zl5U5ZNxv4dOKbB6/WdzPwRcN8qLx9N+QtP9fi66vkxwBmZ+Xhm3kT5H+6AXtTZlsy8\nPTOXVM8fBq4Ftsc2elpmPlI93YAS5u/D9gEgInYAXg2cCoxPYrdtVrbq5H7bB4iI5wCvyMzTATLz\nicx8ANunk8Mp3+UrqKl9mgwe2wMrJv18c/XaqNsmM++ont8BbFM9347SRuNGqr0iYidK79Al2EZP\ni4j1ImIJpR0uyMyrsX3G/V/gL4CnJr1m20xI4PyI+ElEvL16zfYpdgbuiojPRcSiiPhsRGyM7dPJ\n7wFnVM9raZ8mg4ezVtcgSx/VVO00Em0YEZsAZwLvzcyHJr836m2UmU9l5jxgB8pOwYes8v5Itk9E\nvAa4MzMX88x/1QOj2zaTvCwz9wGOAv40Il4x+c0Rb5+ZwL7AKZm5L/BL4IOTPzDi7QNARGwAvBb4\n2qrvddM+TQaPWyhjQuNms3IiGlV3RMQsgIjYFrizen3V9tqhem2oRcT6lNDxxcw8q3rZNlpF1Q38\nHWA/bB+A3wCOjoifUf41dmhEfBHb5mmZeVv1eBfwTUrXt+1T3AzcnJmXVT9/nRJEbrd9VnIUcHlO\n3Iutlr8/TQaPnwBzImKnKjX9LvDtBq83KL4NvKV6/hbKRmzjr/9eRGwQETtT7oVzaQv19UxEBHAa\ncE1mnjzpLdsIiIitxmeNR8SzgSOAxdg+ZOYJmTk7M3emdAX/IDPfjG0DQERsFBGbVs83Bo4ElmL7\nAGV+GbAiyuaXUOYxXA2cje0z2RuYGGaBuv7+NDwb9ijKSoXlwIfanJnbxlH9D3YrZYO1FcBbgS2B\n84EbgPOAzSd9/oSqra4DXtV2/T1on5dTxueXUL5QFwMLbKOn/6wvBhZV7XMl8BfV67bPyu10MBOr\nWmyb8mfdufp7swS4avz/f22fldpob+Ay4ArgG5RVLbbPxJ93Y+BuYNNJr9XSPm4gJkmSeqbRDcQk\nSZImM3hIkqSeMXhIkqSeMXhIkqSeMXhIkqSeMXhIkqSeMXhIkqSeMXhIkqSe+f+FnT36iGaBpQAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f188c044c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(log10(En_array-En_array.min()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convergence is much slower in this setting."
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}