{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "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": { "deletable": true, "editable": true }, "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": { "collapsed": false, "deletable": true, "editable": true }, "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": { "deletable": true, "editable": true }, "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": { "deletable": true, "editable": true }, "source": [ "We initialize the random number generator for reproducibility." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "random.seed(0)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Dimension of the problem." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "N = 400" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Number of measurements." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "P = round(N/4)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "We create a random Gaussian measurement matrix $A$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "A = randn(P,N) / sqrt(P)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Sparsity of the signal." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "S = 17" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "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": { "collapsed": false, "deletable": true, "editable": true }, "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": { "deletable": true, "editable": true }, "source": [ "We perform random measurements $y=Ax^\\sharp$ without noise." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "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": { "deletable": true, "editable": true }, "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": 9, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "def prox_gamma_g (x, gamma) :\n", " return x - x/maximum(abs(x)/gamma,1) # soft-thresholding" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Display the 1-D curve of the thresholding." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(-1.0, 1.0, -0.80000000000000004, 0.80000000000000004)" ] }, "execution_count": 10, "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": [ "" ] }, "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": { "deletable": true, "editable": true }, "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": 11, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "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": { "deletable": true, "editable": true }, "source": [ "We set the values of $\\gamma$ and $\\rho$.\n", "Try different values to speed up the convergence." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "gamma = 0.1 # try 1, 10, 0.1\n", "rho = 1 # try 1, 1.5, 1.9" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Number of iterations." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "nbiter = 700" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "__Exercise__ \n", "\n", "Implement nbiter iterations of the Douglas-Rachford algorithm.\n", "Keep track of the evolution of the $\\ell^1$ norm." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "s = zeros(N)\n", "En_array = zeros(nbiter)\n", "for iter in range(nbiter): # iter goes from 0 to nbiter-1\n", " # put your code here\n", " En_array[iter] = norm(x, ord=1) \n", "x_restored = x" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "\n", "We display the original and the recovered signals." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 15, "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": [ "" ] }, "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": { "deletable": true, "editable": true }, "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": 16, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/condatl/anaconda/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": [ "[]" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFwCAYAAACIBGAMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmYXFWd//H3l4QlYQ9bCElI2EPYdxiQgIAMiIjI4Mq4\nOzqK+nNHlKgzbijq47iMwzCPiuCKyiJgBBpB2UkgK0nIQkI2si90k4R8f3+cc7tv37pV3bV1VXd9\nXs+Tp6rOvbfqdBGST87yvebuiIiIiDTSdo3ugIiIiIgCiYiIiDScAomIiIg0nAKJiIiINJwCiYiI\niDScAomIiIg0XN0CiZmNMrMHzGy6mU0zs6vr9VkiIiLSv1m96pCY2XBguLtPMbNdgKeAN7r7zLp8\noIiIiPRbdRshcfdl7j4lPt8IzARG1OvzREREpP/qkzUkZjYGOB54rC8+T0RERPqXugeSOF3zO+Bj\ncaREREREpJvB9XxzM9se+D1ws7v/Mee4bqQjIiIygLi7VXJdPRe1GvAzYJW7f6LIOV5pxwc6M5vo\n7hMb3Y9mpe+nNH0/pen7KU7fTWn6fkqr5u/1ek7Z/BPwDuAcM5scf11Yx88TERGRfqpuUzbu/jAq\nvCYiIiK9oMDQvNoa3YEm19boDjS5tkZ3oMm1NboDTayt0R1ocm2N7sBAVbc1JL36cK0hERERGTCa\ndQ2JiIiISK8okIiIiEjDKZCIiIhIwymQiIiISMMpkIiIiEjDKZCIiIhIwymQiIiISMMpkIiIiEjD\nKZCIiIhIwymQiIiISMMpkIiIiEjDKZCIiIhIwymQiIiISMM1PJCYMbTRfRAREZHGanggAd7S6A6I\niIhIY5m7N+7DzRx8OXCCO0sa1hERERGpmpm5u1sl1zbDCMmPgJvNGNTojoiIiEhjNEMg+U9CPz7f\n6I6IiIhIYzR8ysbdzYyRwJPAu925u2EdEhERkYr19ykb3FkMvBn4XzMmmrFro/skIiIifacpAgmA\nOw8DpwBHAPPNeFdjeyQiIiJ9ZXCjO5AWR0reYsZ44E4zdnDnp43ul4iIiNRXUwWShDvTzTgfeNCM\n1e78rtF9EhERkfppykAC4M5cMy4CJpmxzp1Jje6TiIiI1EfTrCHJ484zwOXAL804tdH9ERERkfpo\n6kAC4M5DwLuBP5lxZKP7IyIiIrVX10BiZhea2Swzm2Nmn630fdy5C/gkcI8Z+9WuhyIiItIM6lYY\nzcwGAc8B5wEvAk8Ab3X3malzyiqgYsZXgdOAC915tcZdFhERkSo0a2G0U4C57r7A3bcAvwIurfI9\nvwzsAFyTd9CMS824pMrPEBERkT5Wz102BwCLUq8XQ3ULU93ZasZbgafMeMidtuSYGccCPwUGmXGy\nO/Or+SwRERHpO/UcIanLXJA7S4B3Ab8wY2jq0CXALcDPgavq8dkiIiJSH/UcIXkRGJV6PYowStKN\nmU1MvWxz97ae3tide814DPgw8O3Y/M/AdcDW2PblinotIiIivWJmE4AJNXmvOi5qHUxY1PpaYAnw\nOFUuau3+/pwK/AI4PDatA8YAG4BVwBh3Vlf8A4iIiEhZmnJRq7tvBT4C3AvMAH6dDiM18Hh8PAXY\nF9jizmp3tsTPG1/DzxIREZE6qmvpeHe/G7i7Pu+Nm3Ez8E7CDp45qcPTCYHkoXp8toiIiNRW097L\nppduBh4DnqV7IJkGHNWQHomIiEjZmr50fCnuzCOsU/kYMDt1aDqozLyIiEh/0a8DSXQzMBL431Tb\nPMICVxEREekH6rbLplcfXsVq3K73YHtgrHvXCIkZOxF23QxViXkREZG+0ZS7bPqKO1vSYSS2dQBr\nQDfiExER6Q/6fSApYSFwYKM7ISIiIj0byIHkBWB03gEz3m7Gz8w4qI/7JCIiIjkGciBZABxkxnZm\nDDfDAMw4HbgeWAn83owdGthHERERYWAHkgeBS4H7gbnAH8wYAnwH+DzwKUKJ+X9tWA9FREQEGNiB\n5D5CcbRlwN5ABzAZGALc7I4DXwK+oFESERGRxur3235Lvz+XAA+7syZuD74MeDIWVEvOuRe43Z0f\n1qsfIiIiraCqm+YO5EDSuz5wFNAGvMGdfzSyLyIiIv1ZS9chqZY70wjrSP5gxvsb3R8REZFW1PIj\nJAkzDgb+ThgpebzR/REREelvNEJSA+48D3wW+E6yRVhERET6hgJJdzcDw4FTG90RERGRVqJAkhJv\nxPcL4G1JmxlXm7HajBs0ciIiIlIfCiSFbgGuNGOwGccDXwDOBS4kbBsWERGRGlMgyXBnLqHs/GsJ\nFV2/7s4UwvqSTzWwayIiIgOWdtnkMOMDwP8DdgcOc2eDGYOBRcDZ7sxuaAdFRESakHbZ1N6NwF3A\nFe5sAHBnK/An4PXJSWYcZ8aXzRjVmG6KiIgMDBohKYMZVwDvcudiM0YQ7o3zV+BY4AR3Nje0gyIi\nIg2kEZK+8wBwVrwZ38eBX7nzdmAF8OaG9kxERKQfUyApgzsrgbnAWcB7gO/FQz8E3teofomIiPR3\nmrIpkxnXE3bgrHPnnNi2C7AM2M+dTY3sn4iISKNoyqZv/QnYBfhG0uDORsJ6kjMb1SkREZH+TIGk\nTO487M5h7tybOdRGmMrBjCFm3GnGk/GmfSIiIlKCAkntPA0cH59fA7wC/BH4ccN6JCIi0k9oDUmN\nmDEG+AcwGlhMGC1ZCMwB3ujO5NS55k7jvngREZE6aLo1JGZ2vZnNNLNnzOw2M9u9Hp/TZBYCOwFX\nAovdmRPrkvwOeENykhnDgZlmPGvGXo3pqoiISHOp15TNX4Dx7n4sMJtwT5gBLY54tAHfAX6VOnQX\ncHHq9Zdj2wPAN/uqfyIiIs2s7lM2ZnYZcLm7vyPn2ICZsgEw483Ab4ED3Xkhtu0IrAH2jKctAY4j\nrDGZDRygrcIiIjIQNN2UTcZ7gD/3wec0gzuAtyRhBMCdV4DngXHAGcAcdxa5swJ4FLioIT0VERFp\nIoMrvdDMJgHDcw5d4+53xHO+AGx291tKvM/E1Ms2d2+rtE+NFsPHr3MOTQWOIYySPJVqfwA4nTCq\nIiIi0q+Y2QRgQk3eq15TNmb2LuD9wGvdvaPIOQNqyqYYMz4H7APsATzuzn/H9vOAL7nzmtS5g+Od\nhUVERPqVppuyMbMLgU8DlxYLIy1mOmHK5mjCaEniKeB4MwYBmHE5sNGMm81UI0ZERFpHvf7S+wGh\nvPokM5tsZj+q0+f0F88BRwFHAtOSRnfWAGuB0fEOwt8n7Mg5Gq0tERGRFlLxGpJS3P3QerxvPzYf\n2J+woHV95tgi4ID4a5k795nxXeDDwJ3JSWYMJiyKfUq7ckREZKDRtEAfcGcLMA94Iufwi4QwMoFQ\nxwTg98Br4pbhxI+AW4G/JFM8IiIiA4UCSd+ZCTyW0/4iMJJQav5BAHc2EKZ5TgYw4wDgzcARgANv\n6YP+ioiI9BkFkr7zb8BNOe3JCMnBwKxU+0PAmfH5lcAfYlD5IVBQZE5ERKQ/UyDpI+4scydvx9Fi\nwgjJyPg8MYWwCBZCMPlLfH47cKYZu9arryIiIn1NgaTxXgROBNa7055qnwccFJ+PI2wdJi5onRqv\nERERGRAUSBpvGjAWusrNR/OAg+J24DHAnNSxx4FT+qR3IiIifUCBpMHcWQ08Q2EgWQIMA44FFsay\n9InOQGKGmfFTM+aZcUZf9FlERKTWFEiaw/3AgnSDO9ti2/l0Hx0BmAEcHp+/kXA/nM8Dt8R6JSIi\nIv2KAklz+DLw1Zz2+YSRkGU57WPNMMKOm++482tgIXB5+kQzLjLjRrPO9SgiIiJNR4GkCbizPpaR\nz1pKmLJZkTl/HfAKYWfOeXRVdP05YcQEgBhCfgZsAm5XQTUREWlWCiTNbRlhQeuKnGPzCOFjnjsr\nY9tfgdembsz3duBWdz4GdAAX1Le7IiIilVEgaW5L42NeIJkPXAo8mzS4sxBYT1f9ksuB38bnNwJX\n1aebIiIi1VEgaW7J2pG8QPIM4f43UzPtU4CjzdiZsPD1H7F9EuH+OFaHfoqIiFRFgaS5lRoh+SUw\niMJAMp0wQnIoMNedV2P7PMAIU0AiIiJNRYGkuRUNJO4sAD4JPJI5NB0YTxgdeS51vhNGS05P2swY\nlLmjsIiISEMokDS3ZYTdNKvyDrpzgzvrM80zyAkkqWNHAJixF6FK7EIzTqhlp0VERMqlQNLE4n1r\nDnFnaxmXzQUOJEzbZAPJc8Bh8flnCXcU/iLw3ezaEjOONOPkijouIiJSJgWSJufe7Q7AvTm/gzCi\ncg4wK3N4NnBYDB/vBK4H/o9wL53xyUlmjCOElTvNuLjy3ouIiPSOAsnANA/Yl8IRktmEEZJxwGbC\notethMJqr0ud90Xga8C7ga9qZ46IiNSbAsnANB9YHiu6doqvNwJvAx6IC10B7gUuhLDQlVBA7dfA\nPcDuwPF91G8REWlRCiQD03wKR0cSs4G3Ao+l2h4FTowjIScQwszieIO/PwPn1rOzIiIiujPswPQw\nYSQkz2zgLML24MQKwAnTPMcAT6SO/QO4og59FBER6aQRkgHInUnuXF/kcDJyMj11vgMzCWtLxgAL\nUuc/ApxR+16KiIh0USBpPbOBZe4FtU1mELYKHwgsTLUvBHY2Y3cAM7Yz40ozju2T3oqISEtQIGk9\nSe2RrDnAwWRGSOLoyQvA6Nj0GeBa4H4zDqlnR0VEpHUokLQYd1a7c2POoeWENSTZERIIgeRAM3YA\nribs0vkemWBjxmlmTDXj2tr3XEREBjItapXECmAEsD8UFGNbSBghORlY6s5UM1YAs8zY3p0tcYfO\nTcAPgU+b0ebOw33YfxER6ccUSCSxAjiOsOV3c+ZYMmUzmLgDx53lZjxPWPD6IHA0sKM7P4wjKR8E\nBRIREemduk7ZmNknzWybmQ2r5+dITbwEDKNwugbilA1wIvBUqv0+4DXx+XnApPj8ZuCNZgq8IiLS\nO3ULJGY2Cjif/L/gpPm8FB/z/nstAYYTKrZOTrVPJ2wVhlDb5AEAd14CXiR1fxwREZFS6jlCcgNh\nR4b0A/GmfOvpXoMksRLYGxiZOZ7ULoEwgvJ86tgjwGm17qeIiAxMdQkkZnYpsNjdn63H+0vdrCB/\nhGQlYbHrrsDqVPsswt2DtwNG0X0x7OOERbAiIiI9qjiQmNkkM5ua8+sNwOeB69KnV91T6QtL6D7K\nkVgF7AWsive3AcCdDcAawh2EdyMEmsQ8Qk0TAMy4woxZZlxVh36LiEg/V/GiQ3c/P6/dzI4CxgLP\nmBmEYf6nzOwUd1+Rc/7E1Ms2d2+rtE9StTfRfQQEAHdeMWMD3QNHYilhd86SdFgBFhH+22PGUOAn\nwIeBH5pxv3vXaIoZw4AjgEdSdyAWEZEmZ2YTgAk1eS/3+v75b2bzgRPdveAvOjNzd9foST8Qt/jO\nd+e8TPskwvTMa9w5K9W+GyGs7AJcBbzZnUvM+D5hpOUr8bwdCDt39gL+x73byJqIiPQj1fy93heV\nWvUv3oFhJfkjJGsJNUgWpRvdWQ9sBfYkpOc74qHbgEtTp76VEFyOBz5mxt417bWIiPQLdQ8k7n5Q\n3uiI9DulAsk4QqjIWkRY7HoMkCxw/jtwkBn7xNcXAr92ZzlwO/COWnZaRET6B93LRnqrVCAZS7gX\nTtYiwsLWcYSaJbizFZgKHBV357yWroJqf4LuU0IiItIaFEikt24F7slpXwsMIj+srAJOAZbFHTmJ\nacBRhK3E29x5IbY/Apweg4qIiLQQ/cEvveLOPe48nXNoTXzMGyHZQNg9k61tMo2w7mQUdIYR3FlC\nKM52WNUdFhGRfkWBRKq1Nj7mjZBsAA4ghIy06YSy8qPILIYlrDVJqr9ixtFmvM2M7WvTXRERaUa6\n+ZlUKwkkxUZIRgBzMu1LgP0IdUoWZ44ldxbGjNGEG/gtAY4Erq1Nl0VEpNlohESqlQSSl3KObSCs\nE8mOkKwi3Bsnb4RkITGQAP8K/Aa4CPhILKAGgBmDzPhy/KXfxyIi/Zz+IJdqrQXWuvNKzrGNhFG4\nbCBZSyiYNpb8EZID4/O3AjfHtSUPAJekznsv8M+EXTmfqOYHEBGRxlMgkWrNI9y7KE+ys6ZbIIkl\n5lcDx1I4QvICMDqWmx8LPBnbswXVPgx8GvgocLWZph9FRPozBRKpijsd7vykyOENmce0VYTAMT/T\nvpAwQnI4MCfWLYGwJfhEADP2Ag4C/hF3/qwCTq/4hxARkYZTIJF6yh0hiVYCrwDLMu0rCOtLxgMz\nUu3zgGFm7AmcDTzszpZ4rA267qMjIiL9jwKJ1NPG+JgXSFYBCzJ3CMadVwlB5SRSgSSe9yyhDP14\nYErqsoeAM2vXbRER6WsKJFJPPY2QzCty3UZCQbUXMu1TCRVeDwBeTLU/TggwnbTzRkSkf9Ef2lJP\npQLJKgrXj6SvOwBYl2lfAgynMJAsBXY1Y2cAMw4DFpvxlBm7VNh3ERHpQwokUk+lFrX+AvhpietG\nUhhIVgD7kgkkcTonWQwL8B/A9wkl6r+YfXMzdtEIiohIc9EfylJP7cA2ckZI3JnmzjNFrtsA7EFh\nIFlOqPCaHSEBWACMMWNX4HXAjcDXgHeaMSg5yYy3EEZn7tJWYRGR5qFAInXjjgM3EQJAOZLFsHkj\nJCOBPSksVb8AGENY3PqkO6vceY5QQfYMADN2AK4nFFPbHbi8zH6JiEidKJBIXbnzfnc2l3lZMsWT\nF0iOB1bE3ThpCwh1TY4g3LwvcR9dNUpOjtc+BHwL+Pcy+yUiInWiQCLNqNhi2OWE37NP51yT1C85\nAngu1T4FOC4+P5OwRRjgXuAELXoVEWkOCiTSjDYAr+TcHycJKk/kXLMJGEqo8NpjIHGnnRBszqhR\nn0VEpAoKJNKMNlI4XZOsSQF4LOeaTcDOFAaSmcBYM4YAh9G9+utDwD/VosMiIlIdBRJpRhvICSTR\n/u78Jaf9ZcIdhPcl1CsBIJaXX0HYnbMf3UvVPwcckrww4wAzvmLG+Oq6LyIi5VIgkWZUNJC4F9z7\nJrGJEEY6cha8rgBGA0OAtan2eYSb9GGGAX8klKa/M24fFhGRPqJAIs2o1AhJMZsIVVzzirCtAI4G\nlqemfSAVSIBjCYtiLyOsO3lr+g3MONGMb5ixT5n9EhGRXlAgkWY0G3i0zGs2EeqTFAskx1BYu2QZ\nsFvcafNm4LcxsPwC+JfkpFiS/k7CtuLb42iKiIjUkAKJNB13HnPn2jIvezk+lhoh6TbdE0vOzycE\njROBv8VDdwNnxIWwAFcBjxBGTfYEJpTZNxER6YECiQwUm+JjOSMkECq5DiPswJkDnVuCnyOEGICz\ngT/FAHMTqvAqIlJzCiQyUPQ0QrIz5C6I7QB2I5SkT999OF2/5DS6ppDuA86ttrMiItKdAokMCHH0\nop38QJLsrPllzrEOYBywOFPifjJwvBnDCYFlTmyfAgyP7SIiUiN1CyRm9lEzm2lm08zsm/X6HJGU\nl8kPJPcBF7ozM+dYB2FqZk6mfTZwMHAgMDcGHuKW4mmEEvVA2DKsha4iItWpSyAxs3OANwDHuPtR\nwLfr8TkiGZvICSTuvOzOvUWu6QBGENaSpC0n1DXZB1iZOfY8saCaGYOB24GVZpxWeddFRFpbvUZI\nPgR83d23ALh79g97kXrIDSQ9aCcsan05076CEEj2pjCQzCWMngBcGa//MPB/2ZESM7Y3Y6cy+yQi\n0nLqFUgOBV5jZo+aWZuZnVSnzxFJKzZlU0oHsBeFgeQlwujIvhSOnsylq+T85cD/AL8BDDg1OcmM\nvYHpwDyzzh07IiKSY3ClF5rZJMhd2PeF+L57uvtpZnYy4Q/rg3LOFamlSkZIOsgZIXFnsxkbCNuB\nn89c8zxwiBmDgAuA97vjZvwKuIKuHTkfJdQ2mQp8HXh9mX0TEWkZFQcSdz+/2DEz+xBwWzzvCTPb\nZmZ7ufuqnHMnpl62uXtbpX2SlldpINmZwhESCNM24ymsGruCMKoyCljtTvL7+kHgq6nzriKspZoN\nXGvGge4sLLN/IiJNy8wmUKNikRUHkh78kVCr4UEzOwzYIS+MALj7xDr1QVrPD4Bny7ymIz7mBZLl\nhFok2TUkLwNDCdM26dGTycCxceRkGKGq63R3tpkxCTgfuLHM/omINK04iNCWvDaz6yp9r3qtIbkJ\nOMjMpgK3Ev6lKFJX7tztzotlXtYeH4uNkOxG6UAyN/X5a+M1hwInA08m24WhM5CIiEiOuoyQxN01\n76zHe4vUWKkRkkcJN93LLmrtAHYiE0iiqcBRhGJrT6baHwGuqbazIiIDlSq1SqsrFUhuAM4irAHp\nFEc9OgjrS+ZlrllCWOy9P/BCqn0eMMqM7aGzmNpbzHhN1T+BiMgAoEAira5oIHHH3XnYHc+5rp1w\n/5vsdM4yQiDZC1ideq/NwIuEOwsDvBf4CvBbs66twiIirUqBRFpdqRGSUl4G9gPWZ9qXx/ZhQHYh\n9xzCNmKAjwPvBq4lM5Vjxq5mTDTjvDL7JCLSbymQSKsrtai1lJcJVVyz24yTEZJhpEZIojnAoWYc\nTBhBeQT4FXB2LKKW+A5hUeyvzDoDjIjIgKZAIq2umhESozCQpEdIsoFkOaH66wRgkjvb3NkAPBTb\nMGMvQnG1dwPfBT5fZr9ERPolBRJpddUEEiicssldQxJtAHYhVC1O3134CSC5vcKZwKPurAB+DlwS\n65qIiAxoCiTS6qoJJNtyrltOWOw6lMKwspGuQJLenfMkcGJ8fgLwFIA7i+L7nYiIyACnQCKtrppA\nsjG7A8eddmArMChnd84GYFdCIJmfan+aEESIj0+njj1AGDURERnQFEik1bVnHnvrZQpHQBJPF2kv\nNkKyHNjZjKHAMcAzqWMzgMPL7JuISL+jQCKtrgPY7M7WMq9rp/iN/O4r0r6RsL5kZ0IIAUK9E2Ap\noZjaPoR1KInngCOSF2a83ow5ZnykzP6KiDQ1BRJpdesI914q18sUDyQTyR/V2EAoN78sZzpnKXAg\nsD3dp49mJe9lxhDgf4EvAV80Y1z2A8zYufc/gohI81AgkZbmzmZ3PlTBpUWnbNzZ4t693HyUTNlk\nq7tCCCTjgDWZsLKEMJ2zB2Fr8HPu3Ar8H/Ce9BuY8QNgnRmfKPNnERFpOAUSkcqUGiEpZmN8zN6s\nD0IgORJYm26M4WQJobbJRcCf46GfA1cm58UCalcCxwLXmnFAmX0TEWkoBRKRylQSSJLzi42QFASS\naD2wG6F668OxbSZh5GREfP1u4P/cmQ78FnhbmX0TEWkoBRKRyswFppV5TbI2pNxAsoEQSEYDC6Bz\n5ORxQkghPrbF53cSRlNERPoNBRKRCrhzmzvXl3nNNmAT+VM2C4B9gTU5x9YTdt/sRQguiceBU8ww\nwlTNlNj+AHBy3EYsItIvKJCI9K0N5I+QJKXki03ZHAkscefVzDVjgRGAE7cLu7OJUHhN9UtEpN9Q\nIBHpWxvJDyQvEmqiFAskRwMLM+1J7ZJxwPTM7pyZsR0AMw424w1m+n9eRJqT/nAS6VvrgBXZxjid\n8zzFA8lRwAuZ9iWE0ZE9gFWZYzOJBdXM2BN4ELge+EoVfRcRqRsFEpG+9Rbg0SLH5lB8DcnBFB8h\n2Z3CmijpEZL3EarHng18JNY06WTG5WZ81Izte/tDiIjUmgKJSB9yZ24cDcnzeeB3Oe3rAaMwkKwj\nVHbdn8JAshAYFZ9fANzmzjLgXkIoAsCMCcD3gHcAH+/1DyIiUmMKJCJNwp1Z7rnrS5Kw8ULm/OQe\nOIdTGEjWArvHcvOnEXbeQNgSfE7qvA8B3wDeCXzOjB2r+iFERCqkQCLS/JKwkR0hgbCO5AgKA8k6\nwtqSg4DF7p3HHwFOBzBjEHAxcEssdf8cYVpHRKTPKZCINL8kTCzKOZaMkKzLtK8jrC0ZlbnueWBI\nLC1/MLDCvXPdyp3A62vVaRGRciiQiDS/9cBL7t3uApxYCuxK4QjJy4T1JQeTCiRxmmcqYcHrMcCz\nqWseBk6qXbdFRHpvcKM7ICI9eg64tsixpHJrt0DijpuxDhgPLM5cswgYSSiqNjXVPgMYb4ZlapqI\niNSdRkhEmpw7m9z5aZHDS+JjdoQEwsLWoymc6llMmMo5BJid+pyVhOJsIwDM2N2Me834tRa7iki9\nKZCI9G+5IyTROkJBtWwgWUQIJLtSuPZkBmFUBeDLhHL0Q4HPZN88LooVEamJugQSMzvFzB43s8lm\n9oSZndzzVSJSgSSQZIMFhBGSPYB5mfbFhCmbXQil7NPmAmPjDfveBlwHfAH4oFnXFK8ZZwGrzbhb\nwUREaqFeIyTfAr7o7scDX4qvRaT2So2QJEHh+Ux7eoQkG0hWA3sSRknWu7PAnWcJ999J/8PiekIN\nkyHA+yvuvYhIVK9AspSw5RDCv9BerNPniLS6VcDthLsIZ42FzvvkpL0IHEAYIcletwYYRqhH0pZq\nfxg4A8CMI4DRwG+A/wDeXc0PICIC9dtl8zngYTP7NiH0nF6nzxFpaXE3zKVFDu8OuVuF18Rj7RSO\nkKwhbBUeTNjdk/g7cAXwHeBM4C/ubDXjQeBQM0a4dy6wFREpW8UjJGY2ycym5vx6A/C/wNXuPhr4\nBHBTrTosIr12KmGXTTfuvEoIIvtRfMpmBN1HNp8ATojPxwPT43ttAe6nezl6EZGyVTxC4u7nFztm\nZje7+3nx5e+AG0ucOzH1ss3d2yrtk4h0ce82wpG1BtgN2JTTviewA90DySJgRFzAOh6YlDr2DGE3\nj4i0GDObAEyoxXvVa8pmrpmd7e4PAueSqnWQ5e4T69QHESluDbC/O5tz2vckrP3qnIJx5xUzVhJG\nTsYTtgcnpgHvSV7EmiXnAY+6s6o+3ReRZhAHEdqS12Z2XaXvVa9A8gHgh2a2I2Ge+gN1+hwRqcxa\nCqdroCuQ7E/hYvSFhJv17Uf3Ow9Po/sIyXeBC4AOM05yp6NWnRaRgasuu2zc/Ul3P9Xdj3P30919\ncj0+R0Qqtob8nTlrgDFAe869c14g3ANnU2bnzjzCdM6O8aZ9VwKnEALNlek3MGOUGZ80Y//a/Bgi\nMlCoUqunpv9tAAAgAElEQVRIa1pD/ghJUmAte/8bCCMk48js3ImLZNcAexEWtz7gzmrgv4D3JefF\nwmr3ARcCvzfTnz8i0kV/IIi0ptxAkhr5uCXnmhcJ97/J20q8iq76JQ/GtvuAE8wYEl+fS5gqeh2w\nIyGYiIgACiQirarYCAnAZYR6I1kbgH0pHkj2Ak4j1CwhTvlMJWw/hjB9c0sMPb8ALq+08yIy8CiQ\niLSmYotaceePObtviOcXCySrCYFkJGFqJ/F3QkgBOBZ4JD7/I3BJvGeOiIgCiUiLmg48XeY1pQLJ\nKsKW4KGEcJKYAxwUg8fhxOqv7iwAnFDCXkREgUSkFbnT5s5XyrxsI7ATxQPJkcDyWM4+sYCwa2ck\nsNGdtaljzwDHlNkHERmgFEhEpLeSKZ5iUzZH0XX34cQCQiA5ApiVOfYsqUBixufMeMSssNy9iAx8\nCiQi0lulAskqQgXXZZn2hYQ7Ax9I97UlEAqqjQcw41jg44RbTfxSa0tEWo8CiYj0Vk+BZBiZERJ3\n2gkLaA/KuW4pYU0KhNLz/wXcABhhizAAZuxkxj1mPG/G4dX+ECLSnBRIRKS3SgWSpBrzlpxjSwgj\nJNkS8iuBfeLzM4H74/qTW4DXp857d3z8EfDjMvssIv2EAomI9FZyZ+CCQOLOfOBEwghHVjvh/jiv\nZNpfAvY2Y1fCDpynYvuDhAJribcQwsj3gaPNGFNh/0WkiSmQiEivxBLx7eSPkODO0zGYZHUQ7h5c\nbITkKGCme2dgeRI41Iw9Ylg5EbjXna3AbcAVVf8wItJ0FEhEpBwbKRJISmgnJ5DESq5OWF+yJNW+\nmVC/5ND4a24qrPwVOKOinotIU1MgEZFyVBJIio2QQJi2GQ+syLTPIwSVw4DZqfYpwHFlfr6I9AMK\nJCJSjmpGSLJrSCBM24wHlmfa5wFjKQwkzxPWneyZNJgx3oxxZfZJRJqMAomIlKPSEZIhFB8hOZLS\nIyRzksZ4Y76pxIJqZpxBWAT7NzNOL7NfItJEFEhEpBwzgcVlXtMeH/MCyXLgEIoHkpHAC5ljLxDu\nmwPwJeCzwOeA69InmbGdGVeZcUqZ/RWRBlAgEZFec+e97jxR5mUdmce0x+JjdspmCbA/4d457Zlj\nK4B9zRgCnAXcCtwMnGzWGVQgBJVPAXeZcUSZfRaRPqZAIiL1lgSRvDUkD8bH7AhJB7AjIZBkr1sB\n7AecAMxw5+W4C+dB4BwIoyPAR4Arge8QwomINDEFEhGpt1JTNjPi45JMewchjOyYc90KQsn504FH\nU+3301Vy/hRgrTszgZ8Dl5qxU0W9F5E+oUAiIvVWdMrGHXfH3FmTOfQKXYEkO0KynBBIDifcoC/x\nOHB8fH4moWYJ7iwh3GlY9UtEmpgCiYjUW6kRkmLSUzbFRkj2JuzSScwFDo53Cj6E7tuFnyRM8YhI\nk1IgEZF6K7WGpNQ1xUZIkkCyD6lA4s5qYBshqBxMqFmSeAoFEpGmpkAiIvVWyQjJVmAQsDPFF7Xu\nQyisljaXEEYOoXsgeZpUIDHjBDP+asZ7y+iTiNSRAomI1Fupbb+53HG6pm2y120EtgBj6D5lAyGE\nHEmoU7Iw1T4XGGOGxSmdnxB25XxdVV5FmoMCiYjUWyUjJMn52+JdfjvFsDIXGAwFi2EXEe4OvCLe\npC+5pp2ue+qcEh//E/g+8P/Sb2DGe8xYpNETkb6lQCIi9VbJGpLk/GIhZi6w2p1XM+0vE9aQbMq5\nZhmh2NoZwF9jGfpfAxfHURPM2Bv4NvDvwDfNGFVmn0WkQgokIlJvyQhJuYGko8Q1cyicroEQSPYi\n/347S4HhwMmELcIQpng2Q+e0zcXAA+7cDvweeHuZfRaRClUcSMzsCjObbmavmtkJmWOfN7M5ZjbL\nzC6ovpsi0o91AFuzUy+9vK5YIJlL4YJWKB1IkhGSzkASp38eJNQtAbgEuCM+/yWh0quI9IHBVVw7\nFbgM+O90o5kdSfif+EjgAOCvZnaYu2+r4rNEpP/qoPz1I1B6yuYeYFVOezthyiYvrCwlLHYdQ/cd\nOLOAQ+PzY4AvxuePAYebMdS97Dsci0iZKh4hcfdZ7j4759ClwK3uvsXdFxD+JaO7bYq0rnbKn66B\nEiMk7ix3586cQz2NkBwKdMR73yTmAIfG+9+MBhbEz3gFmE5X9VcRqaN6rCEZQffbky8mjJSISGta\nCnyggutKTdkU8zIwhOJrSMYDqzPtcwhBZQRhoWz67sJPEKZ4ADBjqBmnJotgRaR2Sk7ZmNkkwiKw\nrGvc/Y6c9mK8rF6JyIARd8LcVsGlpaZsimnPPKYtIwSSeZn2ucBBhGJq8zPHZsRrMGMQYRpnH+B/\n6JraEZEaKBlI3P38Ct7zRei2VW5kbMtlZhNTL9vcva2CzxSRgafSEZL0Y9pSYHcyIyTubDJjA2Fq\nORtIlgDnxeeXAeuBC4DpZvw43rgPCNVfCetT/hAXy4oMeGY2AZhQi/eqZlFrWnr48nbgFjO7gTBV\ncyhdW+wKuPvEGvVBRAaWShbDlgoky+Jj3mLYDYStvy9k2pcSduYAXAHc6M5SM/4IXA78AMCMg4B7\nCYtpDyLUMhEZ8OIgQlvy2syuq/S9qtn2e5mZLQJOA+4ys7tj52YAvyEMdd4NfNjd9a8FESnXK5Q/\nQlJqymY1oeR8dg0JhHL0I4B1mfYldAWSEwlTNhB2+VyYOu9DwI2EbcOfM2OnMvst0vKq2WXzB3cf\n5e5D3H24u/9z6tjX3P0Qdz/C3e+tTVdFpMXUdMomVmZdTv4IySbCerkNmfZlwHAz9iTc0O+52P5X\n4DVmnaPMFxGmauYCzxCCiYiUQZVaRaRZ1XrKBsIUTLERkuHxsVPc+ruBsI7k2aRUvTur4/scaMYI\nQlh5Ml52O3Bumf0WaXkKJCLSrGo9ZQMhkBQbIdmHwhESCNM2J1G4vuQ54HDCOrlZcQQGwpq5U8vo\ns4hQu0WtIiK1VsmUTRJEio2QfIowbZO1kbA4f2POsWXAERTuFkwCyUpgYap9MnCEGUMyNU1EpASN\nkIhIsyp7yibeL2czRQKJO3PcWZ9zKAkieSMkqwg1StZm2pNAciCpQOJOB6G2yTgAM3Y047tm/MCM\nIWX8OCItRYFERJrVncBvK7iuneJTNsVsio95IySrCFt512TanwfGkgkk0SK6KlS/l1B+/jDg42X2\nS6RlKJCISFNy5xl3Hqng0pcpPmVTTE8jJDtROEKyDtiNEEiy60tepCuQfBiYCHwU+GRqZw5mjDTj\ncTN+ESvBirQsBRIRGWgqCSQ9jZBA4QjJemBXQmXq3EBixr6EYPKgO7PjeaelzvsecH98j38rs88i\nA4oCiYgMNL+jcAqlJz2NkEDhCMkGQiAZRljYmpaMkBwLTEmVkv8zoWYJZgwFXgd8DfgK8EHdtE9a\nmQKJiAwo7nzOvSAg9GQTsMWdzTnHio2QbCBM2exBYYXXJJAcB0xJtf+DUPEV4CxCWFlPKL09jLCF\nWKQlKZCIiIQRkrzRESg9QrIHsC3urElLAskxwLOp9pnE3TfAmcR7gMQaJn+n+3SOSEtRIBERCYEk\nb/0IdFV27RZI3NlC2JacDSoQap3sG38tSbW/AOxtxi6ExbDzUsceAU4vu+ciA4QCiYhImLLpaYQk\nOy0DYWFrXiBZD+xOmNLprHsSS8/PJtQvGUnYHpx4HDg5/SZmHG7G9r3ov0i/p0AiIgJTgW8VObYW\neGMsupa1gfxA0kGo/LovFBRim0MotDaK7oFkHjA6eWHGR4EngMcUSqQVKJCISMtzZ707Py9yzN35\nU5FLc0dI4q6adYRRkGwgWU1YwDoSWJxqXwHsbsZOZuwAfIFwT5zVwHvSb2DGLma8wYzdevzhRPoJ\nBRIRkcptoHD3TWIdsAOFUz1rCZVf2907658kC1uXEBbDXgDMcWcmcAPw9sx7/Az4DnCXmf4cl4FB\nv5FFRCpXbMoGukZGsotl1wFH0X26JrGIMJVzIvBQbLsfOM6MYQBmHENYa3IMsDMhvIj0ewokIiKV\nK7aoFULw2BBHPrLthxPuIpyVBJKjCetakpv1/YNQtwTgfODOeCfhHwPvq+YHEGkWCiQiIpUrNUKy\njsL1I0n7gcBLOccWEwLJUcRAEk0FjozPzwXui8/vBs5WhVcZCBRIREQq9xNCqfo868jfKryW8Gdv\nXiBZCQwHxgDPpdqn0xVITgQeBXBnMeG+ParwKv2eAomISIXcmezerbhZ2nqKj5BA4f1vINRDGUGY\n6tmSap8BjI+7b4bRfbrn78AZZXVcpAkpkIiI1EepKRvIHyHZCOxH4ULYWcARhLCyLBZYS8wkNUJi\nxllm/JcZYyvtuEgjKJCIiNRHsSmbUoFkE2HKplsgiTfg244QPBZnrnmBWFDNjD2BPwFDgV9qS7D0\nJ/rNKiJSH0vJ30mTLIItFkjyRkgg1Ds5inDjvrTOQAJcCUwi7LzZBZiQPtGML5hxhxn796L/In1K\ngUREpD5+CXw8p30D4BSfstmV0oGk6AgJ8EbgV3Gr8W+AS5OTzDgL+CBhsexNvf4pRPqIAomISB3E\nkvPZGiRJRdbZhBGUrKRyazkjJIuBEWYMAsYDk2P77cDFqfPeC1wPXAucYKadOdJcFEhERPqYO0fE\ndSFZpQLJWkLg6DZC4s4rhPvdHA7sSRgxgbBVeH8zdo2vzwb+Ggut/ZYwmiLSNBRIRESax8bMY9oa\nQqn47AgJhEByOjArGZWJO3FmAkeaMYqwpmRWPP9vdFV+FWkKCiQiIs2jpykbKFxDAtABHE9X4EhM\nI0zzjAcmx7sQQ7hPzj9pF440E/1mFBFpHr0JJEtyjrUTapRki61NJwSS4aTWrLizlBBiRiZtZphK\n0EsjVRVIzOwKM5tuZq+a2Ymp9vPN7EkzezY+nlN9V0VEBjZ3tgKv0BVM0tYAL8U1I1kdhPUj7Zn2\nRYQwMpzCLcgL6apfsgfwCDDDjAMq/gFEqlDtCMlU4DLCfKSn2l8CXu/uxwD/Cvyiys8REWkVmyi+\nqDVvugZCIBlGuK9N2mpgL4oHkgPj848A8wlF1b6ZfXMzdjJj+950XqRSg6u52N1nAZhZtn1K6uUM\nYIiZbe/uWxARkVKKBZIFhDUhedrJHyFZQwgqw4HHMscWAgfGaZq3Ah8g/CNzsRm7JbuAzDgOaIvt\nE9xz78EjUrW+WENyOfCUwoiISK9sJCeQuPOAO1cVuabYlM1qugJJsRGSEcA+wCMxhDxE9/olXwW+\nSLjD8MfK+klEytDjCImZTSL8Zs66xt3v6OHa8cA3gPNLnDMx9bLN3dt66pOIyABWbISklA7Ctt5i\ngWQz+YHkUsJN+2amirjdS6hZcqsZOwPnAv8CjAIeNmNi5uZ+0sLMbAKZWxRUqsdA4u5Fw0QpZjYS\nuA14p7vPL/H+Eyt5fxGRAer7wJQez+quPfOYWEcIKgcAyzPHVhLCymGEcvKJZwkBBOAk4Fl32oHZ\nZqwGjq6gfzJAxUGEtuS1mV1X6XvVcsqmcyGJme0B3AV81t0fqeFniIgMaO7cXME6jY742C2QxFGP\ndYQRkjWZa7YA2xMqvM5OtU8HjoprS84A/pE69hAqqCZ1Uu2238vMbBFwGnCXmd0dD30EOBi4zswm\nx197V9lXERHJlxtIotXA9FRRtMRWwij54aRGSNx5ia4aJeMJIyaJh4F/qlGfRbqpdpfNH4A/5LT/\nB/Af1by3iIj0WrEpG4iBJKc9GSHZm8LpnOeAQwiLXdPHZgEfTZ9oxlD3gu3GImVTpVYRkf6vpxGS\nGTntSSDZIT5PW0moX7IPoa5UIl27BDPeAKw143ZVeZVqKZCIiPR/pQLJ98gZySZM2Wwff2UDSVJQ\nLRtIVgC7xN03EIqoXUxYNPvm7AeYMUr3y5He0m8UEZH+r+iUjTv3uudWeN1CmLbPCySryAkkcZHs\nImC0GQcTap/cRwg9b0+/gRkfJVR//UEFP4+0IAUSEZH+LxkhKWctRzJlsz1hF07aasLUzNa45Tct\nmbY5H7gnhpS7gHPNGAoQR1C+SNjw8CYzxpfRL2lRCiQiIv1fqSmbYkpN2awiFEx7KXsR8AKhSNrh\nxFL27qwmLIQ9Pp5zATDFnSeBnwD/Vka/pEUpkIiI9H+VBJJSUzarCYEjL5BsBIYCBwHzUu3TgXHx\n+UmEuwdDGD2ZUEa/pEUpkIiI9H/thDuuZ6deSklP2eSNkOxHfiDZTNiZkxdIjozPTwCeis+fAcaY\nsUcZfZMWpEAiItL/dQDtOcXPSulplw3AkznXbQZ2BMYSFq0mZtA9kDwN4M6W+D6nltE3aUEKJCIi\n/V8H5U3XkLpB3o4UDyR/yrl0M2Gbb4c761LtzwGHm7ED4R45L6aOzSDcM0ekKAUSEZH+r50yA0mU\nrCPJBpIVwC3A5JxrkkCyIueavYA9gLWZ0ZoFwJjkhRmjzbhMNUokTb8ZRET6vwXApyu4bivgqdES\nANzZ6s7bi0wBbQZ2BV7JtG8AhhBql6zN6d+Y1OsfAP9D2BosAiiQiIj0e+5sdudXFVy6hcLRkZ7k\nBpIYXtYQFrsWDSRmjANOJtQoudqM/crttAxMCiQiIq2rZoEkWk240/uaTPsCukZIrgR+7c5c4Hbg\nHclJZmxnxt1mLDXTVuFWo0AiItK6tlL7QJI3QrIS2DlWcr2Yrnvr/JLu98A5ExgN/DvwU7Pq7kgv\n/YsCiYhI66pkhGQLsAulA0m3EZI4nbOJUFBtFDA3HpoCjEvdKfgtwM/duY1QC+X8Mvsm/ZgCiYhI\n66p0yqZUIDmYwhESCLuAdiHsxEl26KyKj8Pi49HA4/H5/cDpZfZN+jEFEhGR1lXplI1R3pQNhFop\nBwKr3dkKnSMnc4FD4jmHAnPi80cJC1+lRSiQiIi0rkpHSCA/kKwilJXPLmqFMEIyFliaaZ8LHGLG\nrsBuwJLY/hhwSmo6RwY4BRIRkda1hfLufwOlA8ni+Lg+51gH+YHkecI0zyHA8+5si+0vAYMIIUVa\ngFYwi4i0rq1Q9ghEqUDyc2An4MGcY8UCyUrCVM7BdC12xR03YzEwErqVqJcBSiMkIiKtq5opm4KR\nFXdedefH7t3uY5NIpmyWZdrXEsrNDyOEk7QkkABgxu5mXGPGoDL7LP2ARkhERFpXuWEESo+QlNJB\nWPC6OtO+BtgT2J3CkZDOQBLXmPwZOJYwkvKbMj9fmpxGSEREWlelu2ygskCyD+GeN2nJCEleIFlE\n1wjJhwmjK+8G3pd9czP2N2OPMvskTUSBRESkddV6l00p7YQFqqUCSXYxbHrK5hDgL4TtwEenTzJj\nN+BvwJ2q7tp/KZCIiLSuvgwkHfExG0iSKZvdKBwhWQYMj8/HAvMJIWWIGXunzruEUL9kR+C8Mvsl\nTUKBRESkdfXllE17fNyYaS81ZZMcgxBIFsRiatOAo1LnHQc8TFhXcmmZ/ZImoUAiItK6mmGEZAPh\nHjd7UThlswbYM+6qGQksjO3PEha3Jo4n3BfnduCiMvslTUKBRESkdTU8kMRCaOsJd/ktNkJyAPCS\ne+dnPg6cChAruR4HTAZmA0PNOKDMvkkTqDiQmNkVZjbdzF41sxNyjo82s41m9snquigiInXSiCmb\n7AgJhOBRKpDsCyxPtafvczOBsNZkWZzO0T1w+qlqRkimApcRVjbnuQG4q4r3FxGR+qpkhCQ5v9IR\nkuwaEuiqTZINJBsJlV/3yhybTZjK2Q/4APDfMYxAuAfOqWX2TZpAxduj3H0WgFlh1WEzeyMwD9hU\ncc9ERKTeyi4d786rZrxKbQPJo8CJZNaQxPLxyejJ+lT7NjMuIoy2fIOutSUAM4F3ltk3aQI1X0Ni\nZrsAnwEm1vq9RUSkpioZIYEwbVPJlM2m1M3z0u4FcO8MLWkFgSSe+4g7L7vzjDtrU4fmE3bk9MiM\n4zLbh6WBSo6QmNkkuvaAp13j7ncUuWwi8F13f9nyhk9ERKRZbKH8m+tBZYGkg/zREQiB5Ooix9YQ\nbr6XdwfhPPOBsWZYahqnmMnAfwP/1sv3ljoqGUjc/fwK3vMU4HIz+xZhMdI2M2t39x/lnWxmE1Mv\n29y9rYLPFBGR8m2FHv/SzlPpCEneglbc2Qz8oMh1awmB5O+9+RB31pixjVBsLXvfnE6piq5revO+\nks/MJhAWFletViV2OxO2u7+ms9HsOmBDsTASz59Yoz6IiEh5ttB3gaSDIoGkB2uBk4C7y7hmPnCm\nGXcAo4DlyZZhM/YEBtFVx2QHM8bGtg2EpQzr3cMaSDN2dmeTGTvG91rrXnBXYuK5wwDcu4JQnBLa\nBViYjNgkfSj2Ppn33I5wU8JN7ixN+pM6PgZw927raDBj5/gz7eLOEjNGE3YiFdyluRru3mbGtPCc\nlfHv/YpUs+33MjNbRNhedZeZlfObRUREGm8aML2C6+6j+zbc3pgPPFjBZ60kBIHsDpxSJhF2em5P\nKJaWrur6CeA/CQFnIfB6wndwD2H3zncJO0gxY3dgUSzMNjqe8/4Sn/sFCm/8915gFnBMqu3jwEd7\n+bPsGD/3i2aMJ+wiSrsp/uoUA9ZU4GbgxXin5L8Q7gdUD08DT1T7JtXssvkD8Icezvlype8vIiL1\n5c6tFV73rxVcM5cQBsq1kPAv/d6uIcGdzxA2V2DWucj1qXh4PvAO4Enga8AdwC3uvN2MlYRlB8n0\n0UGEqZ/93ZlDz3+hzwfGZ/ryTTNOAw4Fnkmdd24vf5b25HPNeBNwmBnbJYuD3XPf57D4MyeDDgcR\npr0W5pxblThyNBJwM7av5r1UqVVERJrZ/PjY60CSsYDuu26S1/OBlzKfkYSXBfH12Mxjbz5rTE57\ndudPtk+9NZYw6jOiF+dBCCEPAKeTmoaqsdGEn2cxYSSrYgokIiLSzKoNJPPpHhLS4WNlTtsrwNL4\nutxAUmzLcba911uTM3rbn+T4CsJ01Dl0/Yy1loS7Sn+mTgokIiLSzGoRSNJ/US4GXiX8qz4ZIVmQ\nelyYqpUylrB2ZUwvP2sBcGBciNpTH/aJ0x3l6G1/xsTzFsRf59D1M9bamPjeCiQiIjKgrSRU/a4m\nkJxoxhUA7mwFXojtGwg7htIjJOmRhLGEhbhvMuu5HH2cEllPYf2ubn9Zu/MqIZQcmD7JjCFmfCv+\nOinnI5L+vNeMb2SufW+87lLgjHhe8vPsQ84IiRlXm/VcxNSMD5nxtfh8ePycS8wYB1xFjUZIarXt\nV0REpOZi+fh/ofJ/4c8i7KpJF2X7IPCP+N5vSr33bYQdI4lbgBnAeZBbRTbP56Fga+184L8ybRMp\n3Aa9ja5ppLxt1T8EHgZeR+F27XXx2o2EHUZ3ELYbLwA+FV9nrct5nzxrgSHx+dbU52wh7GL6LTCM\nsAD4ml68Xy5zr2QLem2Ymbu7qrmKiIgMANX8va4pGxEREWk4BRIRERFpOAUSERERaTgFEhEREWk4\nBRIRERFpOAUSERERaTgFEhEREWk4BRIRERFpOAUSERERaTgFEhEREWk4BRIRERFpOAUSERERaTgF\nEhEREWk4BRIRERFpOAUSERERaTgFEhEREWk4BRIRERFpOAUSERERaTgFEhEREWk4BRIRERFpOAUS\nERERaTgFEhEREWm4igOJmV1hZtPN7FUzOyFz7Bgze8TMppnZs2a2Y/VdFRERkYGqmhGSqcBlwN/S\njWY2GPgF8AF3Pwo4G9hSxee0JDOb0Og+NDN9P6Xp+ylN309x+m5K0/dTPxUHEnef5e6zcw5dADzr\n7lPjeWvcfVuln9PCJjS6A01uQqM70OQmNLoDTW5CozvQxCY0ugNNbkKjOzBQ1WMNyaGAm9k9ZvaU\nmX26Dp8hIiIiA8jgUgfNbBIwPOfQNe5+R5HLtgfOBE4C2oH7zOwpd7+/qp6KiIjIgGXuXt0bmD0A\nfNLdn46vrwT+2d3fFV9fC3S4+7dzrq3uw0VERKSpuLtVcl3JEZIypD/8XuAzZjaEsJj1bOCGvIsq\n7bSIiIgMLNVs+73MzBYBpwF3mdndAO6+lhBAngAmA0+5+9216KyIiIgMTFVP2YiIiIhUqyGVWs3s\nQjObZWZzzOyzjehDo5nZTWa23MymptqGmdkkM5ttZn8xsz1Sxz4fv69ZZnZBY3rdN8xslJk9EAvv\nTTOzq2O7vh/AzHYys8fMbIqZzTCzr8d2fT8pZjbIzCab2R3xtb6fyMwWxKKVk83s8dim7wcwsz3M\n7HdmNjP+/3WqvpvAzA6Pv2eSX+vM7OqafT/u3qe/gEHAXGAMYUfOFGBcX/ej0b+As4Djgamptm8B\nn4nPPwt8Iz4/Mn5P28fvbS6wXaN/hjp+N8OB4+LzXYDngHH6frp9R0Pj42DgUcLONn0/3b+j/wf8\nErg9vtb30/XdzAeGZdr0/YSf92fAe+LzwcDu+m5yv6ftgKXAqFp9P40YITkFmOvuC9x9C/Ar4NIG\n9KOh3P0hYE2m+Q2E/xmIj2+Mzy8FbnX3Le6+gPAf9ZS+6GcjuPsyd58Sn28EZgIHoO+nk7u/HJ/u\nQAj5a9D308nMRgIXATfStej+/7d39z42RHEYx7+PsIUlRCReN7EFnWB1WEIIm6BFIaJQKVSbsP8E\nlUZQbKHwGqIRofYSu14jIlGst93ViGhE/BTn3Luzl4jE5cjO80lu9syZW8w8d7Lzu2dmznU+k7U+\nVFD7fCTNAXoj4gxARHyNiI84m5/ZSjqXj9CmfEoUJEuAkcry69xnsCAiRnN7FFiQ24tJOTXUJjNJ\ny0gjSXdwPk2SpkkaJuVwOyKe4nyqjgP9QHWWaOczIYCbku5LOpT7nA90A+OSzkp6IOmUpE6czc/s\nBc7ldlvyKVGQ+C7a3xBpvOtXWU35HCXNAi4CRyLiU3Vd3fOJiG8RsRpYCmyUtLllfW3zkbQTGIuI\nIX4cBQDqnU+2PiLWAH3AYUm91ZU1zmc60AOcjIge4DNwtPqGGmfTJKkD2AWcb133J/mUKEjekK45\nNcKxtZQAAAGcSURBVHQxuYKqs1FJCwEkLQLGcn9rZktz35QlaQapGBmMiCu52/m0yMPJ14G1OJ+G\ndcBuSa9I3+C2SBrE+TRFxLv8dxy4TBpGdz7pXPQ6Iu7l5QukAuW9s5mkjzSlx3hebsuxU6IguQ8s\nl7QsV1l7gKsFtuN/dBU4kNsHgCuV/r2SOiR1k34v6G6B7fsnJAk4DTyLiBOVVc4HkDS/cRe70gSE\n20hz/jgfICIGIqIrIrpJw8q3ImI/zgcASTMlzc7tTtIPoj7G+RAR74ERSSty11bgKXCNmmfTYh8T\nl2ugXcdOobtz+0hPTrwEjpXYhtKv/GG+Bb6Q7qk5CMwDbgIvgBvA3Mr7B3Jez4Htpbf/L2ezgXTt\nf5h0oh0Cdjif5r6uBB7kfB4B/bnf+fyY1SYmnrJxPmlfu/OxMww8afwPdj7NfV1FmtjzIXCJ9JSN\ns5nY307gAzC70teWfDwxmpmZmRVXZGI0MzMzsyoXJGZmZlacCxIzMzMrzgWJmZmZFeeCxMzMzIpz\nQWJmZmbFuSAxMzOz4lyQmJmZWXHfAcKv4b9WXDVtAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(log10(En_array-En_array.min()))" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "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": { "deletable": true, "editable": true }, "source": [ "__Exercise__ \n", "\n", "Test the recovery of a less sparse signal.\n", "What do you observe?" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true, "deletable": true, "editable": 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" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6YAAAG0CAYAAAAhAzpSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X+QHWd95/vP17KwJCR+eNlrgxnhRTIEGDaBJAaRBeQN\nWJpDAkkqMfHu3r1Foolrs2hcldSFTSDXzr347rKVzY1GENal/LjO7lYMtXCzzjIaOWwymMqVCU7A\nWMFOLAVfjX/gkBCCBBKW8Pf+0ef49PR0n+7Tv7vn/aqamjl9nvP0t59+zunznX76aXN3AQAAAADQ\nlEuaDgAAAAAAsLGRmAIAAAAAGkViCgAAAABoFIkpAAAAAKBRJKYAAAAAgEaRmAIAAAAAGkViCgAA\nAABo1KVNBwBgPTObkWSSrnL3403HAwAAAFSJM6ZAO10i6UWS3tp0IAAAAEDVSEyBdnqrpA9I+mbT\ngQAAAABVM3dvOgYAEWb2jyQ9LenF7v7HTccDAAAAVIkzpkA7Pd10AAAAAEBdSEyB9tosabbpIAAA\nAICqMSsv0DAze7uk70h6o6QHJO2X9LCkH5X0PxoMDQCAzks4zt7m7g81GhiANbjGFGiQme2U9Cx3\nP2lmfyrpByX9E0l/K+lRSS9y9882GSMAAF014Tj7h+7+rZpieKWkGXc/Vsf6gK4iMQVawMyukPRR\nd9/bdCwAAPRN0eOsme1y91M5X7tJ0n929xvzvB7YKLjGFGiQmX2XmX23pIGke4bLBs1GBQBAP5Rx\nnB3OlP+OvDG4+3cknc/7emCj4BpToFnXS9oh6QlJW8zsRxQM4QUAAMXFHmfN7KcVXHP695J+XtIu\nSd8r6fkKZsa/X9Lr3f1XJV0r6fuHCe43Je1RcMnNt4avG9XzYUnfN6zjDyVdGK7/85K21bGxQJeR\nmAINcvfFpmMAAKCvko6zZvZVSa90958fPv5FSb+uIIl9rYIzrEvD4vdK+i53v9/MflvSAXf/jpl9\nVNJ7RvWY2e8oSE53SJqR9M8k/St3P2dmP1PdVgL9wFBeAAA6xsx+y8yeNLMHEp7fa2Z/b2afH/68\nv+4YgZYzSd8IPXZJp9z9U5L+WNJ/kfSu0HMys1dJerakZw2XP3v4e1TP06E67h+u4ztVbQDQN0x+\nBABAx5jZGyWdlfQ77v7qmOf3Svo5d3973bEBXWBm/1LST0j6GXd/wsx2SfpxBcNuv09Bcnqlu3/U\nzDZL+veSPi7pSUnXSfqKgqG6/3BUj4LhuqM6/kLSFkn7FNyi5oOSbnD3R+raRqBrSEwBAOggM7ta\n0u9PSEx/3t1/uOawAADIhaG8AAD0j0t6g5ndb2ZLw/soAgDQWrVNfmRmnJoFAJTK3a3pGFrqzyTN\nuPu3zGxO0u9Jelm0EMdmAEDZ8h6ba52Vty9fIMzsVne/tek4iurLdkhsS1v1ZVv6sh1S77aFpCqB\nu58J/X3UzH7dzC5396/FlOXY3CJ92Q6JbWmrvmxLX7ZD6t225D42M5QXAICeMbMrzMyGf1+rYE6J\ndUkpAABtwX1MAQDoGDP7XUlvlvQCM1uVdIukzZLk7rcrmBn0X5nZRUnfkvSTTcUKAEAWJKb5rDQd\nQElWmg6gRCtNB1CilaYDKNFK0wGUZKXpAEq00nQAKM7db0x5/sOSPlxTOG2x0nQAJVlpOoASrTQd\nQIlWmg6gRCtNB1CSlaYDKNFK0wG0QW23izEz78t1LACA5nFcKY42BACUqchxhWtMAQAAAACNIjEF\nAAAAADSKxBQAAAAA0CgSUwAAAABAo0hMAQAAAACNIjEFAAAAADSKxBQAAAAA0CgSUwAAAABAo0hM\nAQAAAACNIjEFAAAAADSKxBQAAAAA0CgSUwAAAABAo0hMAQAAAACNIjEFAAAAADSKxBQAAAAA0CgS\nUwAAAABAo0hMAQAAAACNIjEFAAAAADSKxBQAAAAA0CgSUwAAAABAo0hMAQAAAACNIjEFAAAAADSK\nxBQAAAAA0CgSUwAAAABAo0hMAQAAAACNSk1Mzey3zOxJM3tgQplFM3vYzO43s9eUGyIAAAAAoM+y\nnDH9bUn7k540s4Gk3e5+jaSfkfSRkmIDAAAAAGwAl6YVcPfPmNnVE4q8XdIdw7KfNbPnmdkV7v5k\nnoDMZgfSzIK0Y4v0le3SFpOed0Y6c15aXXQ/sZSn3uR1lFfvuO7L/3fpOS+VLpN05svS47+Up/4q\n48xT/zTl08qOn//Oi6RNL5S2PiGde3xUrkg7Fmm3uNcGzyTHmqXe9HUl9/XsbVlNP5l+W4q3edZt\nLyO+ouuY3Gfq2Scpfem4NLOnyf4BAACQyt1TfyRdLemBhOd+X9IbQo8/Jel7Y8p5+npeNZAOPCy5\nS5926Rc9+Hv0c+Bh6VWDLDFnW0d59Y7r/tHH1sf9rsenrb/KOPPUP035tLLj5xP38S1527FIu8W/\n9kcfC9Zbbn/M2tezt2U1/aTqvjrptWVtW5XrmNxn6tknk/vSp1068FST/aPavidvOoau/9CG/PDD\nDz/8lPlT5LiSdQVpiekPhB5/StJr8wQp7V8ef3l6X+SL+uhn7mixxgqvo7x6x3WXE3eVceapf5ry\naWXHzye11f6v5m3HIu0W/9r3+eRY8+2PrH09e1tW00+q7quTXlvWtlW5jsl9pp59Mrkv1RtL3T8k\nVbQhP/zwww8/7fopclxJHcqbwWOSZkKPXzxcto6Z3Rp6uOLuK2tL7Ngy/vsDCavbvnXqCBPXUWa9\no7qTmnTa+quMM0/905RPKzt6PrGtNudvxyLtFvfaSyO/89Sbtq48fT3altO8tmxlt3naa+t4L2Vd\nx6Q+k7fOaU3qS3XHUi0z2ytpb8NhAACACpSRmN4l6d2S7jSz10v6uidcX+rut06u6sz59NWdPTdt\ngNnWUbTeUd0XE56btv4q48xT/zTl08qOnk9sqwv527FIu8W99mLkd556s64rrm63yeutup9kUXab\nj16btu1ZVbmOSX0mb53TmtSX6o6lWsN/Zq6MHpvZLY0FAwAASpXldjG/K+n/lfRyM1s1s58ys5vM\n7CZJcvclSX9lZicl3S7pZ/OHs7oozZ8cP35f5PkDp6TTh/PXH7eOsuod1f2lx9fH/a4npq+/yjjz\n1D9N+bSyo+evj1nPgVPS6ofyt2ORdot77Z8/If3UE8mx5t0fWft61rYsK648ym7zrNteRnxF1zGp\nzxSNO6tJfel6SfMX6osFAAAgp7aNNw4m8pg7GlwH9cbPSW+9b3w9WJkT/4zWUV6947rf+LnxtVxv\nva/YBDnVxJmn/mnKp5Vd+/z+r0bLFWnHIu0W99q0WMtp/+S+Pl1blt9PmmjzsretynWk95nq90lK\nX7ql6f5R3XbLm46h6z+0IT/88MMPP2X+FDmu2LCCypmZuycNm4srH8QX/bvcmKqpd1S3JJVRf5Vx\n5ql/mvJpZUfPJ5Ur0o5F2i3utWmx5pW1r2dty7LiyqPsNi+j3rrWManP5K2zSAzRdbehf5Rt2uMK\n1qMNAQBlKnJcSR3KCwAAAABAlUhMAQAAAACNIjEFAAAAADSKxBQAAAAA0CgSUwAAAABAo0hMAQAA\nAACNIjEFAAAAADSKxBQAAAAA0CgSUwAAAABAo0hMAQAAAACNIjEFAAAAADSKxBQAAAAA0CgSUwAA\nAABAo0hMAQAAAACNIjEFAAAAADTq0qYDAAAAQDeYzQ6kmQVpxxbpzHlpddH9xFLTcQHoPhJTAAAA\npAqS0j2HpCO7x0vnd5nNiuQUQFEM5QUAAEAGMwtrk1IpeLzzYDPxAOgTElMAAABksGNL/PLtW+uN\nA0AfkZgCAAAggzPn45efPVdvHAD6iMQUAAAAGawuSvMn1y47cEo6fbiZeAD0CYkpAAAAUgUTHB2/\nWRosB0sGy9K9C0x8BKAM5u71rMjM3d2yl5e7y6J/lxtTNfWO6pakMuqvMs489U9TPq3s6PmkckXa\nsUi7xb02Lda8svb1rG1ZVlx5lN3mZdRb1zom9Zm8dRaJIbruNvSPsk17XMF6tCHy6uNnCoDiihxX\nOGMKAAAAAGgUiSkAAAAAoFEkpgAAAACARpGYAgAAAAAaRWIKAAAAAGgUiSkAAAAAoFEkpgAAAACA\nRpGYAgAAAAAaRWIKAAAAAGgUiSkAAAAAoFEkpgAAdIyZ/ZaZPWlmD0wos2hmD5vZ/Wb2mjrjAwBg\nWiSmAAB0z29L2p/0pJkNJO1292sk/Yykj9QVWBuZzQ7M5pbNblgJfs8Omo4JALDWpU0HAAAApuPu\nnzGzqycUebukO4ZlP2tmzzOzK9z9yTria5MgCd1zSDqye7x0fpfZrNxPLDUXGQAgjDOmAAD0z1WS\nVkOPH5X04oZiadjMwtqkVAoe7zzYTDwAgDicMQUAoJ8s8thjC5ndGnq44u4rVQXUjB1b4pdv31pv\nHADQP2a2V9LeMuoiMQUAoH8ekzQTevzi4bJ13P3WOgJqzpnz8cvPnqs3DgDon+E/M1dGj83slrx1\nMZQXAID+uUvSv5QkM3u9pK9vxOtLA6uL0vzJtcsOnJJOH24mHgBAHM6YAgDQMWb2u5LeLOkFZrYq\n6RZJmyXJ3W939yUzG5jZSUnflPSu5qJtlvuJJbNZSYOD0tJ+abAsnT7MxEcA0C7mHnvJSfkrMnN3\nj17vMqG83D24Pib8d7kxVVPvqG5JKqP+KuPMU/805dPKjp5PKlekHYu0W9xr02LNK2tfz9qWZcWV\nR9ltXka9da1jUp/JW2eRGKLrbkP/KNu0xxWst9HasI/vg6bQlgDiFDmuMJQXAAAAANAoElMAAAAA\nQKNITAEAAAAAjSIxBQAAAAA0isQUAAAAANAoElMAAAAAQKNITAEAAAAAjSIxBQAAAAA0isQUAAAA\nANAoElMAAAAAQKMubToAAADQbmazA2lmQdqxRTpzXlpddD+x1LV1AADai8QUAAAkChLGPYekI7vH\nS+d3mc2qrMSxjnUAANqNobwAAGCCmYW1CaMUPN55sFvrANKZzQ7M5pbNblgJfs8Omo4J2Cg4YwoA\nACbYsSV++fat3VoHMBln7oFmccYUAABMcOZ8/PKz57q1DiANZ+6BJpGYAgCACVYXpfmTa5cdOCWd\nPtytdQBpOHMPNInEFAAAJAqGMB6/WRosB0sGy9K9C2UObaxjHUA6ztwDTTJ3r2dFZu7ulr283F0W\n/bvcmKqpd1S3JJVRf5Vx5ql/mvJpZUfPJ5Ur0o5F2i3utWmx5pW1r2dty7LiyqPsNi+j3rrWManP\n5K2zSAzRdbehf5Rt2uMK1itybK4upmqPzX17HzSlj20Zf43pgVP8kwTIrsixmTOmAAAA2PA4cw80\nK/WMqZntl/RrkjZJ+g13/2Dk+RdI+s+SrlQwy++vuPv/HVMPZ0wL1MUZU86YcsaUM6bTxMAZU2TB\nGVPk1fe27Pv2AVWp7IypmW2S9CFJ+yW9UtKNZvaKSLF3S/q8u3+PpL2S/oOZcRsaAAAAAEAmaUN5\nr5V00t0fcfcLku6U9I5ImSckPWf493Mk/a27Xyw3TAAAAABAX6Wd2bxK0mro8aOSXhcpc0TSH5rZ\n45J2SLqhvPAAAAAAAH2XdsY0y5S9vyjpC+7+IknfI+nDZrajcGQAAAAAgA0h7YzpY5JmQo9nFJw1\nDXuDpNskyd1PmdmXJb1c0n3Ryszs1tDDFXdfmTJeAMAGZWZ7FcxlAAAAembirLzDSYz+QtIPSnpc\n0p9IutHdHwyV+VVJf+/uv2xmV0j6U0n/2N2/FqmLWXkL1MWsvMzKy6y8zMo7TQzMyossmJUXefW9\nLfu+fUBVihybJ54xdfeLZvZuSccU3C7mN939QTO7afj87ZL+T0m/bWb3Kxga/J5oUgoAAAAAQJLU\n+5iWtiLOmBaqizOmnDHljClnTKeJgTOmyIIzpsir723Z9+0DqlLZfUwBAAAAAKgaiSkAAAAAoFFp\ns/ICAAAA6Biz2YE0syDt2CKdOS+tLrqfWGo6LiAJiSkAAADQI0FSuueQdGT3eOn8LrNZkZyirRjK\nCwAAAPTKzMLapFQKHu882Ew8QDrOmAIAAGwADO3cSHZsiV++fWu9cQDZkZgCAAD0HEM7N5oz5+OX\nnz1XbxxAdgzlBQAA6D2Gdm4sq4vS/Mm1yw6ckk4fbiYeIB2JKQAAQO8xtHMjCc6CH79ZGiwHSwbL\n0r0LnB1HmzGUFwAAoPcY2rnRDJPQJTO5+9Jc0/EAaThjCgAA0HsM7QTQbiSmAAAAPcfQTgBtZ+5e\nz4rM3N0te3m5uyz6d7kxVVPvqG5JKqP+KuPMU/805dPKjp5PKlekHYu0W9xr02LNK2tfz9qWZcWV\nR9ltXka9da1jUp/JW2eRGKLrbkP/KNu0xxWsV+TYXF1M1R6b+/Y+yKOMduh7W/Zp+/q0LWi/Isdm\nzpgCAAAAABrF5EcAAABoteA+rDMLwezCZ85Lq4sMQwb6hcQUAAAArRUkpXsOrb0P6/wus1mRnAL9\nwVBeAAAAtNjMwtqkVAoe7zzYTDwAqkBiCgAAgBbbsSV++fat9cYBoEokpgAAAGixM+fjl589V28c\nAKpEYgoAAIAWW12U5k+uXXbglHT6cDPxAKgCiSkAAABaK5jg6PjN0mA5WDJYlu5dYOIjoF/M3etZ\nUYGbeFd1Y+Cqb+ItSWXUX/WNkaetf5ryaWVHzyeVK9KORdot7rVpseaVta9nbcuy4sqj7DYvo966\n1jGpz+Sts0gM0XW3oX+UrchNvBEocmyuLqZqj819ex/kUUY7NNWWda23T32lT9uC9itybOaMKQAA\nAACgUSSmAAAAAIBGkZgCAAAAABpFYgoAAAAAaBSJKQAAAACgUZc2HQAAAACAdjKbHUgzC9KOLdKZ\n89LqIrfqQRVITAEAAACsEySlew5JR3aPl87vMpsVySnKxlBeAAAAADFmFtYmpVLweOfBZuJBn5GY\nAgAAAIixY0v88u1b640DGwGJKQAAAIAYZ87HLz97rt44sBGQmAIAAACIsboozZ9cu+zAKen04Wbi\nQZ+RmAIAAABYJ5jg6PjN0mA5WDJYlu5dYOIjVMHcvZ4Vmbm7W/bycndZ9O9yY6qm3lHdklRG/VXG\nmaf+acqnlR09n1SuSDsWabe416bFmlfWvp61LcuKK4+y27yMeutax6Q+k7fOIjFE192G/lG2aY8r\nWK/Isbm6mKo9NvftfZBHGe3QVFvWtd4+9ZU6jqHASJFjM2dMAQAAAACNIjEFAAAAADSKxBQAAAAA\n0CgSUwAAAABAo0hMAQAAAACNurTpAAAAAABgxGx2IM0sSDu2SGfOS6uL3KKm/0hMAQDoGDPbL+nX\nJG2S9Bvu/sHI83sl/TdJfzVc9HF3/0CtQQJADkFSuueQdGT3eOn8LrNZkZz2G4kpAAAdYmabJH1I\n0lskPSbpc2Z2l7s/GCn6aXd/e/nr50wGgCrNLKxNSqXg8eCgJD5reozEFACAbrlW0kl3f0SSzOxO\nSe+QFE1Mc93gPB1nMgBUaceW+OXbt9YbB+rG5EcAAHTLVZJWQ48fHS4Lc0lvMLP7zWzJzF5Z3urj\nzmTsPFhe/QA2tjPn45efPVdvHKgbZ0wBAOgWz1DmzyTNuPu3zGxO0u9JellcQTO7NfRwxd1Xpg+J\nMxkAyrK6KM3vWvtPsAOnpNOHm4sJSYZzGuwtoy4SUwAAuuUxSTOhxzMKzpo+w93PhP4+ama/bmaX\nu/vXopW5+63FQ+JMBoByuJ9YMptVcE3p0n5psCydPszlAu00/Gfmyuixmd2Sty4SUwAAuuU+SdeY\n2dWSHpf0Tkk3hguY2RWS/trd3cyulWRxSWk+8yc5kwGgSsMkdMlM7r4013Q8qAfXmAIA0CHuflHS\nuyUdk/QlSR919wfN7CYzu2lY7MclPWBmX1BwW5mfLC+C4zcHZzCk4Pe9C5zJAAAUZe5ZLlUpYUVm\n7u6ZZwgM/kMSzCgY/rvcmKqpd1S3JJVRf5Vx5ql/mvJpZUfPJ5Ur0o5F2i3utWmx5pW1r2dty7Li\nyqPsNi+j3rrWManP5K2zSAzRdbehf5Rt2uMK1tuIx+a+vQ/yKKMdmmrLutbbp75SxzG0Kn3aDxtF\nkWMzQ3kBAEDncX9VAOg2ElMAANBpQVLK/VWrtjb5/5jMZge0L4CycI0pAADouJkF7q9arXHyf3Sf\n9LE3B0v3HAqWA0BxJKYAAKDjdmyJX879VctD8g+gWgzlBQAAHXfmfPxy7q9anuaS/64OIea6Z2A6\nJKYAAKDjVhel+V3cX7VKzST/8dcP7znU9uuHue4ZmB5DeQEAQKcFX/S5v2q1Vhel+ZNrl9WR/Hd1\nCHFX4waawxlTAADQecMkdCm47+HSXNPxdFnSEFSzWUmDg9LS/qBkHcl/V68f7mrcQHNITAEAACAp\n0xDUJTO5pJqGpHb1+uGuxs21sWgOQ3kBAAAw1LYhqE0NIS6qm3Gvvy3Q0X3cFgh1ITEFAADAULuG\noK6/fljqwvXD3b3uuW3/mMjHbHZgNrdsdsNK8JvEugsYygsAAICh9g1BbWYIcXHdvO65Xf+YyIMZ\nkbsr9Yypme03s4fM7GEze29Cmb1m9nkzO2FmK6VHCQAAgBp0cwgqytK+f0xMrx9nfTeiiYmpmW2S\n9CFJ+yW9UtKNZvaKSJnnSfqwpB9291lJP15RrAAAAKhQd4egohx9+MdE9KzvPZLeL+nZexjW225p\nQ3mvlXTS3R+RJDO7U9I7JD0YKvPPJH3c3R+VJHf/mwriBAAAQA26OQS1nYIk6ITMbljpwgy3628L\nNFiWTh8OljcdXVbRs77HJN0mSc+VtI9hve2VNpT3KkmrocePDpeFXSPpcjP7IzO7z8z+5zIDBAAA\nALpmfK2j1KUZbt1PLI3+IeG+NNe9BC561ve2yPMM622rtDOmnqGOzZJeK+kHJW2TdNzM7nX3h4sG\nBwAAAHRT0rWOg4OSOpbsdcf6s75xujOZUxtVda/btMT0MUkzocczCs6ahq1K+ht3PyfpnJndI+m7\nJa1LTM3s1tDDFXdfmTZgAMDGZGZ7Je1tOAygUlV94UMTuj/DbVdFZ3Jer/rJnPr6Xq5y1uO0xPQ+\nSdeY2dWSHpf0Tkk3Rsr8N0kfGk6UdJmk10n61bjK3P3WArECADaw4T8zV0aPzeyWxoIBKsBtLqbT\n/i/+fZjhtg/mT659T1U/mVO/38vVjQSYmJi6+0Uze7eCq4Y3SfpNd3/QzG4aPn+7uz9kZsuSvijp\naUlH3P1LRYICAADYeKr5wrc+gTtaMM7m1fHFv3jiu7ooze+qOylC1PGb4yZzqnadfR7GXd1IgLQz\npnL3o4p8grn77ZHHvyLpV4oGAwAAsHGV/4UvPoELlnf7zE0dX/yLJb7jax31yWBJXUkRwpqZZbrP\nw7irGwmQNisvAAAAalHFF764BE7q/qykdXzxj0t8p2u3cBLazRlukU+fh3FXd69bElMAAIBWqOIL\nX1/P3DT1xb/r7Ya8zGYHZnPLZjesBL8n3fanuuStacE/V47fHIwAkILf9y7UMSsvAABAK7R/spti\n1t/mooyhn309c9PU9ZtdbzfkMe01zdW8l9ujquHR5p7lVqUlrMjM3d2yl5e7y6J/lxtTNfWO6pak\nMuqvMs489U9TPq3s6PmkckXasUi7xb02Lda8svb1rG1ZVlx5lN3mZdRb1zom9Zm8dRaJIbruNvSP\nsk17XMF6XTo2J3wxPCkdvzn8Za+Oz5I6lPe5F3+NqTT7trQvyVUcm4uIrjfYtp2VfPEP1hU3m+v0\nZ4aaaq/w+sv4HtXE+6fMY3Ox7yhzy9LRfeufGSynJWZd+tyZVvx3n/zHZobyAgCADkia7Kbr10pW\nK37Y3dprH7vK/cTSKCmo5vrNaoYroov6OiS+XRjKCwAAOoAvhnlFh92NzuBhsmZmc0U79XVIfLtw\nxhQAAHQAXwwBNKW/kxm1CYkpAADoAL4YAmhGlTPRYozJj5j8iMmPMr6WyY/SMfkRkx/VicmPiuva\nsTnLZDfBcWP2bXGz93bpfVBFrNMex9o++VF4eVVtVUb9TH6UX1smPyoSU5c+d6ZV9uRHXGMKAAA6\nIfs1f/G3dZBOVB0iACAnElMAANAzcbP3Djo3e29wi4p+3rMV8fp+r15gEhJTAACwAXRn9t4gOTmh\ntfdNDM76kqT0V8K9etnvqEUb/ilCYgoAADaALs3eO7OwftkzZ31JUHor6V697HdUqy3/FGFWXgAA\n0DNdn72Xe7ZuTOx3NCXpnyI7a70EgsQUAAD0TNdv68A9Wzcm9vtGYzY7MJtbNrthJfg9O2gmknb8\nU4ShvAAAoFeyz97bVquLkvatXda1s76Y3uqiNL9r7Zkr9ntftWX4bKAd/xThjCkAAECLjL+Udvms\nL6YV7N+un+1Hdu0YPhtYXWzDJRAkpgAAAC00PtvrJr3yPc0O9UMd3E8sjfa7+9IcSWk/Be/jbdfG\nP1v/NcVt+acIQ3kBAEAt2nA7gi7htjHtRV+uV5/aezyE94rnx5do5priNlwCQWIKAAAq167rqbqC\n28a0EX25Xv1r79EQ3nskvU/SbaHnNvY1xQzlBQAANWjT9VRd0Y6ZMhFFX65X39p79L5+k4Kk9JeG\ny3/oaxv9mmISUwAAUAOSrOm1Y6ZMRNGX69W39o6+r/+P4e+n/2QjJ6USiSkAAKgFSdb0VhfXL9vY\nQ/3agb5cr761d9wMuBLvaxJTAABQi3bcjkBq003tJ+O2MW3Vnr68MfSrveNnwFVHr5ctl7l7PSsy\nc3e37OXl7rLo3+XGVE29o7olqYz6q4wzT/3TlE8rO3o+qVyRdizSbnGvTYs1r6x9PWtblhVXHmW3\neRn11rWOSX0mb51FYoiuuw39o2zTHlewXt3H5iAB3HlQWtoffBk7fTj6ZSxLvVnfy3HvA2n2bTET\nqZyUjt/cpi+G0WNOGe/haeuq4thcRNJ6q/h8S+vrWfpyWtxZ119U3rracBwZ99ns7R33+ujfZcRU\nRvmqvlsWUfR7fpFjM4kpiSmJacbXkpimIzElMa0TiWlxTR2bi74XiyWmc8fW3n5lZLDc1C0S4pCY\nxscTt96FKADEAAAgAElEQVQmEtNp1k1iml/RGEhMp9dkYspQXgAAsIH0bSIVAOgH7mMKAAA2kL5N\npIK+MptbDv6Rcua8tLrYpqHmaE4wrHlmoY99g8QUAABsIKuL0vyutdeYdnciFfRPkHic0Noh5/O7\nzGYLT5DT56RmWl1siyDmddfIl9I32oDEFAAAbBjuJ5bMZiUNpp5IBajHzML6ZUd2B31WuftpWUlN\nFxO6qO4meDMLa2OWyugbbUFiCgAANpThF8+lYOKO6SY86sOXcrRd9DroeyTdLenZe4LhvXn7XPGk\nprsJXVS0Le6RdMVu6eX/yWzuc+19X/f7GnkSUwAAgAzK/lJOkot40eugj0m6TZKeK2lf/j5XRlLT\nlzN20bZ4po0vV6E2rlq/r5FnVl4AANB7QRIomd2wYja3PHo8naQv5TsP5otnz6HgOsKPvTn4vedQ\nvrjQL6uLax/fFnk+X58rJ6npyxm7aFuU1cZVW10M7rsc1p9r5ElMAQBAr42TQCmcBE5fU5lfystL\nctEv47N0g+XkUnn6XBlJTV/O2MW1RVT7ku2gbxy/edw3BsvSvQvtO7ObD0N5AQBAzyUlgWPZhtWW\n+aW8L2eeUBX3pTkzefyz0/e5cib+6ses1uvbIk47k+0i18i3HYkpAADouaQkMJD92tEyv5RHk9zR\nBDffeXUwwc3RTLUUuU6Va1y7Yv5kWYlg0aSmiVmt1/bTj8lsdlDG+sJtUWYbowB3r+UnWNU05d3j\n/i43pmrqHdVdVv1Vxpmn/mnKp5UdPZ9Urkg7Fmm3uNemxVrGuibVnbUtm/wpu83L3rYq1zGpzzTR\n/tF1t6F/lL+96t02tb0Nyzo2F30vZn39+HN7//LoeBL9mfz83NH19b9qIM0dHT0vvWqQrw1eNZAO\nPDxe1y/GrP9Vg0nHnvV1uAeP42MK15XltdMem4M69y9LP7ES/M7XNpP2bV2ftWUdmyfFnWX9432V\nv8/lbbO040gdx5Vxv/RMfTzPfpy2jcv83lzku2Ub8qP4vqXccdmwgsqZmbu7ZS8fxBf9u9yYqql3\nVLcklVF/lXHmqX+a8mllR88nlSvSjkXaLe61abHmlbWvZ23LsuLKo+w2L6PeutYxqc/krbNIDNF1\nt6F/lG3a4wrWa+rYXPS9mPX148/tuDOiB05Jv7EreP6GleDa06gbPu3+sb1540wTxLVzwlDCwbK0\ntD/p2DM8s7ov7nVxZ8PCxzFp7ljaa6c5NiecdT4pHb+5rLNpSd8Jqvh8K+vYPCojTfddJu47R97t\nzHt8SjuO1HFcSR7KHN/H416f1n7TtnGZ35uLfLdsQ34U37fyH5uZ/AgAAPTaeMKQkWDCkPHjZiZ0\ncT+xNPnLddr1pkWuUy37Glcmc0KduBa7j0hMAQBA74XP2rkvzcVcO9rCWzCkJcZFEuqyk3Emc0Kd\n2jkxEYohMQUAABtaO27BEHfrirTEuEhCXXYyXu9ZZ7PZQXA/2htWRo+rWA/aoI3/NEIVSEwBAMBU\nRkmA2Q0rQXLQ/aQgPKx2/RnVOkQT47VneeMUSajLT8brO+s8vp716L7xtcF7DvWhHyJOuJ9Kfbpv\nJ9Zi8iMmP2Lyo4yvZfKjdEx+xORHdWLyo+LyHJul2bcVneSm7smP1sY/PqYUed+W+Z6Kn7Ap/neR\neJLqyjIpTHS5FDcJUXgyp/JvIzKeCCd94qYy1sXkR+2Y/Ci8/VK+doz+PW2ZpPLTxpD0HJMfBbiP\nKQAAmELSJDeDg5I4i7GBFb1HZnZczwq0RTA7+PheyEXqIjEFAABTIClA05qZRRnAWDBC4oTWjl6Y\n3yX9ee46SUwBAMAUSArQtNXF4Atw9L60013PGnyxnlkIn+2ZNPQ4rnzwxRzYiGYW1i87slv6jdw1\nkpgCAIAplJMUAHm5n1gym1UwfHxpf7B0uglxxhMorblWepfZbOykU0nlc28E0HlJo2fyIzEFAACZ\njZMCfTJYUv4kN0Ca8PWsocdTmPZa6aTywEaVNHomP24XAwBAx5jZfjN7yMweNrP3JpRZHD5/v5m9\npsz1h5OAZm6tAhQ17bXS5ZwdCt9/tS+3WsJGFTfR0YFTRWrkjCkAAB1iZpskfUjSWyQ9JulzZnaX\nuz8YKjOQtNvdrzGz10n6iKTXNxIw0ErTXitd/OzQtMOHgTYLRs9IwaiZ8S2i9MxomumRmAIA0C3X\nSjrp7o9IkpndKekdkh4MlXm7pDskyd0/a2bPM7Mr3P3JaGVm152RtEWK3rdum6QLT0nnn5Z2PC1d\nNOmYzN72DWnLRenjKa+PM6rzbpld901Jzxq/Lvzc274RrO+cjddtT0vb/j/pb56S7gmtd5tJl52T\nzv6l9Phd0gmZ3bAifWX7uNy2y6TLvi19YhjzD31NurA52J5RHNssObboNlx2TvqEzN74sHTplUH9\nFy4GbbXJkl87qa7r/1TS5mD9b/my9KnQ7x/7YlBm9hbp8h+WNr183OZxbTralksuSBe/Jm3eLl04\nu7bO674pbdscxB3d5qT9NI43aMfrvyhtemFQ9+bt0tYn1u+fafpFUvslrXfrE+Nlc8vSd140jufi\n8yXbvLY9wvV/26R3Kfh5k6S3STonaev1Ztd/a20/PCbp79+8tvxPSvr609LyJUGbjtY3atPR+2b8\nXpGuuks6skm6R9LvSHpC0rd3S//Tfze77um1bRHuB+M41q5rUj8btdlXvyJ9JvJ+mbRvv/oVaes3\nxn1x9H7/678d1zPq72nvlSKinwfhz5zrvyg9tTO+DUavG/X9aPvFfbYci5RJ+my65ELwGRR+70/q\n41k/78Lv97h9G943kvTcf7D+82tS+2d5f0Xb7C1fHn9uXHy+9OxLxu22ycKfHdJT3xXUYz8gvfx1\nzMoLAMDGcZWk1dDjRyW9LkOZF0tal5hKL92+ftmVkvZJumNr6O/hc5/cEXyxnvT6OOE6Jeml25Kf\ne++OYH2j5cck3SbpnleP43jp9uD52yRpu3TPa6X/9N3Bc+9+87jcG7YPX7t5vL73PD9Uz7bJsUW3\n4TZJ9wy3+WW7Q6/dHP+aJNG63v/aYDsl6XVXr/39iVcHv1/+ful5oe9ucXE/02ZbpX1bpWPPGT5+\nQVDmJcM637BtbdyjbU5qi2i8kvT9rx7XfZuke16wdv9kbYfw+uKej1vvbZL0gvGyX9gXtN8+SXe8\nYO3rr5K0slXaofXLP6Qg0bxyFMsl6/uhJF23ZW35n5V0bHhJ3KhNw33hmbhD67ximJTeIek3Jf30\n6AmTtGntvgvX+94d474xWjZJuM3uGJ6dHb1f0vbtsd3x7/c7nh88fsP2+H5Tpmic0c+cmVdPft2o\n79/xnGB5eD9GP1uiy6X4z6ZjW6XbtgafQVLw3p9mG5LqlMbv87h9G92X0bjT2j/L+yv6eTGKafRe\niq7zmddsXhv3J4dvsPz/oyAxBQCgWzxjuei3g4TXvTBm2QckvX/4XPjvkbsVfDlPen2caD0vnPDc\n3ZF13xZZrtDz9wyXn5R056b15cKvfdOEepJii25DuK6k9skiWtfdoVijv+8ZlnlF5HtbXNzhNov+\nlta2S9w2J7VFNN5wHUn7J4u09ktab9So/aJ1XS/pw5JelrDe0fK4fRntR+Hy4f0V13Zxce/U5DaK\n9vdwHNFlk4TXnaWfh8vHtWFcO0y7n6eR9zMn2vfj9mP0syWuTNxnU/RzJG27s3zeRd/nSZ/F9yjo\nw9ek1JsljqTno+/lcBtG60jqp8WZe9bjW8EVmbm7Z06hzeTuwUE1/He5MVVT76huSSqj/irjzFP/\nNOXTyo6eTypXpB2LtFvca9NizStrX8/almXFlUfZbV5GvXWtY1KfyVtnkRii625D/yjbtMeVvjCz\n10u61d33Dx//gqSn3f2DoTL/UdKKu985fPyQpDdHh/Kamd8Serx3+AMAQBYrw5+RX5aU99jMGVMA\nALrlPknXmNnVkh6X9E5JN0bK3CXp3ZLuHCayX4+7vlSSflnvC/09MvqPeNzfCj2eRrieSc/F/R1e\nb7SeD2i9cLlozHFxTIotWi5aV1rdWeoKb0M09rjfk+KObndY2mvTlsfFF61/mjaIri9peyatN2lf\nj5ZdKuniFOtKej6pX07qC+HHo2U/K+lyxUuqN2lZkqT3adZ9m/R3nv2bR1K7Zt1/WV6TtT2n+RxJ\n2oZJdWbZR0kxZDEaMRA+2zoplixtOGm5VGQoL7eLAQCgQ9z9ooKk85ikL0n6qLs/aGY3mdlNwzJL\nkv7KzE5Kul3Bt+EET8T8vE/BF5onhmWuDy0PP57m530TXve+yDqif4fXG475dMI2hctFY46LY1Js\n0XLhuhTTVtO0x6iu90XqCv++NFQuWv+kbQnXGY1z0jZPWh6ON24dcTFO0y/iXjtpvaPl8xfWLg/X\n9WCo/p/O2HZxzyf1y0l9QTF1/GSGtsiyrmn66TT7diRu+TTvlSI/0e0fxZClr0T7QdpnS9o2xX0G\nZd0HWeq8XsHnWFK9o8+48L6Zpv0/LOlfZ2yzuPfY6HVx/TxuPxXDUF6G8jKUN+NrGcqbjqG8DOWt\n00YdylsmM3Np7xnJtkrRtnxm1sit0tvOSBcknb9kONPohWCWTimYmTLu9XGis0Pas8avi84CHF7f\nBUl375B+7IFgIqDRjJjP2i4di6z3lyR96WvSxy8fl3v2FulZ35a+Ovxm9dx/ID21Oah/FMdWS44t\nug2XnQvq+sw10nVngvqfuhi01XXfSn5tWl3XP7B2ltsLZ6XNL5KOPisoH57N9ZxL/nR8m4625e6t\n0lsekT51dfA7PMumbZaevTmIO7rNSfspHO/Wb0h6lnT3q8friO6faftFUvslrffHHpDOPyadvlfa\n+Xrp4lVr47nuW9JTl0kv2xTMpvtM22nYdp7cdtF+uMnWttfouUueli7+3Xh9ozKj9030vTLqK0n9\n7JkZUp8K6h3NiBq3rkn9LNxmo3Xb1uR+/szMx8NJpq5/IJj59tmXBNsQfT9N2oYyRD8PRu0Y3v9J\nfSXc98PtF/48mfR30mfT3VuDPvc335bu+b7gvT+pj2f5vLvkqWB7lrdGZmuW9LRLmzX+jBtda36j\npL+9KJ3/dnr7P3eT9Huh174tVLdHZoIe9blnZvG+evxeGh0D7t6xtp9H++mWC9InLs97bCYxJTEl\nMc34WhLTdCSmJKZ1IjEtLkhM9x+TVheT7qPYlr4T9xkcf1/IA6ekexekE5+sOu4y32fJx8HkbUy7\n92VVx61J66jqOJC1bHI/2XlQ2r41uE/p6cNl3ze0zs/cJr6/tuVzYKSJ70RV5SaTP8ekvO//oO65\nZenovvXPDJbdl+YmvzbuvZS+3UWOzVxjCgDAhnZ0nzS/y2xWZX9Zr1pwg/dZSYNh0vGxN4++sFlr\nvkJPFnwpnVmQjg6/RK79J8H6bawmseqzYVvRXmiltPd40mdcttpXF6X5XesT29OHy9+S4jhjyhlT\nzphmfC1nTNNxxpQzpnXijGlxwRnT0feA+P+gt6XvZPkMruO7Q9L68qwz/kzJ/Enp+M1lJJ4b6Yxp\nKMHfJ80dk47ua+rzN+5xlesqu26JM6ZpdTbVHnnWm3fEAGdMAQBAQ7Yn3IAd1ZlZWJuUSsHjwUFx\nhi+z9Ql+MHTRbHbAmWVsdF0aMcCsvAAAQMF/0tvHbHYQDHEdXS/VJzu2xC/nnwTTiUvwpeAsEbII\nv89Gj5uMBxtTamJqZvvN7CEze9jM3juh3Peb2UUz+7FyQwQAANVq5zVH4zNho8k7xmfCmowrKpo8\nZ4/vzPn45e38J0F7keAXsf59Jkl7DrXtfYb+m5iYmtkmSR+StF/SKyXdaGavSCj3QUnLUnvGnwMA\ngDSD5ekm06hT+8+ExSfPWb/Ury4G15SGtfOfBO1Ggl9M0pDy9rzPsDGkXWN6raST7v6IJJnZnZLe\noeBuxWEHJf1XSd9fdoAAAKA6abcMaFYXzoTlv06UGXfLEjfzqESCn1UX3mfYCNIS06skrYYePyrp\ndeECZnaVgmT1nypITOuZ5hcAAPRcF86EFftS36WJSdIM75kYe9ubKsUn+Ev7SfCzSn6fpd3OCChT\nWmKaJcn8NUn/xt3dzEwM5QUAAKXowpmwLiTP1QqSlxNaO5y53nvjRhP80W1PkEXivS7vXT/bcTfv\neYxuSEtMH5M0E3o8o+Csadj3SrozyEn1AklzZnbB3e+KVmZmt4Yerrj7yrQBAwA2JjPbK2lvw2Gg\nRt04E9atG9hXY2Zh/TJue9MV699nH3uzdO8CtzNC3dIS0/skXWNmV0t6XNI7Jd0YLuDuLx39bWa/\nLen345LSYdlbC8QKANjAhv/MXBk9NrNbGgsGtWn7mTCuE5W4RrH7wu8zM3nQr294T3xp9iuqMTEx\ndfeLZvZuScckbZL0m+7+oJndNHz+9hpiBAAAaK0+XSeaD8OZ+2D99aTf2RFfkv1at41yrW/qfUzd\n/ai7v9zdd7v7vx0uuz0uKXX3d7n7J6oIFAAAAG20urh+2UYbztxt8bc92n6l9FNPrC3Z7v2a/57C\n7VXsllTdYu71jIgxM3f3zBMjBcMIgomUwn+XG1M19Y7qlqQy6q8yzjz1T1M+rezo+aRyRdqxSLvF\nvTYt1ryy9vWsbVlWXHmU3eZl1FvXOib1mbx1Fokhuu429I+yTXtcwXpZ2rBtfSfr+7iOuNvWNlFV\nHbeS1hHcEzf/cOYqvl9Ms/4i6vzMraLu4YzK+9Y/86b7pO1/07Zh6vHH3FECF74udv6kdPzmMmKu\n+/NlvN6kfTNYrvJ2X3GfH1m2u8ixOe0aUwAAgM7YKEPe2qjd98TFZEnXCV/5TfePdWS/9nWypo1z\nDXfqUF4AAIDu2BhD3oBy9eE64WoSuOaHB/dh32RDYgoAAHok7ozJzoPNxAJ0xepiMOw1rN3Xk65X\nfgKXdH1n3vryqXffRBPxKtaRhKG8AACg5/o35A0oUz9ue1TFPYWThgfXp859s/463SAhN5sd1NEX\nSEwBAEDP9W/IG1C2rt/2qJoELml4cL3q2zdxibg0HHVCYgoAAJDd/Mlyz5gA6IryE7ik4cF91exE\nS1xjCgAAeuT4zcFtS274dPD73oVuDUcE0B5J13f2VbMTLXHGFAAA9EbXhyMCaI+k4cGSPtl0bNWI\nu05XqmvUibl7HeuZ+marddzEtuqbH0tSGfVXfRPfaeuv4gbYSeWKtGORdou/cXM1NyrP2tfbdjPx\nsmMosu1tWMekPpO3ziIx1Hmz96YUuYk3AlnasG19p47PijJiaYOqjltVraOK7xdF4plGnZ+5be93\ndWi6DZpef5WCCZB2hhLxpf3T5GRFjs2cMQUAAAAArBt1MjpJVAeuMQUAAAAANIrEFAAAAADQKBJT\nAAAAAECjSEwBAADQGLPZgdnccvD33HIw+QqAjYbEFAAAAI0IktA9h6Sj+4IlR/dJew6RnAIbD4kp\nAAAAGjKzsP6eiUd2B7erALCRkJgCAAD0SHRobPXrKzIMd8eW+OXbtxaNC0C3kJgCAAD0RPzQ2NHy\nKtY1Xke+YbhnzscvP3uuWHQAuobEFAAAoDeiQ2PfP1p+R/nJ6czC+mXTDsNdXZTmT65dduCUdPpw\nsdgAdM2lTQcAAACAsoyGxt4j6U2SPjBcfvQF0vwhs1m5n1gqd11R2Yfhup9YMpuVNDgYvO7sOen0\n4fJiBNAVJKYAAAC9MRoae7eCxDTsyO4gAVRJSV85w3CHSSiJKLDBMZQXAAC0Hve6zGo0NDbp3EOZ\nkwqtLq5ftrGH4dJP0Ud19WvOmAIAgFYbT+gzunby6D5pfle5w1L7YTw0duYOSS9YX6K8SYWCdUnS\nYJlhuPRT9Fk9/drcvcz6kldk5u5u2cvL3WXRv8uNqZp6R3VLUhn1VxlnnvqnKZ9WdvR8Urki7Vik\n3eJemxZrXln7eta2LCuuPMpu8zLqrWsdk/pM3jqLxBBddxv6R9mmPa5gvSxt2Ia+E/ynfjTza9hg\n2X1pblyuufdc26xPkqTgbOa9C2V+mWyqDdp2TAy+r8wdy9JPy1hXW/tdXZpug6bXX6fRd/H14vt1\nkWMzZ0wBAEDLca/LaTGpUBPop30X/MNnZkE6OvyH2erixn1Pld+vSUwBAEDLca/LPJhUqG700z5j\nqHZU+f2ayY8AAEDLca9LdAH9tN+i9wiWpr9vb1fV0685YwoAAFqNYanoAvpp323kodrHb66jX5OY\nAgCA1mNYKrqAftpnG3eodl39mqG8AAAAADARQ7WrxhlTAAAAAJiAodrVIzEFAAAAgBQM1a4WQ3kB\nAAAAAI0iMQUAAAAANIrEFAAAAADQKBJTAAAAAECjSEwBAAAAAI0iMQUAAAAANIrEFAAAAADQKBJT\nAAAAoOXMZgdmc8vB33PLZrODpmMCykRiCgAAALRYkITuOSQd3RcsObpP2nOI5BR9QmIKAAAwBc5c\noX4zC9KR3WuXHdkt7TzYTDxA+UhMAQAAMuLMFZqxY0v88u1b640DqA6JKQAAQGacuUITzpyPX372\nXL1xYKNoYmQIiSkAAEBmnLlCE1YXpfmTa5cdOCWdPtxMPOizpkaGXFpl5QAAAP3CmSvUz/3Ektms\npMHB4J8gZ89Jpw+7n1hqOjb0UdLIkMFBSZX1ORJTAADQacF/8WcWpKMKhp6tLlb3hX11UZrftfZL\nG2euUL1hnyYRRQ2aGRlCYgoAADprPORslCge3SfN7zKbVRXJKWeuAPRfMyNDSEwBAECH1T/kjDNX\nAPqtmZEhJKYAAKDDmIwIAMrU1MgQElMAANBhTEYEAGVrYmQIt4sBAAAdxm00AKAPOGMKAAA6i8mI\nAKAfSEwBAECnMRkRAHQfiSkAAB1iZpdL+qikl0h6RNIN7v71mHKPSPqGpO9IuuDu19YYJgAAU+Ea\nUwAAuuXfSPoDd3+ZpP8xfBzHJe1199eQlAIA2o7EFACAbnm7pDuGf98h6UcmlLXqw8FGZDY7MJtb\nDv6eWzabHTQdE4BuIzEFAKBbrnD3J4d/PynpioRyLulTZnafmc3XExo2giAJ3XNIOrovWHJ0n7Tn\nEMkpgCK4xhQAgJYxsz+QdGXMU+8LP3B3NzNPqOYH3P0JM/uHkv7AzB5y98+UHSs2opkF6cjutcuO\n7A5mRmYSKgD5kJgCANAy7v7WpOfM7Ekzu9Ldv2JmL5T01wl1PDH8/VUz+38kXStpXWJqZreGHq64\n+0qR2LER7NgSv3z71nrjANA0M9sraW8ZdZGYAgDQLXdJ+l8kfXD4+/eiBcxsm6RN7n7GzJ4t6XpJ\nvxxXmbvfWl2o6Kcz5+OXnz1XbxwAmjb8Z+bK6LGZ3ZK3rkzXmJrZfjN7yMweNrP3xjz/z83sfjP7\nopn9sZn947wBAQCAif6dpLea2V9K+qfDxzKzF5nZJ4dlrpT0GTP7gqTPSvrv7n53I9Gih1YXpfmT\na5cdOCWdPtxMPAD6wNyTLk0ZFjDbJOkvJL1F0mOSPifpRnd/MFRmj6Qvufvfm9l+Sbe6++sj9bi7\nZ54d0EzuHswmGP67TFXVO6pbksqov8o489Q/Tfm0sqPnk8oVacci7Rb32rRY88ra17O2ZVlx5VF2\nm5dRb13rmNRn8tZZJIboutvQP8o27XEF62Vpwz72HRQXTHS082AwfPfsOen0YfcTtV1f2rZjIu+T\netHe7VXk2JxlKO+1kk66+yPDld0p6R2SnklM3f14qPxnJb04TzAAAKAdgsRjZkE6quC2IKuLdSYe\naLdhX6A/AChNlqG8V0laDT1+dLgsyU+LDyoAADqL24EAAOqWJTGdPNY3xMyuk/RTktZdhwoAALoi\n6XYgOw82Ew8AoO+yDOV9TNJM6PGMgrOmawwnPDoiab+7/11cRUxJDwDIq8wp6ZGG24EAaB8uMei3\nLInpfZKuMbOrJT0u6Z2SbgwXMLOdkj4h6V+4+8loBSNMSQ8AyKvMKemRhtuBAGiX8SUGo9EcR/dJ\n87vMZkVy2g+pQ3nd/aKkd0s6JulLkj7q7g+a2U1mdtOw2P8m6fmSPmJmnzezP6ksYgAAUDFuBwKg\nbbjEoO+ynDGVux+VdDSy7PbQ3wckHSg3NAAA0AT3E0tms5IGjd0OBADW4hKDvsuUmAIAgI2F24Gg\ni7gGsc+4xKDvsszKCwAAALQatznqOy4x6DvOmAIAAKAHkq5BHBwUZ/87j0sM+o/EFAAAAD3ANYh9\nxyUG/cZQXgAAAPQA1yACXUZiCgAAgB7gGkSgyxjKCwAAgM7jGkSg20hMAQAA0Atcgwh0F0N5AQAA\nAACNIjEFAAAAADSKxBQAAAAA0CgSUwAAAABAo0hMAQAAAACNIjEFAAAAADSKxBQAAAAA0CgSUwAA\nAABAo0hMAQAAAACNIjEFAAAAADSKxBQAAAAA0CgSUwAAACAHs9mB2dxy8PfcstnsoOmYgK4iMQUA\nAACmFCShew5JR/cFS47uk/YcIjkF8iExBQAAAKY2syAd2b122ZHd0s6DzcQDdBuJKQAAADC1HVvi\nl2/fWm8cQD+QmAIAAABTO3M+fvnZc/XGAfQDiSkAAAAwtdVFaf7k2mUHTkmnDzcTD9BtlzYdAAAA\nANA17ieWzGYlDQ4Gw3fPnpNOH3Y/sdR0bEAXkZgCAAAAOQyTUBJRoAQM5QUAAAAANIrEFAAAAADQ\nKBJTAAAAAECjSEwBAAAAAI0iMQUAAAAANIrEFAAAAADQKBJTAAAAAECjSEwBAAAAAI0iMQUAAAAA\nNIrEFAAAAADQKBJTAAAAAECjSEwBAAAAAI0iMQUAAAAANIrEFAAAAADQKBJTAAAAAECjSEwBAAAA\nAI0iMQUAAAAANIrEFAAAAADQKBJTAAAAAECjSEwBAAAAAI0iMQUAAAAANIrEFAAAAADQKBJTAAAA\nAECjSEwBAAAAAI0iMQUAAAAANIrEFAAAAADQKBJTAAAAAECjSEwBAAAAAI0iMQUAAAAANIrEFAAA\nAJtftOgAAAoqSURBVJ1gNjswm1sO/p5bNpsdNB0TgHKQmAIAAKD1giR0zyHp6L5gydF90p5DJKdA\nP5CYAgAAoANmFqQju9cuO7Jb2nmwmXgAlInEFAAAAB2wY0v88u1b640DQBVITAEAANABZ87HLz97\nrt44AFSBxBQAAAAdsLoozZ9cu+zAKen04WbiAVCm1MTUzPab2UNm9rCZvTehzOLw+fvN7DXlhwkA\nAMzsJ8zsz83sO2b22gnlUo/dQNe4n1iSjt8sDZalGz4d/L53IVgOoOsmJqZmtknShyTtl/RKSTea\n2SsiZQaSdrv7NZJ+RtJHKoq1Ncxsb9MxlKEv2yGxLW3Vl23py3ZI/dqWDeoBST8q6Z6kAlmO3X3U\nl77dl+2QqtkW9xNL7ktz7h/bG/yuJyllv7RPX7ZD6te2FHFpyvPXSjrp7o9IkpndKekdkh4MlXm7\npDskyd0/a2bPM7Mr3P3JaGVm112QZJNXuU3Shaeku2X2tm9IF006JrPrvinpWemvz2KbpMvOSZ+Q\n2Rsfli69UtKW7HW/5BKz655Oq1uSzK47M13d0bpGbVHm9o/svsTs+vNr699mwTrPPy1tsrXrnCae\ntLLh59/yZelTkXLTtmN4nxRpt7jXbrNxf5m9RTpRcL/GrWvc16WX3G12nSlzu1fdT7JswyUXpItf\nW78fJ71XonXExV/WthVdR9p7Pvz6bZuly74tffUr0mdK6itpJn1ubtssXbg4ju8lW8yu8+piSYsx\n7rMlT12Xbchrytz9IUkym9h0WY7dfbRX0krDMZRhr/qxHRLb0lZ71Y9t2at+bIfUr23JLS0xvUrS\naujxo5Jel6HMiyWtS0yll6as70pJ+yTdMZxd7b07hjmvpJduS4k1oysl3Sbpnu3B45ftnlQ63tcl\nvXRTet2S9NLt68tlEW2LsrY/XP/jki4J1R9dZ9540sqOnj82fP4lV68tl6cdR/ukSLvFvXYUi4br\n/65fyhbPtOsK9/Xnb5ZeOiGmtJjrEt6P+7ZKdzxnfQxJ75VoHXHxl7VtZaxj0ns+/Po3bBv23c3S\nHcPPlqJ9Jc2kvvSGbcPnNg9j2TbclmpDSo2xaF2j92RN/3/pnizHbgAAWiUtMc36X/Xot4OE170w\npZoPSHp/qNzdob/TXpvVB0J1vylnvTsSXhetWznrH9UVbouytj9c/5u0Ns7oOvPGk1Z29Pxtkedf\nGHpemq4dR/ukSLvFvfYDkTLftWn6erOsK9zXw/0ra1tW1U8mCe/HpBiS3ivROuJeW9a2lbGOSe/5\n8OtHfbqKz64kk/pS3L5J2ydVmPTZkqeufjOzP1CQgUf9orv/foYqGjgjDgBAMeaefPwys9dLutXd\n9w8f/4Kkp939g6Ey/1HSirvfOXz8kKQ3R4fymhkHSgBAqdx9Q542NbM/kvTz7v5nMc+lHrtDZTk2\nAwBKlffYnHbG9D5J15jZ1QrGfb5T0o2RMndJerekO4cHw6/HXV+6Ub88AABQkaTjapZjtySOzQCA\n9pg4K6+7X1SQdB6T9CVJH3X3B83sJjO7aVhmSdJfmdlJSbdL+tmKYwYAYEMysx81s1VJr5f0STM7\nOlz+IjP7pJR87G4qZgAAspg4lBcAAAAAgKpNPGNahq7f5NvMHjGzL5rZ583sT4bLLjezPzCzvzSz\nu83seU3HGcfMfsvMnjSzB0LLEmM3s18Y7qeHzOz6ZqKOl7Att5rZo8N983kzmws918ptMbMZM/sj\nM/tzMzthZgvD5Z3bLxO2pVP7xcy2mNlnzewLZvYlM/u3w+Vd3CdJ29KpfRJmZpuGMf/+8HHn9ksb\ncWxuTl+OzX05Lkscm9u4LRyb27ktI5Udm929sh9JmySdlHS1pM2SviDpFVWus4Jt+LKkyyPL/r2k\n9wz/fq+kf9d0nAmxv1HSayQ9kBa7gpuwf2G4n64e7rdLmt6GlG25RdLPxZRt7bYomGnze4Z/b5f0\nF5Je0cX9MmFburhftg1/XyrpXkn/pIv7ZMK2dG6fhGL8OUn/RdJdw8ed3C9t+hHH5qZj78WxOWE7\nOvlZM+F41sX9wrG5ZdsxYVs6t09CMVZybK76jOkzN/l29wuSRjf57pro5BBv1/hGgXdI+pF6w8nG\n3T8j6e8ii5Nif4ek33X3Cx7clP2kgv3XCgnbIsVP/tHabXH3r7j7F4Z/n1Vww/ur1MH9MmFbpO7t\nl28N/3yWgi/tf6cO7hMpcVukju0TSTKzF0saSPoNjePv5H5pGY7NDerLsbkvx2WJY7Pauy0cm1u4\nLVUem6tOTONu8n1VQtm2ckmfMrP7zGx+uOwKH888/KSkK5oJLZek2F+kYP+MdGVfHTSz+83sN0PD\nBjqxLRbMmPkaSZ9Vx/dLaFvuHS7q1H4xs0vM7AsK2v6P3P3P1dF9krAtUsf2ydD/Jel/lfR0aFkn\n90vLcGxunz716y5+1jyDY3N7toVj8zNatS2q8NhcdWLah5mVfsDdXyNpTtK/NrM3hp/04Dx1J7cz\nQ+xt366PSPpHkr5H0hOS/sOEsq3aFjPbLunjkm529zPh57q2X4bb8l8VbMtZdXC/uPvT7v49kl4s\n6U1mdl3k+c7sk5ht2asO7hMz+yFJf+3un1fCbVG6tF9apg/twrG5nTr3WRPGsfkZrdgWjs1rq6g8\nyAyqPjb//+3dv2oUURTH8e8pDKgIIoIKpkhhaaWVjSD4rxHs0kjQd7AwT+EL2BghpSFgI+IDWBiN\nIiKClYqxsrXwWMwdXWU3Rsgw94bvB5Ykk2z2/Dize7hkcnfohelHYH7i63n+XDVXLzM/l49fgYd0\nf37+EhHHASLiBLA1XoX/bVbtf/fqZDlWrczcyoLucoL+0oCqs0TEPrrBt5KZa+Vwk32ZyPKgz9Jq\nXwAy8xvwCDhDoz3pTWQ522hPzgHXIuIDsApciIgVGu9LJZzN9dkT53WjrzWAs3ni7lVlAWczdWUZ\ndDYPvTD99SbfETFH9ybf6wM/5q6JiAMRcah8fhC4BLyiy7BUfmwJWJv+G6o0q/Z1YDEi5iJiATgF\nPBuhvh0rJ37vOl1voOIsERHAPeBNZt6d+FZzfZmVpbW+RMTR/vKZiNgPXAQ2aLMnU7P0w6KovicA\nmbmcmfOZuQAsAk8z8wYN9qVCzub67InzurXX/56zub4szuY6sww+m3P4XZuu0u0I9h64M/Tj7XLt\nC3Q7Sb0AXvf1A0eAJ8A74DFweOxaZ9S/CnwCvtP9P9HN7WoHlkuf3gKXx67/H1luAfeBTeBleQIc\nqz0L3S5sP8o5tVFuV1rsy4wsV1vrC3AaeF5ybAK3y/EWezIrS1M9mZLrPL93/muuLzXecDaPWf+e\nmM1TcjQ5l0ttzubKsmwzz1rsibN5h1mi3EGSJEmSpFEMfSmvJEmSJEnbcmEqSZIkSRqVC1NJkiRJ\n0qhcmEqSJEmSRuXCVJIkSZI0KhemkiRJkqRRuTCVJEmSJI3KhakkSZIkaVQ/AeTK6yUBnozGAAAA\nAElFTkSuQmCC\n", "text/plain": [ "" ] }, "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": 19, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/condatl/anaconda/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": [ "[]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiEAAAFwCAYAAABuLatIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeUZVWVgPFv2w0okiSnVhglR9FBJDYCiiJRkgiSRIEh\nyIiAqCRzwoQgkkRyUrQV0VZpRFByhpYwIKElSBAFJO7549zWoqj87nv3VdX3W4tF1U1nV13ot/uE\nfSIzkSRJ6rTXNB2AJEkan0xCJElSI0xCJElSI0xCJElSI0xCJElSI0xCJElSI1pOQiJiUkRcEhG3\nRsQtEbFfHYFJkqSxLVqtExIRCwMLZ+YNETEHcC2wRWbeXkeAkiRpbGq5JyQzH8rMG6qv/wncDiza\n6nMlSdLYVuuckIhYAngrcGWdz5UkSWNPbUlINRRzPrB/1SMiSZLUr4l1PCQiZgEuAE7PzAt7nXNz\nGkmSxpDMjDqeU8fE1ABOBR7LzAP6OJ91Bat6RMQRmXlE03Go8H10H99Jd/F9dJc6P9frGI5ZC9gR\nWD8irq/+2biG50qSpDGs5eGYzPwDFj2TJEnDZPIwPk1rOgC9wrSmA9CrTGs6AL3CtKYDUHu0PCdk\n0AacEyJJ0pjRbXNCJEmShs0kRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIk\nNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIk\nRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIk\nNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIk\nRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNcIkRJIkNaLlJCQi\nTo6IhyPi5v6vYaNW25EkSWNLHT0hpwAbD3LNVjW0I0mSxpCWk5DMvAx4YpDLlmu1HUmSNLZ0ak7I\ngh1qR5IkjRKdSkIW6FA7kiRplJjYmWYOny/i80fCywlMy8xpnWlXkiS1IiImA5Pb8uzMbP0hEUsA\nUzJzpT7OJeTjwIqlPWa03KAkSWpERGRmRh3PqmOJ7lnAFcDSEXF/ROzax2WPAt8BHmy1PUmSNDbU\n0hMyYAOlJ2QasDowOzAxk5fa2qgkSWqLruoJGaL7KQkIwFIdalOSJHWxTiUhD1f/vhT4rw61KUmS\nulinkpCZq3DuBhbrUJuSJKmLdSoJuQKYATwALN6hNiVJUhfrSBKSyXmZLEZZHWMSIkmSOtYTMtP9\nwBs73KYkSepCnU5C/gws2+E2JUlSF+p0EvIXYL4I5upwu5Ikqct0NAmpipTdDqzQyXYlSVL36XRP\nCMDlwLoNtCtJkrpIE0nIb4H3NNCuJEnqIh3ZO6ZnjfkIXgf8FVg2k4fa2rgkSarVaNw75t8yeRaY\nAmzd6bYlSVL3aGI4BuBs4IMNtS1JkrpAx4djyjFmpZRxXy2T+9oagCRJqs2oHo4ByOR5YCqwfhPt\nS5Kk5jU1HANlU7s1G2xfkiQ1qOkkZK0G25ckSQ1qZE5IOc4swOPApEyebGsQkiSpFqN+TghAJi9Q\nqqdu0lQMkiSpOU0OxwAcB+zVcAySJKkBjQ3HlHMu1ZUkaTQZE8Mx8O+luhcCO0UwT5OxSJKkzmp6\nOAbgRGA/YEaEVVQlSRovGh2OeeV1vAv4diYrtTUgSZI0YmNmOKaXacD8Ebyl6UAkSVL7dU0SksnL\nwCXA2k3HIkmS2q9rkpDKtcBqTQchSZLarxuTkLc1HYQkSWq/rpmYWq7lDcBfgHmq4RlJktRFxurE\nVDJ5AngSWLLpWCRJUnt1VRJSuRFYpaqmKkmSxqhuTEJuAjYFnotgz6aDkSRJ7dGNSciNwC7AH4CP\nNxuKJElql25MQqYAnwS2ARaKYOGG45EkSW3QdUlIJs9m8vVMHqL0hqzbdEySJKl+XZeE9HIpsF7T\nQUiSpPp1exLyeyzjLknSmNRVxcpefS+zAY8DC2TyTL2RSZKk4Rqzxcp6y+Q54FbgHU3HIkmS6tXV\nSUjlLGDfCGrJuiRJUnfo6uGYcj+vp6ySeRTYNpMnawtOkiQNy7gZjgHI5GlgdcrckL0aDkeSJNWk\n63tC/vMcVgfOBJZ2h11JkpoxrnpCergaeBaLl0mSNCaMmiQkkwROBD7SdCySJKl1o2Y4pjyL+YC7\ngSUzeaKOZ0qSpKHrquGYiNg4IqZHxJ0RcXAdQfUnk8eA04GfRjBPO9uSJEnt1VJPSERMAP4MbAg8\nSJm38cHMvL3HNbVlTOV5vAY4DlgY2KIappEkSR3QTT0hqwN3Zea9mfkCcDaweeth9a9aGbMvsCjw\n0Xa2JUmS2qfVJGQx4P4e3z9QHWurTJ4HdgE+H9H+9iRJUv1aTUIaGwrJ5Fbgu8BJlnSXJGn0mdji\n/Q8Ck3p8P4nSG/IKEXFEj2+nZea0Ftud6UvA5ZRKqseWtlgMeCSTF2pqQ5KkcSsiJgOT2/LsFiem\nTqRMTN0AmAFcRZsnpr46BpahJCJrA58GtgYuzOSD7WpTkqTxqmsmpmbmi8A+wK+A24BzeiYgnZDJ\nn4HDgVuBN1FWzWwUwRKdjEOSJA3PqCpW1n8bBLAaMD2TpyM4Ebg5k2+3s11JksabrukJ6RaZZCbX\nVjvuAvyONo1fSZKkeoyJnpBXt8kilOGZ+d1xV5Kk+tgTMohM/go8AqzSdCySJKlvYzIJqUzDIRlJ\nkrrWWE5CLgXWnflNBBHBhyK4KIKNGoxLkiQxtpOQS4DJEcxWff9JSh2RnwNnRrBSY5FJkqSWK6Z2\nrUweiuBmYJMIXgL2A9bI5IEIJgBfBDZtNEhJksaxMbk65j9tsxXwfUqPz/syuao6/lrgLmDzTK5t\nIjZJkkajOj/Xx3oSEsAOwN2Z/KnXuQOBlTP5cBOxSZI0GpmE1CCCeYG7gaWBF4A9gGmZXN1oYJIk\ndTHrhNQgk8eBc4FTgeuANYCLIliv0cAkSRonxm0SUjkUmA7snskHgO2B8yPYsdmwJEka+8btcEx/\nqqW7PwO+mMkJTccjSVI3qfNzfcwu0R2pTG6uipldGsGDmVzUdEySJI1F9oT0I4KNgW8BK2byYtPx\nSJLUDZyY2hm/Av4K7DLQRRG8JYLPRVj4TJKk4TAJ6UcmCRwEHBXBfH1dE8HewJXAQsDxEby3gyFK\nkjSqORwziAiOBhbLZLtex18H3Aesk8n0CNYCfgwsncnfGwhVkqS2czimsw4F1olgxV7HNwRuzWQ6\nQCaXA78GPtrh+CRJGpVMQgaRyb+A44F9ep3aHPhpr2PHA7tU5eIlSdIAHI4ZgggWBv4MLJTJvyJ4\nDTADWDOT/+txXVBKwW+dyXXNRCtJUvs4HNNhmTwE3AJMrg69A3i0ZwJSXZfAaeCmeJIkDcYkZOgu\npswDgb6HYmY6DfhgBLN0JCpJkkYpk5ChuxRYrxpy2Qa4oK+LMrkLuAt4TwdjkyRp1DEJGbqrgOWA\ndYEEbhjg2h8BO3UiKEmSRiuTkCGqVslcA5wEnFfN/+jPucBGESzZkeAkSRqFXB0zDBFsD5wFLFJN\nVh3o2sOApTJf2SMSweuBjwCzAZcA1wyS0EiS1DXq/Fw3CRmmCGIoSUMEc1OW664+cxVNBHMAUyl7\n0txPmTcyP3AesLfJiCSp27lEt0FDTRSq0u3fp+w/QwSzU1bU3AZ8IJP9KXNM3gGsA+47I0kaX+wJ\naaMIFqAkHXsAewOPADtn8lKv67YDPk4pfmZviCSpa9kTMkpk8iiwLfBNyhDMLr0TkMr5wCLAKh0M\nT5KkRtkT0iUi+AIwWyYHNh2LJEn9sSdkbDoD2CGCCQNdFMEWERwbwdIdikuSpLYwCekSmdxGGbKZ\n3N81ERwCHA38E7i4WoEjSdKoZBLSXc4AduzrRARbAHsBa2VyEPBr4MsdjE2SpFo5J6SLRLAIZTXN\nGzP5R4/jbwBuBbbO5Irq2PzAHcCKmcxoIl5J0vjjnJAxKpO/Uno49uh16qvAT2YmINW1f6PsUfO/\nnYtQkqT62BPSZSJ4KzAFeHMmz0WwLPB7Sgn4v/e6dhJwI7Bk73OSJLWDPSFjWCbXA7fwn114Dwa+\n21eSkcn9wO8otUgkSRpV7AnpQhFMBk4ANqf0giydyeP9XLspcHAma/dxbg5gO2AlYBJlv5rvzNzL\nRpKk4bInZOy7lDLp9FbgqP4SkMrFwNIRLNnzYATLAX8GNqEkH+cA/wCuiGCpoQQRwSwRfKiqTeJ/\nK5KkWtkT0qWqDe9WBK4ebD+ZCE4Abs/k6Or71wBXA8dn8oNe1+4PbJLJuwd5ZgBnU3pQZgcurTbd\nkySNY/aEjAOZPJPJVUPc0O7HwFY9vt8USMqQTm/HActXE2AHciDwRmB9YF1g2whWGkIskiQNiT0h\nY0AEs1Gqra4MPAj8CfhqJhf0c/2BwGqZ7NDP+YUp9UpWzeS+6tihlBU7u7fhR5AkjRJ1fq6bhIwR\nEXwPeAq4BPg2sEImL/dz7dzAPfRIMnqd/zowayb79Tg2P3AnsEwmj7ThR5AkjQImIXqVqmbINcCs\nwK6ZXDjI9UcDmckneh1fgDKhdeVMHuh17gfA/Zl8rtbgJUmjhkmI+hTBosCimVwzhGtnFjpbuqq+\nOvP4F4E3ZLJXH/esAPwGWCKT5+qLXJI0WnTFxNSI2CYibo2IlyJitTqCUWsymTGUBKS69n7gfOCA\nmcciWBD4GP1sjJfJrZRlwx9oPVpJ0njXyuqYm4EtKcW0NDp9HtgzgpWr7z8DnJHJXwa45/uURKVP\nESwfwT4RLF1jnJKkMWjESUhmTs/MO+oMRp1VTUrdE/h1BN8FtqEkJgP5KaU42nK9T1RLeC8F3gZc\nHsEWNYcsSRpDrBMyzmVyHqU8/KPAuwZb+ZLJC8DJlOTl36oJrVOA/TPZlVKp9fgI5mpL4JKkUW/A\niakRMRVYuI9Th2bmlOqaS4BPZOZ1/TwjgSN7HJqWmdNGHLEaF8HilEmty2fycHXsDOChnqttIjgV\n+EsmhzUTqSSpVRExGZjc49DhXbM6ZihJiKtjxp4Ivg1MyGSfCDakVGddMZOne1zzJuA6eiQrkqTR\nrStWx/RikjH+fA7YPIIDgBOB/+mZgABUE1x/BHy2gfgkSV1uxD0hEbEl8B1gfuDvwPWZ+d4+rrMn\nZIyKYA3KXjQnZHJsP9fMD0wH3pnJnZ2MT5JUP4uVaVSJ4CDKeOImQ9yQT5LUpUxCNKpEMCtwPXAu\n8F3KXJJHm41KkjQS3TgnROpXJs8DGwNvpezye1cEF0QwR7ORSZKaZBKijsjk/ky2yOR1wILAE8Cl\nEXwwgjc0HJ4kqQEmIeq4avO7PYBTgO2B+yK4LYJdIvxvUpLGC+eEqHHVnJH1gKMolVv3z+SeZqOS\nJPXFOSEaUzJ5PpOpwPqUCazXRHBMBG8c7N4I1o/g1Ai+MJTrJUndwyREXSOTf2VyOLAs8CxwfQRn\nRDB3X9dH8FHgTOBKYFbgqghW7FjAkqSWOByjrhXBnMCXgXcBm2dyR3X8NcDhwIeAjTO5qzq+Q3X9\nqpk8PoTnzw0cAKxOSXq+X/XISJL6YZ0QjStVj8c3gIeBm4DlgMeBLXvv+htRqvhmssMgz3w7cB5w\nCXAh8AZKKfqTMzmi7p9BksYKkxCNOxFMBJYCVgH+AlyZyct9XDc7ZdO8wzM5p59nrQFMAfbM5IIe\nxxcA/gh8IpOf1v9TSNLoZxIiDSCC/wZ+DiybyRO9zi0AXAvs21eiEcHalMquy2fyZCfilaTRxCRE\nGkQExwNPZ/K/PY5NAH4JXJfJIYPc+1wm+7U/UkkaXUxCpEFEsBBwO7ByJg9Ux44C1gE2yuTFAe6d\nD7gNeE8mNwyxvVmAPYH3AVcDRw3UhiSNVtYJkQaRycPAycBBABFsAuwGbD9YcpDJY8CBwDkRzDNY\nWxEEcD6wKXAi8E7gpOq4JKkfE5sOQGqjrwO3RHAP8CnKMt+Hh3JjJqdVK2jOiuD9mbw0wOVbAm8C\n/juTFyL4JfAnYGvKChxJUh/sCdGYlclDwF7Ae4GtMvnjMB9xIPBa4Mj+LqhW7XwBOCSTF6p2nwH2\nBb5WlaSXJPXBOSHSACJYkLLkd6dMLunj/G7Ah4H1M8le5y4Bjs20N0TS2OGcEKlDqmJoewA/7F0+\nPoLXAUdQekH6yuZPqO6VJPXBnhBpCCI4DnhdJrv0OHYwsEYmW/Zzz2uBGcCKmcwYYjuzARMzebr1\nqCWpfvaESJ13ILBWREk4IlgS+GR1vE+Z/Av4FbDJUBqI4EPAQ8AjEdwcwacjWLTlyCWpS5mESENQ\n9Ux8GDgugmWB44CvZ3L3ILdOAd4/2PMj2Bj4GrA2MA9lGGcScGsEp5uMSBqLHI6RhiGCXYEfUPaY\n2WDmipgBrp8XuAdYOJNn+7lmInAzcEAmF/c6NydlefEewA7u8iupaVZMlRpUzfV4vq8N9Pq5/lLg\nK5lc1M/5HYGP0McKmx7XrANcAGydye9HFrkkta7Oz3WLlUnDVM31GI4plGqqfSYhlATkO/0lIFWb\nl1VzRs6N4N2Z3DTzXAQbABsB/wTuBe4GrhqkwJokNc6eEKnNIliOMkH1TX3UElkSuApYPJPnhvCs\n7YDvAocCE4BtKNVafwTMDiwBrEiZ7/XJ/npfJGmkHI6RRpFqD5m7gC0yubnXuc8CC2WyzzCetzZw\nAPAEcDlwes+5KVV77wO+CdwBfDyTu1r+QSQJkxBp1IngO8BfM/lSj2NBSRI+lMlVbWhzVuDjwCeA\nzTK5su42JI0/1gmRRp9f8Op6IesAzwNXt6PBTJ7P5KuUOSc/iWChdrQjSSNlEiJ1xqXAShHM1+PY\n7sBJA01IrUMmU4CTgVOq3hdJ6gomIVIHVCtqLgE2BohgAWAz4LQOhXAksAiwUzseHsFaEfwsgh9E\n8MZ2tDFA27NW7S/eyXYltc4kROqcC4Ftq6//Bzgvk0c70XA1cXU34OsRLNLfdRFMiGCfCM6P4GPV\nJn0DqnYSPp/y8z0EXB3BusONcSS9NFXicSNwDHBDBHsN9xmSmuPEVKlDIpgD+AslGTgBWCeTP3c4\nhs8DKwBb9bFceAnKUt+XKcM3WwNvo6yyObqv4mwRrA+cCayXyR3VsQ2Bs4BNM/nTIPFMpvw+5qfU\nOvkNsOVQarFEMDvwJ+CMTL5SLXe+glLQ7fLB7pc0Mk5MlUahTP5JWa1yHPD5Ticglc8BbwG+GsH7\nIvivCOatlgpfC/ycUo7+R5lsBrwX2Ao4OeKVf15EMA/wQ2C3mQkIQCa/AXYGLoxg+b6CiGC2CL4B\nnE4pgX8SsCjwFHDsEH+WQ4DpwFerdu+h9DCdXK0MktTl7AmRxplqzsahlGRkeWBuylDKZ6oP8t7X\nz04ptnZlZtk1uBo6ORN4IpO9+2lnR+CLwNqZ3Nfj+GLAT4AZwO6ZPNbj3OuBW4CdBypPX23odzOw\nciYP9jp3EXBxJt8Z7HchafisEyKpFlUyMUsmzw9y3bzA7ynDHwcC+wAfANbK5JkB7vs48DFK78qM\nCNYCzgG+B3y5r5VBEXwUeG8mWw7w3MMpRd5elQBFsBJlWGeJ/jYNlDRyJiGSOi6CuSjzQ3ajVGrd\noWcPxwD3HQwcBrwIPA3smcnPBrh+5tyZVTO5v4/zEyl75Lyv5x46va75BXBuJqcOFp+k4TEJkdSY\nCCZm8uIw75kAzAM8PpS6KFWF2X9k8uk+zm0JHJjJWgPcvwlwWCbvGE6ckgZnEiJpTItgWUqBtzf2\n3tgvgl8DP8rk9AHun0DZr2ebTK5pa7DSOOPqGEljWibTKfU/tul5PIKlgFUpdUkGuv8lyjLoj7Qr\nRkmtMwmR1K2OAfbtdWxv4IdDqSNCWf67dQSz1R6ZpFqYhEjqVr8AFowo8zqqSq87A98ays3VpNlb\nKLVOJHUhkxBJXakaUvkKpdT8ayi1TU7JZMYwHnM6sGM74pPUOiemSupa1QTTyylLcjcClsvkkWHc\nP09175KZPNGOGKXxxompksaFqjdkS+B+4P3DSUCq+58EptJrgquk7mBPiKQxLYItgP0zWb/pWKSx\nwDohkjRE1eqYGcAqmTzQdDzdJoI3A+sDCVwG3DmUgnIav7pmOCYivhYRt0fEjRHx44iYu46gJKku\nVbGznwDbNx1LN4lgwQjOpOxivA4lEbkEuDaCjRsNTuNGq3NCfg2skJmrAHcAn2o9JEmq3RnADk0H\n0S0iWAa4GniAMml350x2BCYBhwPHR/C9ap8eqW1aSkIyc2pmvlx9eyWweOshSVLtfg8sFMFyTQfS\ntAjmo9Rg+XwmB2Xy9MxzmbycyRRgJWAp4MwIZm0oVI0Dda6O2Q24qMbnSVItqlU25wDbDee+COaN\n4NQIbo/gzghmRPBkBFMieEt7om2fCILSK/STTE7o77pMngI2A2YDzq3qtEi1G/Q/rIiYGhE39/HP\npj2u+TTwfGae2dZoJWnkfgq8f6gXR/A24BrgCeADwCbA6sDSwO+AK6qVN6PJVsBiDGHovCqNvzUw\nL/CZNselcarl1TERsQuwB7BBZr5qP4eISODIHoemZea0lhqVpGGKYBbgYWDFwaquVgnIxcDemZzX\nzzVrAGcBNwPTgL8BbwAWAR6l7OI7C/Au4L+AiZTN9yZQeo2PyOTOln+wIYpgduA2YNdMLhnGfYtQ\nkrEPZTKtTeGpi0XEZGByj0OHd8US3YjYGPgGsF5m/q2fa1yiK6krRHAW8NtMThzgmjkoicUnMwfe\nrTeCuSi9C28F5gMepyQ6CwBvoSx7vRS4tfr6JuAlyu6+B1B2Az5qmKXoRySCI4DlM9l2BPduAXwR\nWDWT5+uOTaNL19QJiYg7gVkp/+MB/DEz9+51jUmIpK4QwQ7Atpn9D6NE8GVg0Uw+3OZY5gMOBnYH\nTgWOyeT/2tTW/JQVjKtlcu8I7g9K783vMvlazeFplOmaJGRIDZiESOoSEcwL3AMsVM156H1+UcrO\nuytk8tcOxbQYsD+wK3AB8OlMHqu5jc8DC2by0RaesRSlpsgqmTxYW3AadbqmWJkkjSaZPA7cyCvH\nt3vaDzitUwlIFdODmRxEmfD6InBzBDtXw0ItqxKvvSjDKSNWzV/5PvD1OuKSwCRE0vjzc2DT3ger\nIYuPAN/qeERAJk9ksg+lsuv2wPQI1qrh0fsDF45kGKYPXwTWjHAfHtXD4RhJ40oEywO/BJbouUdK\nBN8Hns9kv8aC6yGC9wE/BDbI5OYRPmNByoqY1euabxLBVsDnKJNUX6jjmRpdHI6RpJG7nbJCZcWZ\nByJ4D6V35PCmguotk4soK2gujGCk+3J9Fji95gmvP6GUe9+3xmdqnLInRNK4E8G3gScyOSKCJYA/\nAdtkclmzkb1aBN8DFqSs6hnyH9jV/jB/AJbLpM8SCi3EtBylFP4KmTxS57PV/VwdI0ktiGBFYCqw\nCmXp6ZmZHN1sVH2L4LWUZOIa4GjgWUrRs42AuYCjM7mp1z0Tgcspk2yPaVNcXwZWA7bI5Jl2tNEp\nVQn+JYCngGurMv/qh8MxktSCTG6hFAq7D7gb+GazEfWvWkq8MfACZefyq4H/oVRovQWYGsExvYZs\nPg38HTi2jaF9FngE+F0192TUiWC+CKYClwGHAicBj0ZwXgSbV/VR1Eb2hEgal6oPmOWB6aP5b75V\n0bMvUXpGTgaWBDYA1mx3PY/qd3gUZVXRncAKlNL1zwLXUfbdeQ44l7JCp2smskawAPAb4FeU2iwv\nVMcXBd4NHATcAOxildhXcjhGkvQKEWxA+fB8mDIM82gH214BWIiyEudRyjDRW6t/v4GSpDxNlwzd\nRLAQ8FvKJNvD+pprUw2DnQ08ncmHOhxiVzMJkSSNGtUclVOB2YEPZPJyg7EsTNkF+exMjhrk2tdR\nekM+099GhuORSYgkaVSJYDZK78OPm5oEXO2kfCkwNXNoy7GrwmzHU1YZjdphuzo5MVWSNKpk8hyw\nG/CpBieyfoUyoffIYdwzjbJJa7+bHmrkTEIkSR2RyR2UeRaf7HTbEXyAkkjsPJzhoGq+yDcp+++o\nZg7HSJI6JoIlKcuM35TJ0x1qcylK3ZT3ZXLNCO6fHZgBLJPJw3XHN9o4HCNJGpUyuYdSl+PDnWiv\nWkZ8LPClkSQgANWKnp8D29QZm0xCJEmd9y1g/4iOfAZtAiwOLVeOvYA+dl9Wa0xCJEmd9nvgeUpR\ntbapan18A/hEDYXSLgHWrJ6pmpiESJI6qprseQywT7vaiGACcCJwI/DLVp+XyZOUYmzvbPVZ3SSC\nJSI4J4JnI7g6gq062b5JiCSpCWcA67VjuW4E81KGTxamlF2vawXGb4ANa3pW4yJYD7gKuAmYBBwG\nfCuCnToWg6tjJElNiOAs4PeZHDfM+/6bsrfLvJS9aS4BHgTeCKwLrA38EPhkVZ+krnjfDRyayeQW\nnzMbcDCwH2Xn3s8CDwBXVhsWtl0EK1Eqx26Xye96HF+OUtBt40yu6/teK6ZKkka5CDYHDhjOh3o1\nXPB9ygf33cCclM375gPuB64FfpHJU22Id56qjTdk8uIIn7Eg8DPKHjsHUPbc+SYwN3BjJtvWFO5A\nMUwE/gh8P5OT+ji/K2W/n7X73lfHJESSNMpVkzwfBt4ylA33IngL5cPzPf39Lb3dIrgN+FAm14/g\n3gnAr4HrKb002ePca4HbKcXUfl9XvP3EsSewPbB+P5v3vQaYDuyeyWWvPm+dEEnSKFcNPUylLKMd\nULXvy1nAkU0lIJU/AWuM8N6DgInAIb0//KvfxWeAL7QW3sAimAs4nNID1WcvRFVR9mRKotJWJiGS\npCb9DNhsCNd9BngE+F57wxnUiJKQCBYFDgR2HGAo5xxgUjXnpV0OAS4eQk/OT4At2l3LxeEYSVJj\nIpifMrdj0f7KuEfwJuA6YMVM/trJ+PqIZSXg/EyWGeZ93wWez+QTg1z3v8DbM9mhhTD7e/Yk4AZg\n5UweHMKAJeXVAAALIElEQVT1t1NWF135yuMOx0iSxoBM/kbpXRhoSOZjwGlNJyCV24BFqmXAQxLB\nnMBOwNeGcPlJwHsiWHyE8Q3kQODkoSQglV/Q5iXJJiGSpKadA2zX14lqOevuMLxlvO2SyUuUAmir\nDuO2bYFpmTw0hOf/HTiNmgu5VatydgKOHsZtl1OWO7eNSYgkqWk/ATasegx62xq4KZM/dzimgdzA\n8JKQnYFThnH9d4Hdai4R/3Hg7GH2Jl0BvLOd80JMQiRJjcrkCcrOuq/YpbbaAXd/Wt98rm43AqsM\n5cJq2GZV4OKhPjyTuymJzqAl1COICN4fwZQIfhfBZtXvrec1C1GGtL461BiqOB6m1DNZYTj3DYdJ\niCSpG3wB+HyvuRZrUqqi/ryZkPo1nJ6QdwOXjqBy6/HAnkO47kuUuSbnUVYOfRU4sVev0uHAGZnc\nO8wYoAzJrDmC+4ZkYrseLEnSUGXyxwjOp8xZ2KX62/xhwNHVPIxuciuwVASzZvL8INduzMg20PsZ\n8N0IVsjk1r4uiGAbSm/JmtUEXyL4FWX+zIwIrgYmAG8C3jaCGACuAVYb4b2DcomuJKkrRDAH5UPv\ngeqfd1CWk77QaGB9qCqnbp/JTQNc8xpgBiVJ+L8RtHEYMCmTPfo4Nyelwur2mfyhj/NzAJOB54Er\nMvnncNuvnrMm8J1M3v6fY5ZtlySNQdVwzFrABpRludc2HFKfql6b8zI5Z4BrVgPOGm5NkR73LwDc\nASxbzc/oee5wYKlMdhzJs4cRw1yURGrOmRVW6/xcdzhGktQ1MnkcmFL9081uA5Yf5JqRDsUAkMmj\nEZxCmfex28zjEcwH7AusPtJnDyOGpyJ4GlgY6q/T4sRUSZKG73YGT0LeyzBWxfTjcGCjCNbtcewz\nwDkjGeIZobuAt7TjwSYhkiQN323Acv2djGAeyjLeS1tpJJN/UHo9LozgmAh2ohQdO7KV5w6TSYgk\nSV3kDuDN1e6+fdkQ+EMmz7baUCYXUpYEP0pJQD6cySOtPncY7gKWaseDnRMiSdIwZfJsBA8Abwam\n93FJHUMxPdu7j872fvR0F7BlOx5sT4gkSSPT57yQqsZJS5NSu8ydOBwjSVJX6W9eyErAs5nc2eF4\n2uVeYIl2PNgkRJKkkelvmW6tQzFd4DHg9TVvqAeYhEiSNFL9LdMdS0MxVEXKHqbUCqmVSYgkSSMz\nHVgmggkzD1QVRt8OTGsqqDaZASxa90NNQiRJGoGqhscjlBUyM20A/DGTp5uJqm3+CixS90NNQiRJ\nGrkrgXf2+H5MDcX08BAOx0iS1FUup2y4N3Np7liblDrTY8B8dT/UJESSpJH7dxJCmaT6Mn0XLxvt\n/kY3JSER8bmIuDEiboiI30bEpDoDkyRpFLgRmFTtbLs58IuZW96PMY8B89f90FZ6Qr6amatk5qrA\nhZSd/iRJGjcyeRH4FbAdZV+Xs5uNqG3a0hMy4r1jMvMfPb6dgxKgJEnjzbeBy4DrgD80HEu7tKUn\npKUN7CLiC5TM7xlgjVoikiRpFMnkDxGsDDw0RodioHQ01J6ERGb/v6+ImErfS3IOzcwpPa47BFgm\nM3ft4xmZmVFHsJIkqfMimBe4K5N56/xcH7AnJDM3GuJzzgQu6u9kRBzR49tpmTltiM+VJEkNiojJ\nMPFdcOg8EUcdUeuzB+oJGSSopTLzzurrfYHVM3OnPq6zJ0SSpFEugmeA+SGe7khPyCC+FBHLAC8B\ndwN71RGQJEnqSk8Bc9X5wFZWx2xdZyCSJKmr1Z6EWDFVkiQNhUmIJElqhEmIJElqhEmIJElqxFPA\n3HU+0CREkiQNxTPA7HU+0CREkiQNxb+A19b5QJMQSZI0FCYhkiSpEf8CZqvzgSYhkiRpKOwJkSRJ\njTAJkSRJjTAJkSRJjTAJkSRJjXgOkxBJktQAe0IkSVIjTEIkSVIjrBMiSZIaYU+IJElqhEmIJElq\nhEmIJElqhEmIJElqxD+Aq+p8YGRmnc97dQMRmZnR1kYkSVJH1Pm5bk+IJElqhEmIJElqhEmIJElq\nhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmI\nJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElqhEmIJElq\nhEmIJElqhEmIJElqhEmIJElqRMtJSER8IiJejoh56whIkiSNDy0lIRExCdgI+Es94agTImJy0zHo\nP3wf3cd30l18H2NXqz0hRwMH1RGIOmpy0wHoFSY3HYBeZXLTAegVJjcdgNpjxElIRGwOPJCZN9UY\njyRJGicmDnQyIqYCC/dx6tPAp4B397y8xrgkSdIYF5k5/JsiVgR+CzxTHVoceBBYPTMf6XXt8BuQ\nJEldKzNr6XgYURLyqodE3AO8LTMfbz0kSZI0HtRVJ8TeDkmSNCy19IRIkiQNV1srpkbExhExPSLu\njIiD29mWioiYFBGXRMStEXFLROxXHZ83IqZGxB0R8euImKfHPZ+q3tH0iHh3/0/XSEXEhIi4PiKm\nVN/7PhoUEfNExPkRcXtE3BYR7/CdNKf6/d4aETdHxJkRMZvvo7Mi4uSIeDgibu5xbNjvICLeVr3H\nOyPi24O127YkJCImAMcAGwPLAx+MiOXa1Z7+7QXggMxcAVgD+J/q934IMDUzl6ZMKj4EICKWB7aj\nvKONgWMjwnL+9dsfuI3/DF36Ppr1beCizFwOWBmYju+kERGxBLAHsFpmrgRMALbH99Fpp1B+nz0N\n5x3MnKh6HLB7Zi4FLBURvZ/5Cu18casDd2XmvZn5AnA2sHkb2xOQmQ9l5g3V1/8EbgcWAzYDTq0u\nOxXYovp6c+CszHwhM+8F7qK8O9UkIhYH3gecyH+Wsvs+GhIRcwPrZObJAJn5Ymb+Hd9JU56i/OVp\n9oiYCMwOzMD30VGZeRnwRK/Dw3kH74iIRYA5M/Oq6rof9binT+1MQhYD7u/x/QPVMXVI9TeMtwJX\nAgtl5sPVqYeBhaqvF6W8m5l8T/X7JvBJ4OUex3wfzVkSeDQiTomI6yLihIh4Pb6TRlSrKr8B3EdJ\nPp7MzKn4PrrBcN9B7+MPMsi7aWcS4ozXBkXEHMAFwP6Z+Y+e57LMRh7o/fjuahIR7wceyczr6aeg\nn++j4yYCqwHHZuZqwNNU3cwz+U46JyLeDHwcWILyITZHROzY8xrfR/OG8A5GpJ1JyIPApB7fT+KV\nGZLaJCJmoSQgp2XmhdXhhyNi4er8IsDMonK939PMwnOqx5rAZlUtnbOAd0XEafg+mvQAZcuJq6vv\nz6ckJQ/5ThrxduCKzHwsM18Efgy8E99HNxjOn1MPVMcX73V8wHfTziTkGsqklCUiYlbKJJaftbE9\nAdXkoJOA2zLzWz1O/QzYufp6Z+DCHse3j4hZI2JJYCngKlSLzDw0Mydl5pKUyXa/y8yd8H00JjMf\nAu6PiKWrQxsCtwJT8J00YTqwRkS8rvrza0PKJG7fR/OG9edU9f/WU9VqswB26nFPnwbcO6YVmfli\nROwD/Ioy2/mkzLy9Xe3p39YCdgRuiojrq2OfAr4MnBsRuwP3AtsCZOZtEXEu5X/6F4G90+Ix7TTz\nd+v7aNa+wBnVX5DuBnal/DnlO+mwzLwxIn5E+Yvry8B1wA+AOfF9dExEnAWsB8wfEfcDhzGyP6f2\nBn4IvI6yAu3iAdv13UmSpCa4tlqSJDXCJESSJDXCJESSJDXCJESSJDXCJESSJDXCJESSJDXCJESS\nJDXCJESSJDXi/wH45KNQYxUmCAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(log10(En_array-En_array.min()))" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Convergence is much slower in this setting." ] }, { "cell_type": "raw", "metadata": {}, "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.4.5" } }, "nbformat": 4, "nbformat_minor": 0 }