{ "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": [], "source": [ "using PyPlot" ] }, { "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": [ "srand(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": [ { "data": { "text/plain": [ "100" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P = Int(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": [ { "data": { "text/plain": [ "100×400 Array{Float64,2}:\n", " 0.0679107 -0.0682832 -0.050376 … 0.29911 0.254127 \n", " 0.0828413 0.0156483 0.00646768 0.0065698 0.202104 \n", " -0.0353007 0.0140714 0.0907637 -0.0912724 0.0744238 \n", " -0.0134854 -0.0126612 0.0720839 -0.0435621 0.166798 \n", " 0.0586617 -0.116898 -0.0686564 -0.103057 0.108181 \n", " 0.0297336 0.0634067 -0.0237696 … -0.0351895 -0.080147 \n", " 0.00649475 -0.099871 -0.103289 0.1587 0.132639 \n", " -0.0109017 0.102912 -0.180055 0.0525735 -0.118446 \n", " -0.051421 -0.115751 -0.116446 -0.183929 -0.00756541\n", " 0.157433 0.0533604 0.121158 0.000143597 0.0405142 \n", " -0.0688907 -0.0308418 0.00927033 … 0.107258 0.054582 \n", " -0.0762804 0.0985859 -0.12736 -0.0345689 -0.0430195 \n", " 0.0397482 -0.118882 0.126059 -0.031277 -0.0195955 \n", " ⋮ ⋱ \n", " 0.0174594 0.161578 0.00648771 0.024462 0.0263169 \n", " -0.151761 -0.0428635 -0.0675281 -0.0366188 0.0180868 \n", " 0.0441573 -0.1596 0.0519537 … 0.0467241 -0.0847827 \n", " 0.167689 -0.099587 0.00711909 0.00442244 -0.192641 \n", " 0.0238481 -0.0608836 0.00945763 0.120106 0.163596 \n", " 0.231379 -0.179251 0.0710683 -0.016659 0.185589 \n", " -0.0081851 -0.114456 -0.0345721 0.102517 0.11927 \n", " 0.13378 -0.228392 -0.0835398 … 0.0431003 -0.0204308 \n", " 0.121963 -0.211102 0.0592341 0.0417325 -0.0403561 \n", " -0.122187 0.121933 0.0543548 -0.104907 -0.019856 \n", " 0.0256914 -0.0222105 -0.0712324 0.0778704 0.0626852 \n", " 0.0799711 -0.100243 0.0141308 -0.04014 -0.144645 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "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": [ { "data": { "text/plain": [ "17" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "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 = randperm(N)\n", "sel = sel[1 : 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*xsharp;" ] }, { "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": [ { "data": { "text/plain": [ "(::#1) (generic function with 1 method)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prox_gamma_g = (x, gamma) -> x - x./max(abs(x)./gamma, 1)" ] }, { "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": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGgCAYAAACwio2MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt4VPWB//HPJCETvGSy3BJSIogXbqIol9lQVqmMhcVb1T6AsnJZhEWFKqGUpF2h7HafiGUVCyhVC9iKQmkBrdooclmEjSABqoWACioRmCDNZgKKATLf3x/+mDrmwkwyZ2bOzPv1POd5OiffM/P9egjz7syZwWGMMQIAALChlFhPAAAAoLkIGQAAYFuEDAAAsC1CBgAA2BYhAwAAbIuQAQAAtkXIAAAA2yJkAACAbREyAADAtggZAABgW4QMAACwrbRYTyDS/H6/jhw5oosvvlgOhyPW0wEAACEwxujEiRPKzc1VSkoYr7OYKFi4cKHp3LmzcTqdZsCAAWbbtm1Njn/hhRfM1VdfbVq3bm1ycnLM+PHjzfHjx0N6rIqKCiOJjY2NjY2NzYZbRUVFWI3hMMbaf/165cqVGjNmjBYvXiy326358+dr1apV2r9/vzp06FBv/NatW3X99dfriSee0K233qrDhw9r8uTJuvLKK7V69erzPp7P51NWVpYqKiqUmZlpxZIAAECE1dTUKC8vT9XV1XK5XCEfZ3nIuN1u9e/fXwsXLpT09Vs/eXl5mjp1qgoLC+uNnzdvnp5++mkdOHAgsG/BggWaO3euPvvss/M+Xk1NjVwul3w+HyEDAIBNNPf529KLfU+fPq2ysjJ5PJ6/P2BKijwej0pLSxs8Jj8/XxUVFXr99ddljFFlZaVWrVql4cOHNzi+trZWNTU1QRsAAEgOlobM8ePHVVdXp+zs7KD92dnZ8nq9DR7z3e9+V8uXL9fIkSOVnp6unJwcZWVladGiRQ2OLy4ulsvlCmx5eXkRXwcAAIhPcffx67179+qhhx7SrFmzVFZWppKSEn3yySeaPHlyg+OLiork8/kCW0VFRZRnDAAAYsXSj1+3a9dOqampqqysDNpfWVmpnJycBo8pLi7WwIEDNWPGDEnS1VdfrQsvvFD/9E//pF/84hfq2LFj0Hin0ymn02nNAgAAQFyz9BWZ9PR09e3bV+vXrw/s8/v9Wr9+vfLz8xs85ssvv1RaWnBfpaamSpIsvi4ZAADYjOVvLRUUFOjZZ5/V888/r/Lyct1///364osvNH78eElfvzU0ZsyYwPhbb71Vf/zjH/X000/r4MGD2rp1q370ox9pwIABys3NtXq6AADARiz/Zt+RI0fq888/16xZs+T1etWnTx+VlJQELgA+evSoDh06FBg/btw4nThxQgsXLtT06dOVlZWlG2+8UXPnzrV6qgAAwGYs/x6ZaON7ZAAAsJ+4/B4ZAAAAKxEyAADAtggZAABgW4QMAACwLUIGAADYFiEDAABsi5ABAAC2RcgAAADbImQAAIBtETIAAMC2CBkAAGBbhAwAALAtQgYAANgWIQMAAGyLkAEAACFbsf2QNn/weaynEUDIAACAkLy0/ZAKV7+v+367Qx8dOxHr6UgiZAAAQAhe2n5IRavflyT9i7uzLmt/UYxn9DVCBgAANOmbEfOv371Uj9zSQw6HI8az+hohAwAAGhXPESMRMgAAoBHxHjESIQMAABpgh4iRCBkAAPAtdokYiZABAADfYKeIkQgZAADw/9ktYiRCBgAAyJ4RIxEyAAAkPbtGjETIAACQ1OwcMRIhAwBA0rJ7xEiEDAAASSkRIkYiZAAASDqJEjESIQMAQFJJpIiRCBkAAJJGokWMRMgAAJAUEjFiJEIGAICEl6gRIxEyAAAktESOGClKIbNo0SJ16dJFGRkZcrvd2r59e5Pja2tr9bOf/UydO3eW0+lUly5dtGTJkmhMFQCAhJHoESNJaVY/wMqVK1VQUKDFixfL7XZr/vz5Gjp0qPbv368OHTo0eMyIESNUWVmp3/zmN7r88st19OhR+f1+q6cKAEDCSIaIkSSHMcZY+QBut1v9+/fXwoULJUl+v195eXmaOnWqCgsL640vKSnRqFGjdPDgQbVp0+a8919bW6va2trA7ZqaGuXl5cnn8ykzMzNyCwEAwCbsGDE1NTVyuVxhP39b+tbS6dOnVVZWJo/H8/cHTEmRx+NRaWlpg8e88sor6tevnx577DF95zvf0ZVXXqkf//jHOnXqVIPji4uL5XK5AlteXp4lawEAwA7sGDEtYWnIHD9+XHV1dcrOzg7an52dLa/X2+AxBw8e1JYtW/TXv/5Va9as0fz58/WHP/xBDzzwQIPji4qK5PP5AltFRUXE1wEAgB0kW8RIUbhGJlx+v18Oh0PLly+Xy+WSJD3++OP64Q9/qKeeekqtW7cOGu90OuV0OmMxVQAA4kYyRoxk8Ssy7dq1U2pqqiorK4P2V1ZWKicnp8FjOnbsqO985zuBiJGkHj16yBijzz77zMrpAgBgS8kaMZLFIZOenq6+fftq/fr1gX1+v1/r169Xfn5+g8d897vf1ZEjR3Ty5MnAvg8++EApKSnq1KmTldMFAMB2kjlipCh8j0xBQYGeffZZPf/88yovL9f999+vL774QuPHj5f09TUuY8aMCYy/55571LZtW40fP1579+7V5s2bNWPGDP3rv/5rvbeVAABIZskeMVIUrpEZOXKkPv/8c82aNUter1d9+vRRSUlJ4ALgo0eP6tChQ4HxF110kdatW6epU6eqX79+atu2rUaMGKFf/OIXVk8VAADbIGK+Zvn3yERbcz+HDgCAXSRixMTl98gAAIDISsSIaQlCBgAAmyBi6iNkAACwASKmYYQMAABxjohpHCEDAEAcI2KaRsgAABCniJjzI2QAAIhDRExoCBkAAOIMERM6QgYAgDhCxISHkAEAIE4QMeEjZAAAiANETPMQMgAAxBgR03yEDAAAMUTEtAwhAwBAjBAxLUfIAAAQA0RMZBAyAABEGRETOYQMAABRRMREFiEDAECUEDGRR8gAABAFRIw1CBkAACxGxFiHkAEAwEJEjLUIGQAALELEWI+QAQDAAkRMdBAyAABEGBETPYQMAAARRMREFyEDAECEEDHRR8gAABABRExsEDIAALQQERM7hAwAAC1AxMQWIQMAQDMRMbFHyAAA0AxETHwgZAAACBMREz8IGQAAwkDExBdCBgCAEBEx8ScqIbNo0SJ16dJFGRkZcrvd2r59e0jHbd26VWlpaerTp4/FMwQAoGlETHyyPGRWrlypgoICzZ49Wzt37tQ111yjoUOH6tixY00eV11drTFjxmjIkCFWTxEAgCYRMfHLYYwxVj6A2+1W//79tXDhQkmS3+9XXl6epk6dqsLCwkaPGzVqlK644gqlpqZq7dq12r17d4PjamtrVVtbG7hdU1OjvLw8+Xw+ZWZmRnYxAICkQ8RER01NjVwuV9jP35a+InP69GmVlZXJ4/H8/QFTUuTxeFRaWtrocUuXLtXBgwc1e/bs8z5GcXGxXC5XYMvLy4vI3AEAIGLin6Uhc/z4cdXV1Sk7Oztof3Z2trxeb4PHfPjhhyosLNQLL7ygtLS08z5GUVGRfD5fYKuoqIjI3AEAyY2IsYfzl0IU1dXV6Z577tGcOXN05ZVXhnSM0+mU0+m0eGYAgGRCxNiHpSHTrl07paamqrKyMmh/ZWWlcnJy6o0/ceKEduzYoV27dmnKlCmSvr6mxhijtLQ0vfnmm7rxxhutnDIAIMkRMfZi6VtL6enp6tu3r9avXx/Y5/f7tX79euXn59cbn5mZqffff1+7d+8ObJMnT1a3bt20e/duud1uK6cLAEhyRIz9WP7WUkFBgcaOHat+/fppwIABmj9/vr744guNHz9e0tfXuBw+fFi//e1vlZKSoquuuiro+A4dOigjI6PefgAAIomIsSfLQ2bkyJH6/PPPNWvWLHm9XvXp00clJSWBC4CPHj2qQ4cOWT0NAAAaRcTYl+XfIxNtzf0cOgAgOREx8SEuv0cGAIB4RsTYHyEDAEhKRExiIGQAAEmHiEkchAwAIKkQMYmFkAEAJA0iJvEQMgCApEDEJCZCBgCQ8IiYxEXIAAASGhGT2AgZAEDCImISHyEDAEhIRExyIGQAAAmHiEkehAwAIKEQMcmFkAEAJAwiJvkQMgCAhEDEJCdCBgBge0RM8iJkAAC2RsQkN0IGAGBbRAwIGQCALRExkAgZAIANETE4h5ABANgKEYNvImQAALZBxODbCBkAgC0QMWgIIQMAiHtEDBpDyAAA4hoRg6YQMgCAuEXE4HwIGQBAXCJiEApCBgAQd4gYhIqQAQDEFSIG4SBkAABxg4hBuAgZAEBcIGLQHIQMACDmiBg0FyEDAIgpIgYtQcgAAGKGiEFLETIAgJggYhAJUQmZRYsWqUuXLsrIyJDb7db27dsbHbt69WrddNNNat++vTIzM5Wfn6833ngjGtMEAEQJEYNIsTxkVq5cqYKCAs2ePVs7d+7UNddco6FDh+rYsWMNjt+8ebNuuukmvf766yorK9P3vvc93Xrrrdq1a5fVUwUARAERg0hyGGOMlQ/gdrvVv39/LVy4UJLk9/uVl5enqVOnqrCwMKT76NWrl0aOHKlZs2bV+1ltba1qa2sDt2tqapSXlyefz6fMzMzILAIAEBFEDBpTU1Mjl8sV9vO3pa/InD59WmVlZfJ4PH9/wJQUeTwelZaWhnQffr9fJ06cUJs2bRr8eXFxsVwuV2DLy8uLyNwBAJFFxMAKlobM8ePHVVdXp+zs7KD92dnZ8nq9Id3HvHnzdPLkSY0YMaLBnxcVFcnn8wW2ioqKFs8bABBZRAyskhbrCTTlxRdf1Jw5c/Tyyy+rQ4cODY5xOp1yOp1RnhkAIFREDKxkaci0a9dOqampqqysDNpfWVmpnJycJo9dsWKF7rvvPq1atSrorSkAgH0QMbCapW8tpaenq2/fvlq/fn1gn9/v1/r165Wfn9/ocS+99JLGjx+vl156STfffLOVUwQAWISIQTRY/tZSQUGBxo4dq379+mnAgAGaP3++vvjiC40fP17S19e4HD58WL/97W8lff120tixY/Xkk0/K7XYHrqVp3bq1XC6X1dMFAEQAEYNosTxkRo4cqc8//1yzZs2S1+tVnz59VFJSErgA+OjRozp06FBg/DPPPKOzZ8/qwQcf1IMPPhjYP3bsWC1btszq6QIAWoiIQTRZ/j0y0dbcz6EDAFqOiEFzxeX3yAAAkgcRg1ggZAAALUbEIFYIGQBAixAxiCVCBgDQbEQMYo2QAQA0CxGDeEDIAADCRsQgXhAyAICwEDGIJ4QMACBkRAziDSEDAAgJEYN4RMgAAM6LiEG8ImQAAE0iYhDPCBkAQKOIGMQ7QgYA0CAiBnZAyAAA6iFiYBeEDAAgCBEDOyFkAAABRAzshpABAEgiYmBPhAwAgIiBbREyAJDkiBjYGSEDAEmMiIHdETIAkKSIGCQCQgYAkhARg0RByABAkiFikEgIGQBIIkQMEg0hAwBJgohBIiJkACAJEDFIVIQMACQ4IgaJjJABgARGxCDRETIAkKCIGCQDQgYAEhARg2RByABAgiFikEwIGQBIIEQMkg0hAwAJgohBMiJkACABEDFIVlEJmUWLFqlLly7KyMiQ2+3W9u3bmxy/adMmXXfddXI6nbr88su1bNmyaEwTAGyJiEEyszxkVq5cqYKCAs2ePVs7d+7UNddco6FDh+rYsWMNjv/44491880363vf+552796thx9+WPfdd5/eeOMNq6cKALZDxCDZOYwxxsoHcLvd6t+/vxYuXChJ8vv9ysvL09SpU1VYWFhv/MyZM/Xaa6/pr3/9a2DfqFGjVF1drZKSknrja2trVVtbG7hdU1OjvLw8+Xw+ZWZmWrAiIHrO1vn14vZD+vj4F7GeCuLQya/OalXZZ5KIGNhfTU2NXC5X2M/faRbOSadPn1ZZWZmKiooC+1JSUuTxeFRaWtrgMaWlpfJ4PEH7hg4dqocffrjB8cXFxZozZ07kJg3EibN1fhX8/i965S9HYj0VxDkiBsnM0pA5fvy46urqlJ2dHbQ/Oztb+/bta/AYr9fb4PiamhqdOnVKrVu3DvpZUVGRCgoKArfPvSID2Nk3IyYtxaEx+V3UOp1r81HfldkX67ZrcokYJC1LQyYanE6nnE5nrKcBRMy3I+ap0dfp+71yYj0tAIhLlv5fvHbt2ik1NVWVlZVB+ysrK5WT0/BfzDk5OQ2Oz8zMrPdqDJBoiBgACI+lIZOenq6+fftq/fr1gX1+v1/r169Xfn5+g8fk5+cHjZekdevWNToeSBREDACEz/I33QsKCvTss8/q+eefV3l5ue6//3598cUXGj9+vKSvr3EZM2ZMYPzkyZN18OBB/eQnP9G+ffv01FNP6fe//72mTZtm9VSBmCFiAKB5LL9GZuTIkfr88881a9Yseb1e9enTRyUlJYELeo8ePapDhw4Fxl966aV67bXXNG3aND355JPq1KmTnnvuOQ0dOtTqqQIxQcQAQPNZ/j0y0dbcz6EDsUDEAMDXmvv8zec5gRghYgCg5QgZIAaIGACIDEIGiDIiBgAih5ABooiIAYDIImSAKCFiACDyCBkgCogYALAGIQNYjIgBAOsQMoCFiBgAsBYhA1iEiAEA6xEygAWIGACIDkIGiDAiBgCih5ABIoiIAYDoImSACCFiACD6CBkgAogYAIgNQgZoISIGAGKHkAFagIgBgNgiZIBmImIAIPYIGaAZiBgAiA+EDBAmIgYA4gchA4SBiAGA+ELIACEiYgAg/hAyQAiIGACIT4QMcB5EDADEL0IGaAIRAwDxjZABGkHEAED8I2SABhAxAGAPhAzwLUQMANgHIQN8AxEDAPZCyAD/HxEDAPZDyAAiYgDArggZJD0iBgDsi5BBUiNiAMDeCBkkLSIGAOyPkEFSImIAIDFYFjJVVVUaPXq0MjMzlZWVpQkTJujkyZONjj9z5oxmzpyp3r1768ILL1Rubq7GjBmjI0eOWDVFJCkiBgASh2UhM3r0aO3Zs0fr1q3Tq6++qs2bN2vSpEmNjv/yyy+1c+dOPfLII9q5c6dWr16t/fv367bbbrNqikhCRAwAJBaHMcZE+k7Ly8vVs2dPvfvuu+rXr58kqaSkRMOHD9dnn32m3NzckO7n3Xff1YABA/Tpp5/qkksuCemYmpoauVwu+Xw+ZWZmNnsNSDxEDADEr+Y+f1vyikxpaamysrICESNJHo9HKSkp2rZtW8j34/P55HA4lJWV1eiY2tpa1dTUBG3AtxExAJCYLAkZr9erDh06BO1LS0tTmzZt5PV6Q7qPr776SjNnztTdd9/dZJkVFxfL5XIFtry8vBbNHYmHiAGAxBVWyBQWFsrhcDS57du3r8WTOnPmjEaMGCFjjJ5++ukmxxYVFcnn8wW2ioqKFj8+EgcRAwCJLS2cwdOnT9e4ceOaHNO1a1fl5OTo2LFjQfvPnj2rqqoq5eQ0/SRyLmI+/fRTbdiw4bzvkzmdTjmdzpDmj+RCxABA4gsrZNq3b6/27dufd1x+fr6qq6tVVlamvn37SpI2bNggv98vt9vd6HHnIubDDz/Uxo0b1bZt23CmBwQQMQCQHCy5RqZHjx4aNmyYJk6cqO3bt2vr1q2aMmWKRo0aFfSJpe7du2vNmjWSvo6YH/7wh9qxY4eWL1+uuro6eb1eeb1enT592oppIkERMQCQPMJ6RSYcy5cv15QpUzRkyBClpKTorrvu0q9+9augMfv375fP55MkHT58WK+88ookqU+fPkHjNm7cqMGDB1s1VSQQIgYAkosl3yMTS3yPTPIiYgDAvuLqe2SAaCNiACA5ETKwPSIGAJIXIQNbI2IAILkRMrAtIgYAQMjAlogYAIBEyMCGiBgAwDmEDGyFiAEAfBMhA9sgYgAA30bIwBaIGABAQwgZxD0iBgDQGEIGcY2IAQA0hZBB3CJiAADnQ8ggLhExAIBQEDKIO0QMACBUhAziChEDAAgHIYO4QcQAAMJFyCAuEDEAgOYgZBBzRAwAoLkIGcQUEQMAaAlCBjFDxAAAWoqQQUwQMQCASCBkEHVEDAAgUggZRBURAwCIJEIGUUPEAAAijZBBVBAxAAArEDKwHBEDALAKIQNLETEAACsRMrAMEQMAsBohA0sQMQCAaCBkEHFEDAAgWggZRBQRAwCIJkIGEUPEAACijZBBRBAxAIBYIGTQYkQMACBWLAuZqqoqjR49WpmZmcrKytKECRN08uTJkI+fPHmyHA6H5s+fb9UUEQFEDAAgliwLmdGjR2vPnj1at26dXn31VW3evFmTJk0K6dg1a9bonXfeUW5urlXTQwQQMQCAWLMkZMrLy1VSUqLnnntObrdbgwYN0oIFC7RixQodOXKkyWMPHz6sqVOnavny5WrVqpUV00MEEDEAgHhgSciUlpYqKytL/fr1C+zzeDxKSUnRtm3bGj3O7/fr3nvv1YwZM9SrV6+QHqu2tlY1NTVBG6xFxAAA4oUlIeP1etWhQ4egfWlpaWrTpo28Xm+jx82dO1dpaWn60Y9+FPJjFRcXy+VyBba8vLxmzxvnR8QAAOJJWCFTWFgoh8PR5LZv375mTaSsrExPPvmkli1bJofDEfJxRUVF8vl8ga2ioqJZj4/zI2IAAPEmLZzB06dP17hx45oc07VrV+Xk5OjYsWNB+8+ePauqqirl5DT8xPf222/r2LFjuuSSSwL76urqNH36dM2fP1+ffPJJg8c5nU45nc5wloFmIGIAAPEorJBp37692rdvf95x+fn5qq6uVllZmfr27StJ2rBhg/x+v9xud4PH3HvvvfJ4PEH7hg4dqnvvvVfjx48PZ5qIMCIGABCvwgqZUPXo0UPDhg3TxIkTtXjxYp05c0ZTpkzRqFGjgj5S3b17dxUXF+uOO+5Q27Zt1bZt26D7adWqlXJyctStWzcrpokQEDEAgHhm2ffILF++XN27d9eQIUM0fPhwDRo0SM8880zQmP3798vn81k1BbQQEQMAiHcOY4yJ9SQiqaamRi6XSz6fT5mZmbGejm0RMQCAaGru8zf/1hLqIWIAAHZByCAIEQMAsBNCBgFEDADAbggZSCJiAAD2RMiAiAEA2BYhk+SIGACAnREySYyIAQDYHSGTpIgYAEAiIGSSEBEDAEgUhEySIWIAAImEkEkiRAwAINEQMkmCiAEAJCJCJgkQMQCAREXIJDgiBgCQyAiZBEbEAAASHSGToIgYAEAyIGQSEBEDAEgWhEyCIWIAAMmEkEkgRAwAINkQMgmCiAEAJCNCJgEQMQCAZEXI2BwRAwBIZoSMjRExAIBkR8jYFBEDAAAhY0tEDAAAXyNkbIaIAQDg7wgZGyFiAAAIRsjYBBEDAEB9hIwNEDEAADSMkIlzRAwAAI0jZOIYEQMAQNMImThFxAAAcH6ETBwiYgAACA0hE2eIGAAAQmdZyFRVVWn06NHKzMxUVlaWJkyYoJMnT573uPLyct12221yuVy68MIL1b9/fx06dMiqacYVIgYAgPBYFjKjR4/Wnj17tG7dOr366qvavHmzJk2a1OQxBw4c0KBBg9S9e3dt2rRJ7733nh555BFlZGRYNc24QcQAABA+hzHGRPpOy8vL1bNnT7377rvq16+fJKmkpETDhw/XZ599ptzc3AaPGzVqlFq1aqXf/e53zX7smpoauVwu+Xw+ZWZmNvt+oomIAQAku+Y+f1vyikxpaamysrICESNJHo9HKSkp2rZtW4PH+P1+vfbaa7ryyis1dOhQdejQQW63W2vXrm3ysWpra1VTUxO02QkRAwBA81kSMl6vVx06dAjal5aWpjZt2sjr9TZ4zLFjx3Ty5Ek9+uijGjZsmN58803dcccduvPOO/U///M/jT5WcXGxXC5XYMvLy4voWqxExAAA0DJhhUxhYaEcDkeT2759+5o1Eb/fL0m6/fbbNW3aNPXp00eFhYW65ZZbtHjx4kaPKyoqks/nC2wVFRXNevxoI2IAAGi5tHAGT58+XePGjWtyTNeuXZWTk6Njx44F7T979qyqqqqUk9Pwk3W7du2Ulpamnj17Bu3v0aOHtmzZ0ujjOZ1OOZ3O0BYQJ4gYAAAiI6yQad++vdq3b3/ecfn5+aqurlZZWZn69u0rSdqwYYP8fr/cbneDx6Snp6t///7av39/0P4PPvhAnTt3DmeacY2IAQAgciy5RqZHjx4aNmyYJk6cqO3bt2vr1q2aMmWKRo0aFfSJpe7du2vNmjWB2zNmzNDKlSv17LPP6qOPPtLChQv1pz/9SQ888IAV04w6IgYAgMiy7Htkli9fru7du2vIkCEaPny4Bg0apGeeeSZozP79++Xz+QK377jjDi1evFiPPfaYevfureeee05//OMfNWjQIKumGTVEDAAAkWfJ98jEUjx+jwwRAwBA0+Lqe2Twd0QMAADWIWQsRMQAAGAtQsYiRAwAANYjZCxAxAAAEB2ETIQRMQAARA8hE0FEDAAA0UXIRAgRAwBA9BEyEUDEAAAQG4RMCxExAADEDiHTAkQMAACxRcg0ExEDAEDsETLNQMQAABAfCJkwETEAAMQPQiYMRAwAAPGFkAkREQMAQPwhZEJUfvSESv7qJWIAAIgjabGegF307uTSr8f01emzfiIGAIA4QciE4XvdOsR6CgAA4Bt4awkAANgWIQMAAGyLkAEAALZFyAAAANsiZAAAgG0RMgAAwLYIGQAAYFuEDAAAsC1CBgAA2BYhAwAAbIuQAQAAtkXIAAAA2yJkAACAbREyAADAtggZAABgW4QMAACwLctCpqqqSqNHj1ZmZqaysrI0YcIEnTx5ssljTp48qSlTpqhTp05q3bq1evbsqcWLF1s1RQAAYHOWhczo0aO1Z88erVu3Tq+++qo2b96sSZMmNXlMQUGBSkpK9MILL6i8vFzTpk3TlClT9Morr1g1TQAAYGOWhEx5eblKSkr03HPPye12a9CgQVqwYIFWrFihI0eONHrc//7v/2rs2LEaPHiwunTpookTJ+qaa67R9u3brZgmAACwOUtCprS0VFlZWerXr1/D1dy9AAAPE0lEQVRgn8fjUUpKirZt29bocQMHDtQrr7yiw4cPyxijjRs36oMPPtD3v//9Ro+pra1VTU1N0AYAAJKDJSHj9XrVoUOHoH1paWlq06aNvF5vo8ctWLBAPXv2VKdOnZSenq5hw4Zp0aJFuv766xs9pri4WC6XK7Dl5eVFbB0AACC+hRUyhYWFcjgcTW779u1r9mQWLFigd955R6+88orKysr03//933rwwQf11ltvNXpMUVGRfD5fYKuoqGj24wMAAHtJC2fw9OnTNW7cuCbHdO3aVTk5OTp27FjQ/rNnz6qqqko5OTkNHnfq1Cn99Kc/1erVq3XLLbdIkq6++mrt3r1b8+bNk8fjafA4p9Mpp9MZzjIAAECCCCtk2rdvr/bt2593XH5+vqqrq1VWVqa+fftKkjZs2CC/3y+3293gMWfOnNGZM2eUlhY8pdTUVPn9/pDnaIyRJK6VAQDARs49b597Hg+ZsciwYcPMtddea7Zt22a2bNlirrjiCnP33XcHjenWrZtZvXp14PYNN9xgevXqZTZu3GgOHjxoli5dajIyMsxTTz0V8uNWVFQYSWxsbGxsbGw23CoqKsLqDYcx4aZPaKqqqjRlyhT96U9/UkpKiu666y796le/0kUXXRQY43A4tHTp0sDbVV6vV0VFRXrzzTdVVVWlzp07a9KkSZo2bZocDkdIj+v3+3XkyBFdfPHFIR8TqpqaGuXl5amiokKZmZkRve94kOjrkxJ/jazP/hJ9jYm+Pinx12jV+owxOnHihHJzc5WSEvolvJaFTCKqqamRy+WSz+dL2D+cibw+KfHXyPrsL9HXmOjrkxJ/jfG2Pv6tJQAAYFuEDAAAsK3Un//85z+P9STsJDU1VYMHD6736apEkejrkxJ/jazP/hJ9jYm+Pinx1xhP6+MaGQAAYFu8tQQAAGyLkAEAALZFyAAAANsiZAAAgG0RMgAAwLYImW/4r//6Lw0cOFAXXHCBsrKyQjrGGKNZs2apY8eOat26tTwejz788MOgMV999ZUefPBBtW3bVhdddJHuuusuVVZWWrGE86qqqtLo0aOVmZmprKwsTZgwQSdPnmzyGIfD0eD2y1/+MjBm8ODB9X4+efJkq5dTT3PWN27cuHpzHzZsWNCYeDmH4a7vzJkzmjlzpnr37q0LL7xQubm5GjNmjI4cORI0Lpbnb9GiRerSpYsyMjLkdru1ffv2Jsdv2rRJ1113nZxOpy6//HItW7as3phVq1ape/fuysjIUO/evfX6669bNPvzC2d9q1ev1k033aT27dsrMzNT+fn5euONN4LGLFu2rN65ysjIsHoZTQpnjZs2bWrw7xOv1xs0zq7nsKG/TxwOh3r16hUYE0/ncPPmzbr11luVm5srh8OhtWvXnveYuPsdDOtfZkpws2bNMo8//rgpKCgwLpcrpGMeffRR43K5zNq1a81f/vIXc9ttt5lLL73UnDp1KjBm8uTJJi8vz6xfv97s2LHD/OM//qMZOHCgVcto0rBhw8w111xj3nnnHfP222+byy+/vN4/5vltR48eDdqWLFliHA6HOXDgQGDMDTfcYCZOnBg0zufzWb2cepqzvrFjx5phw4YFzb2qqipoTLycw3DXV11dbTwej1m5cqXZt2+fKS0tNQMGDDB9+/YNGher87dixQqTnp5ulixZYvbs2WMmTpxosrKyTGVlZYPjDx48aC644AJTUFBg9u7daxYsWGBSU1NNSUlJYMzWrVtNamqqeeyxx8zevXvNv//7v5tWrVqZ999/3/L1fFu463vooYfM3Llzzfbt280HH3xgioqKTKtWrczOnTsDY5YuXWoyMzODzpXX643WkuoJd40bN240ksz+/fuD1lBXVxcYY+dzWF1dHbSuiooK06ZNGzN79uzAmHg6h6+//rr52c9+ZlavXm0kmTVr1jQ5Ph5/BwmZBixdujSkkPH7/SYnJ8f88pe/DOyrrq42TqfTvPTSS4HbrVq1MqtWrQqMKS8vN5JMaWlp5CffhL179xpJ5t133w3s+/Of/2wcDoc5fPhwyPdz++23mxtvvDFo3w033GAeeuihiM21OZq7vrFjx5rbb7+90Z/HyzmM1Pnbvn27kWQ+/fTTwL5Ynb8BAwaYBx98MHC7rq7O5ObmmuLi4gbH/+QnPzG9evUK2jdy5EgzdOjQwO0RI0aYm2++OWiM2+02//Zv/xbBmYcm3PU1pGfPnmbOnDmB26H+/RQt4a7xXMj83//9X6P3mUjncM2aNcbhcJhPPvkksC/ezuE5oYRMPP4O8tZSC3z88cfyer3yeDyBfS6XS263W6WlpZKksrIynTlzJmhM9+7ddckllwTGREtpaamysrLUr1+/wD6Px6OUlBRt27YtpPuorKzUa6+9pgkTJtT72fLly9WuXTtdddVVKioq0pdffhmxuYeiJevbtGmTOnTooG7duun+++/X3/72t8DP4uUcRuL8SZLP55PD4aj39mm0z9/p06dVVlYW9N81JSVFHo+n0f+upaWlQeMlaejQoUHjQxkTDc1Z37f5/X6dOHFCbdq0Cdp/8uRJde7cWXl5ebr99tu1Z8+eiM49VC1ZY58+fdSxY0fddNNN2rp1a9DPEukc/uY3v5HH41Hnzp2D9sfLOQxXPP4Oxv67hW3s3Hu62dnZQfuzs7MDP/N6vUpPT6/3pPHNMdHi9XrVoUOHoH1paWlq06ZNyHN5/vnndfHFF+vOO+8M2n/PPfeoc+fOys3N1XvvvaeZM2dq//79Wr16dcTmfz7NXd+wYcN055136tJLL9WBAwf005/+VP/8z/+s0tJSpaamxs05jMT5++qrrzRz5kzdfffdQf9qbSzO3/Hjx1VXV9fg78++ffsaPMbr9TY4vqamRqdOnVLr1q0bHRPt37fmrO/b5s2bp5MnT2rEiBGBfd26ddOSJUt09dVXy+fzad68eRo4cKD27NmjTp06RXQN59OcNXbs2FGLFy9Wv379VFtbq+eee06DBw/Wtm3bdN1110lq/Dzb7RweOXJEf/7zn/Xiiy8G7Y+ncxiuePwdTPiQKSws1Ny5c5scU15eru7du0dpRpEX6hojYcmSJRo9enS9C9MmTZoU+N+9e/dWbm6ubrzxRh04cECXXXZZix7T6vWNGjUq8L979+6tq6++Wpdddpk2bdqkIUOGNPt+QxWt83fmzBmNGDFCxhg9/fTTQT+z8vyheV588UXNmTNHL7/8clDA5ufnKz8/P3B74MCB6tGjh37961/rP//zP2Mx1bB069ZN3bp1C9weOHCgDhw4oCeeeEK/+93vYjizyHv++eeVlZWlH/zgB0H77X4O403Ch8z06dM1bty4Jsd07dq1Wfedk5Mj6eu3Wzp27BjYX1lZqT59+gTGnD59WtXV1UH/j76ysjJwfEuFusacnBwdO3YsaP/Zs2dVVVUV0lzefvtt7d+/XytXrjzv2AEDBkiSPvrooxY/EUZrfd+8r3bt2umjjz7SkCFDLD+H0VjfuYj59NNPtWHDhqBXYxoSyfPXmHbt2ik1NbXep7+a+u+ak5PT4PjMzEy1bt26yTGR+n0LVXPWd86KFSt03333adWqVfVeov+2Vq1a6dprr9VHH33U4jmHqyVr/KYBAwZoy5YtgduJcA6NMVqyZInuvfdepaenNzk2lucwXHH5O2jJlTc2F+7FvvPmzQvs8/l8DV7s+4c//CEwZt++fTG92HfHjh2BfW+88UbIF4uOHTu23qddGrNlyxYjyfzlL39p9nzD1dL1nVNRUWEcDod5+eWXjTHxcw6bu77Tp0+bH/zgB6ZXr17m2LFjIT1WtM7fgAEDzJQpUwK36+rqzHe+850mL/a96qqrgvbdfffd9S40vOWWW4LG5Ofnx+xC0XDWZ4wxL774osnIyDBr164N6THOnj1runXrZqZNm9bi+TZHc9b4bR6Px9xxxx2B23Y/h8b8/aLmUD6pE+tzeI5CvNg33n4HCZlv+PTTT82uXbvMnDlzzEUXXWR27dpldu3aZU6cOBEY061bN7N69erA7UcffdRkZWWZl19+2bz33nvm9ttvb/Dj15dcconZsGGD2bFjh8nPzzf5+flRXds5w4YNM9dee63Ztm2b2bJli7niiivqfXz322s05utAu+CCC8zTTz9d7z4/+ugj8x//8R9mx44d5uOPPzYvv/yy6dq1q7n++ustXUtDwl3fiRMnzI9//GNTWlpqPv74Y/PWW2+Z6667zlxxxRXmq6++ChwTL+cw3PWdPn3a3HbbbaZTp05m9+7dQR/3rK2tNcbE9vytWLHCOJ1Os2zZMrN3714zadIkk5WVFfgoamFhobn33nsD48999HPGjBmmvLzcLFq0qMGPfqalpZl58+aZ8vJyM3v27Jh+dDec9S1fvtykpaWZRYsWBZ2r6urqwJg5c+aYN954wxw4cMCUlZWZUaNGmYyMDLNnz56or8+Y8Nf4xBNPmLVr15oPP/zQvP/+++ahhx4yKSkp5q233gqMsfM5POdf/uVfjNvtbvA+4+kcnjhxIvBcJ8k8/vjjZteuXYFPNdrhd5CQ+YaxY8caSfW2jRs3BsZIMkuXLg3c9vv95pFHHjHZ2dnG6XSaIUOGmP379wfd76lTp8wDDzxg/uEf/sFccMEF5o477jBHjx6N0qqC/e1vfzN33323ueiii0xmZqYZP358UKgZU3+Nxhjz61//2rRu3TroL9RzDh06ZK6//nrTpk0b43Q6zeWXX25mzJgRk++RCXd9X375pfn+979v2rdvb1q1amU6d+5sJk6cWO87HeLlHIa7vo8//rjBP9Pf/HMd6/O3YMECc8kll5j09HQzYMAA88477wR+NnbsWHPDDTcEjd+4caPp06ePSU9PN127dq33Z9UYY37/+9+bK6+80qSnp5tevXqZ1157zeJVNC6c9d1www0NnquxY8cGxjz88MOB+8vOzjbDhw8P+p6ZWAhnjXPnzjWXXXaZycjIMG3atDGDBw82GzZsqHefdj2Hxnz9Km7r1q3NM8880+D9xdM5PPfKUWN/5uzwO+gwxhhr3rQCAACwFt8jAwAAbIuQAQAAtkXIAAAA2yJkAACAbREyAADAtggZAABgW4QMAACwLUIGAADYFiEDAABsi5ABAAC2RcgAAADb+n85ea8M8FWwcQAAAABJRU5ErkJggg==", "text/plain": [ "PyPlot.Figure(PyObject )" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "(-1.1,1.1,-0.7699999999999999,0.7699999999999999)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "figsize = (9, 6)\n", "t = -1 : 0.001 : 1\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": [ { "data": { "text/plain": [ "(::#3) (generic function with 1 method)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pA = pinv(A) # pseudo-inverse. Equivalent to pA = A.T.dot(inv(A.dot(A.T)))\n", "prox_f = (x, y) -> x + pA*(y - A*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": [ { "data": { "text/plain": [ "1" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "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": [ { "data": { "text/plain": [ "700" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nbiter = 700" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

\n", " \n", "__Exercise: 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": [ { "ename": "LoadError", "evalue": "syntax: incomplete: \"for\" at In[14]:3 requires end", "output_type": "error", "traceback": [ "syntax: incomplete: \"for\" at In[14]:3 requires end", "" ] } ], "source": [ "s = zeros(N)\n", "En_array = zeros(nbiter)\n", "for iter in 1 : nbiter # iter goes from 1 to nbiter\n", " # put your code here\n", " En_array[iter] = maximum(sum(abs(x), 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": [ { "ename": "LoadError", "evalue": "UndefVarError: x_restored not defined", "output_type": "error", "traceback": [ "UndefVarError: x_restored not defined", "" ] } ], "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](L\"$x^\\sharp$\")\n", "subfig2[:set_title](L\"$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": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGgCAYAAABMn6ZGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAH7RJREFUeJzt3X9s1dXh//HXbWtvO8e9hQD3Url8KMJEFOkAqSUy57yxKCpsLANkKggyycDxw0FZFIguARnbAKkjKbq6xA1xGYxYV4LFxh+7FgW7gfyISo1VuBdZ03tLwQLt+f7hlzuvLbWwvvvj9PlI7h89Pe97zzlh9pn7ay5jjBEAAIClkjp6AQAAAE4idgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWI3YAAIDViB0AAGC1lI5eQEdobGzUsWPH1KNHD7lcro5eDgAAaAVjjGpra5WZmamkpNY/X9MtY+fYsWMKBAIdvQwAAHAZqqqq1L9//1bP75ax06NHD0lfHpbH4+ng1QAAgNaIxWIKBALxv+Ot1S1j58JLVx6Ph9gBAKCLudS3oPAGZQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWI3YAAIDViB0AAGA1YgcAAFiN2AEAAFYjdgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWI3YAAIDViB0AAGA1YgcAAFiN2AEAAFYjdgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWI3YAAIDViB0AAGA1YgcAAFiN2AEAAFYjdgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWa5fYKSgo0MCBA5WWlqacnBzt2bOnxfllZWUaOXKk3G63Bg8erKKioovO3bJli1wulyZNmtTGqwYAADZwPHZefPFFLVq0SCtWrNC+ffs0YsQI5eXl6cSJE83Or6ys1IQJE3TrrbeqoqJCCxYs0OzZs7Vz584mcz/++GM9+uijGjdunNPbAAAAXZTLGGOcfICcnBzdeOON2rhxoySpsbFRgUBA8+fPV35+fpP5S5cuVXFxsQ4cOBAfmzp1qmpqalRSUhIfa2ho0Pe+9z09+OCDeuONN1RTU6Pt27e3ak2xWExer1fRaFQej+d/3CEAAGgPl/v329Fnds6ePau9e/cqGAz+9wGTkhQMBhUKhZq9JhQKJcyXpLy8vCbzn3jiCfXt21ezZs36xnXU19crFosl3AAAQPfgaOycPHlSDQ0N8vl8CeM+n0/hcLjZa8LhcLPzY7GYzpw5I0l688039eyzz6qwsLBV61i1apW8Xm/8FggELmM3AACgK+pyn8aqra3Vfffdp8LCQvXu3btV1yxbtkzRaDR+q6qqcniVAACgs0hx8s579+6t5ORkRSKRhPFIJCK/39/sNX6/v9n5Ho9H6enpqqio0Mcff6y77747/vvGxkZJUkpKio4cOaKrr7464Xq32y23290WWwIAAF2Mo8/spKamatSoUSotLY2PNTY2qrS0VLm5uc1ek5ubmzBfknbt2hWfP3ToUO3fv18VFRXx2z333BP/9BYvUQEAgK9y9JkdSVq0aJEeeOABjR49WmPGjNG6detUV1enmTNnSvryJabPPvtMf/rTnyRJDz/8sDZu3KglS5bowQcf1O7du7V161YVFxdLktLS0nT99dcnPEZGRoYkNRkHAABwPHamTJmizz//XMuXL1c4HFZ2drZKSkrib0I+fvy4Pvnkk/j8rKwsFRcXa+HChVq/fr369++vzZs3Ky8vz+mlAgAACzn+PTudEd+zAwBA19Mpv2cHAACgoxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKu1S+wUFBRo4MCBSktLU05Ojvbs2dPi/LKyMo0cOVJut1uDBw9WUVFRwu8LCws1btw49ezZUz179lQwGPzG+wQAAN2T47Hz4osvatGiRVqxYoX27dunESNGKC8vTydOnGh2fmVlpSZMmKBbb71VFRUVWrBggWbPnq2dO3fG55SVlWnatGl67bXXFAqFFAgEdPvtt+uzzz5zejsAAKCLcRljjJMPkJOToxtvvFEbN26UJDU2NioQCGj+/PnKz89vMn/p0qUqLi7WgQMH4mNTp05VTU2NSkpKmn2MhoYG9ezZUxs3btT999/f5Pf19fWqr6+P/xyLxRQIBBSNRuXxeP7XLQIAgHYQi8Xk9Xov+e+3o8/snD17Vnv37lUwGPzvAyYlKRgMKhQKNXtNKBRKmC9JeXl5F50vSadPn9a5c+fUq1evZn+/atUqeb3e+C0QCFzGbgAAQFfkaOycPHlSDQ0N8vl8CeM+n0/hcLjZa8LhcLPzY7GYzpw50+w1S5cuVWZmZpNIumDZsmWKRqPxW1VV1WXsBgAAdEUpHb2A/9Xq1au1ZcsWlZWVKS0trdk5brdbbre7nVcGAAA6A0djp3fv3kpOTlYkEkkYj0Qi8vv9zV7j9/ubne/xeJSenp4wvnbtWq1evVqvvvqqbrjhhrZdPAAAsIKjL2OlpqZq1KhRKi0tjY81NjaqtLRUubm5zV6Tm5ubMF+Sdu3a1WT+mjVr9OSTT6qkpESjR49u+8UDAAArOP7R80WLFqmwsFDPP/+8Dh06pLlz56qurk4zZ86U9OX7ab76CaqHH35YR48e1ZIlS3T48GE988wz2rp1qxYuXBif89RTT+nxxx/Xc889p4EDByocDiscDuvUqVNObwcAAHQxjr9nZ8qUKfr888+1fPlyhcNhZWdnq6SkJP4m5OPHj+uTTz6Jz8/KylJxcbEWLlyo9evXq3///tq8ebPy8vLic/7whz/o7Nmz+vGPf5zwWCtWrNDKlSud3hIAAOhCHP+enc7ocj+nDwAAOk6n/J4dAACAjkbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALBau8ROQUGBBg4cqLS0NOXk5GjPnj0tzi8rK9PIkSPldrs1ePBgFRUVNZnz0ksvaejQoUpLS9Pw4cP1yiuvOLR6AADQlTkeOy+++KIWLVqkFStWaN++fRoxYoTy8vJ04sSJZudXVlZqwoQJuvXWW1VRUaEFCxZo9uzZ2rlzZ3zOP//5T02bNk2zZs3Se++9p0mTJmnSpEk6cOCA09sBAABdjMsYY5x8gJycHN14443auHGjJKmxsVGBQEDz589Xfn5+k/lLly5VcXFxQrhMnTpVNTU1KikpkSRNmTJFdXV1evnll+NzbrrpJmVnZ2vTpk1N7rO+vl719fXxn2OxmAKBgKLRqDweT5vtFQAAOCcWi8nr9V7y329Hn9k5e/as9u7dq2Aw+N8HTEpSMBhUKBRq9ppQKJQwX5Ly8vIS5rdmzletWrVKXq83fgsEApe7JQAA0MU4GjsnT55UQ0ODfD5fwrjP51M4HG72mnA43Oz8WCymM2fOtDjnYve5bNkyRaPR+K2qqupytwQAALqYlI5eQHtwu91yu90dvQwAANABHH1mp3fv3kpOTlYkEkkYj0Qi8vv9zV7j9/ubne/xeJSent7inIvdJwAA6L4cjZ3U1FSNGjVKpaWl8bHGxkaVlpYqNze32Wtyc3MT5kvSrl27Eua3Zg4AAIDUDh89X7RokQoLC/X888/r0KFDmjt3rurq6jRz5kxJX76f5v7774/Pf/jhh3X06FEtWbJEhw8f1jPPPKOtW7dq4cKF8Tm/+MUvVFJSot/+9rc6fPiwVq5cqXfffVfz5s1zejsAAKCLcfw9O1OmTNHnn3+u5cuXKxwOKzs7WyUlJfE3GB8/flyffPJJfH5WVpaKi4u1cOFCrV+/Xv3799fmzZuVl5cXnzN27Fj9+c9/1mOPPaZf/epXGjJkiLZv367rr7/e6e0AAIAuxvHv2emMLvdz+gAAoON0yu/ZAQAA6GjEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrORY71dXVmj59ujwejzIyMjRr1iydOnWqxWuMMVq+fLn69eun9PR0BYNBffDBBwn3OX/+fF1zzTVKT0/XgAED9MgjjygajTq1DQAA0MU5FjvTp0/X+++/r127dunll1/W66+/rjlz5rR4zZo1a7RhwwZt2rRJ5eXluvLKK5WXl6cvvvhCknTs2DEdO3ZMa9eu1YEDB1RUVKSSkhLNmjXLqW0AAIAuzmWMMW19p4cOHdKwYcP0zjvvaPTo0ZKkkpIS3Xnnnfr000+VmZnZ5BpjjDIzM7V48WI9+uijkqRoNCqfz6eioiJNnTq12cd66aWX9NOf/lR1dXVKSUlp1fpisZi8Xq+i0ag8Hs9l7hIAALSny/377cgzO6FQSBkZGfHQkaRgMKikpCSVl5c3e01lZaXC4bCCwWB8zOv1KicnR6FQ6KKPdWHDLYVOfX29YrFYwg0AAHQPjsROOBxW3759E8ZSUlLUq1cvhcPhi14jST6fL2Hc5/Nd9JqTJ0/qySef/MaXx1atWiWv1xu/BQKB1m4FAAB0cZcUO/n5+XK5XC3eDh8+7NRaE8RiMU2YMEHDhg3TypUrW5y7bNkyRaPR+K2qqqpd1ggAADpe697k8v8tXrxYM2bMaHHOoEGD5Pf7deLEiYTx8+fPq7q6Wn6/v9nrLoxHIhH169cvPh6JRJSdnZ0wt7a2VuPHj1ePHj20bds2XXHFFS2uye12y+12tzgHAADY6ZJip0+fPurTp883zsvNzVVNTY327t2rUaNGSZJ2796txsZG5eTkNHtNVlaW/H6/SktL43ETi8VUXl6uuXPnxufFYjHl5eXJ7XZrx44dSktLu5QtAACAbsaR9+xce+21Gj9+vB566CHt2bNHb731lubNm6epU6cmfBJr6NCh2rZtmyTJ5XJpwYIF+vWvf60dO3Zo//79uv/++5WZmalJkyZJ+jJ0br/9dtXV1enZZ59VLBZTOBxWOBxWQ0ODE1sBAABd3CU9s3MpXnjhBc2bN0+33XabkpKSNHnyZG3YsCFhzpEjRxK+EHDJkiWqq6vTnDlzVFNTo5tvvlklJSXxZ2/27dsX/zTX4MGDE+6rsrJSAwcOdGo7AACgi3Lke3Y6O75nBwCArqdTfc8OAABAZ0HsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALCaY7FTXV2t6dOny+PxKCMjQ7NmzdKpU6davMYYo+XLl6tfv35KT09XMBjUBx98cNG5d9xxh1wul7Zv3+7EFgAAgAUci53p06fr/fff165du/Tyyy/r9ddf15w5c1q8Zs2aNdqwYYM2bdqk8vJyXXnllcrLy9MXX3zRZO66devkcrmcWj4AALCEyxhj2vpODx06pGHDhumdd97R6NGjJUklJSW688479emnnyozM7PJNcYYZWZmavHixXr00UclSdFoVD6fT0VFRZo6dWp8bkVFhe666y69++676tevn7Zt26ZJkya1en2xWExer1fRaFQej+d/3C0AAGgPl/v325FndkKhkDIyMuKhI0nBYFBJSUkqLy9v9prKykqFw2EFg8H4mNfrVU5OjkKhUHzs9OnTuvfee1VQUCC/39+q9dTX1ysWiyXcAABA9+BI7ITDYfXt2zdhLCUlRb169VI4HL7oNZLk8/kSxn0+X8I1Cxcu1NixYzVx4sRWr2fVqlXyer3xWyAQaPW1AACga7uk2MnPz5fL5WrxdvjwYafWqh07dmj37t1at27dJV23bNkyRaPR+K2qqsqhFQIAgM4m5VImL168WDNmzGhxzqBBg+T3+3XixImE8fPnz6u6uvqiLz1dGI9EIurXr198PBKJKDs7W5K0e/duffTRR8rIyEi4dvLkyRo3bpzKysqavW+32y23293iugEAgJ0uKXb69OmjPn36fOO83Nxc1dTUaO/evRo1apSkL0OlsbFROTk5zV6TlZUlv9+v0tLSeNzEYjGVl5dr7ty5kr58Zmn27NkJ1w0fPly///3vdffdd1/KVgAAQDdxSbHTWtdee63Gjx+vhx56SJs2bdK5c+c0b948TZ06NeGTWEOHDtWqVav0wx/+UC6XSwsWLNCvf/1rDRkyRFlZWXr88ceVmZkZ/6SV3+9v9pmhAQMGKCsry4mtAACALs6R2JGkF154QfPmzdNtt92mpKQkTZ48WRs2bEiYc+TIEUWj0fjPS5YsUV1dnebMmaOamhrdfPPNKikpUVpamlPLBAAAlnPke3Y6O75nBwCArqdTfc8OAABAZ0HsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrpXT0AjqCMUaSFIvFOnglAACgtS783b7wd7y1umXs1NbWSpICgUAHrwQAAFyq2tpaeb3eVs93mUvNIws0Njbq2LFj6tGjh1wuV0cvp8PFYjEFAgFVVVXJ4/F09HKsxTm3D865/XDW7YNz/i9jjGpra5WZmamkpNa/E6dbPrOTlJSk/v37d/QyOh2Px9Pt/4fUHjjn9sE5tx/Oun1wzl+6lGd0LuANygAAwGrEDgAAsFryypUrV3b0ItDxkpOT9f3vf18pKd3ylc12wzm3D865/XDW7YNz/t90yzcoAwCA7oOXsQAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YidbqC6ulrTp0+Xx+NRRkaGZs2apVOnTrV4jTFGy5cvV79+/ZSenq5gMKgPPvjgonPvuOMOuVwubd++3YktdAlOnHN1dbXmz5+va665Runp6RowYIAeeeQRRaNRp7fTqRQUFGjgwIFKS0tTTk6O9uzZ0+L8srIyjRw5Um63W4MHD1ZRUVGTOS+99JKGDh2qtLQ0DR8+XK+88opDq+862vqcCwsLNW7cOPXs2VM9e/ZUMBj8xvvsDpz493zBli1b5HK5NGnSpDZedRdnYL3x48ebESNGmLffftu88cYbZvDgwWbatGktXrN69Wrj9XrN9u3bzb/+9S9zzz33mKysLHPmzJkmc3/3u9+ZO+64w0gy27Ztc2obnZ4T57x//37zox/9yOzYscN8+OGHprS01AwZMsRMnjy5PbbUKWzZssWkpqaa5557zrz//vvmoYceMhkZGSYSiTQ7/+jRo+Zb3/qWWbRokTl48KB5+umnTXJysikpKYnPeeutt0xycrJZs2aNOXjwoHnsscfMFVdcYfbv399e2+p0nDjne++91xQUFJj33nvPHDp0yMyYMcN4vV7z6aeftte2Oh0nzvmCyspKc9VVV5lx48aZiRMnOr2VLoXYsdzBgweNJPPOO+/Ex/7xj38Yl8tlPvvss2avaWxsNH6/3/zmN7+Jj9XU1Bi3223+8pe/JMx97733zFVXXWWOHz/erWPH6XP+qq1bt5rU1FRz7ty5tttAJzZmzBjz85//PP5zQ0ODyczMNKtWrWp2/pIlS8x1112XMDZlyhSTl5cX//knP/mJmTBhQsKcnJwc87Of/awNV961OHHOX3f+/HnTo0cP8/zzz7fNorsgp875/PnzZuzYsWbz5s3mgQceIHa+hpexLBcKhZSRkaHRo0fHx4LBoJKSklReXt7sNZWVlQqHwwoGg/Exr9ernJwchUKh+Njp06d17733qqCgQH6/37lNdAFOnvPXRaNReTyebvFNqmfPntXevXsTzigpKUnBYPCiZxQKhRLmS1JeXl7C/NbM6U6cOuevO336tM6dO6devXq1zcK7GCfP+YknnlDfvn01a9astl+4BYgdy4XDYfXt2zdhLCUlRb169VI4HL7oNZLk8/kSxn0+X8I1Cxcu1NixYzVx4sQ2XnXX4+Q5f9XJkyf15JNPas6cOW2w6s7v5MmTamhouKQzCofDzc6PxWI6c+ZMi3Mudp+2c+qcv27p0qXKzMxs8se7u3DqnN988009++yzKiwsdGbhFiB2uqj8/Hy5XK4Wb4cPH3bs8Xfs2KHdu3dr3bp1jj1GZ9DR5/xVsVhMEyZM0LBhw8T/pR26mtWrV2vLli3atm2b0tLSOno51qitrdV9992nwsJC9e7du6OX02nZ/zy4pRYvXqwZM2a0OGfQoEHy+/06ceJEwvj58+dVXV190ZeeLoxHIhH169cvPh6JRJSdnS1J2r17tz766CNlZGQkXDt58mSNGzdOZWVll7ijzqmjz/mC2tpajR8/Xj169NC2bdt0xRVXXMZuup7evXsrOTlZkUgkYTwSibR4rs3N93g8Sk9Pb3FOd3051qlzvmDt2rVavXq1Xn31Vd1www1tu/guxIlzrqio0Mcff6y77747/vvGxkZJXz67fOTIEV199dVtvJMuqKPfNARnXXjj7Lvvvhsf27lzZ6veOLt27dr4WDQaTXjj7PHjx83+/fsTbpLM+vXrzdGjR53dVCfk1DlfGLvpppvMLbfcYurq6pzbRCc1ZswYM2/evPjPDQ0N5qqrrmrxDZ3XX399wti0adOavEH5rrvuSpiTm5vb7d+g3NbnbIwxTz31lPF4PCYUCrX9orugtj7nM2fONPlv8cSJE80PfvADs3//flNfX+/cZroQYqcbGD9+vPnud79rysvLzZtvvmmGDBnS5CPR11xzjfnb3/4W/3n16tUmIyPD/P3vfzf//ve/zcSJEy/60fML1I0/jWWMM+ccjUZNTk6OGT58uPnwww/N8ePH47fz58+36/46ypYtW4zb7TZFRUXm4MGDZs6cOSYjI8OEw2FjjDH5+fnmvvvui8+/8FHdX/7yl+bQoUOmoKCg2Y+ep6SkmLVr15pDhw6ZFStW8NFzB8559erVJjU11fz1r39N+LdbW1vb7vvrLJw456/j01hNETvdwH/+8x8zbdo08+1vf9t4PB4zc+bMJv+xkWT++Mc/xn9ubGw0jz/+uPH5fMbtdpvbbrvNHDlypMXH6e6x48Q5v/baa0ZSs7fKysp22lnHe/rpp82AAQNMamqqGTNmjHn77bfjv3vggQfMLbfckjD/tddeM9nZ2SY1NdUMGjQo4cwv2Lp1q/nOd75jUlNTzXXXXWeKi4sd3kXn19bn/H//93/N/ttdsWKF85vpxJz49/xVxE5TLmOMae+XzgAAANoLn8YCAABWI3YAAIDViB0AAGA1YgcAAFiN2AEAAFYjdgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgtf8HyT6sQPPu5OoAAAAASUVORK5CYII=", "text/plain": [ "PyPlot.Figure(PyObject )" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "1-element Array{Any,1}:\n", " PyObject " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(log10(En_array - minimum(En_array)))" ] }, { "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": [ "

\n", " \n", "__Exercise: test the recovery of a less sparse signal.\n", "What do you observe?__

" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "100-element Array{Float64,1}:\n", " 0.11573 \n", " 0.742032 \n", " 1.1382 \n", " 0.0195104 \n", " -0.0564791 \n", " -0.303705 \n", " -0.400856 \n", " 0.131482 \n", " -0.31011 \n", " -0.00970144\n", " -0.0424864 \n", " -0.0766934 \n", " 0.619124 \n", " ⋮ \n", " -0.210507 \n", " 0.502878 \n", " -0.785536 \n", " -0.682072 \n", " -0.991464 \n", " 0.409874 \n", " -0.938657 \n", " 0.257866 \n", " -0.129917 \n", " 0.105005 \n", " 0.036563 \n", " 0.0194049 " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S = 31\n", "srand(0)\n", "sel = randperm(N)\n", "sel = sel[1 : S] # indices of the nonzero elements of xsharp\n", "xsharp = zeros(N)\n", "xsharp[sel] = 1\n", "\n", "y = A*xsharp\n", "\n", "# put your code here" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "ename": "LoadError", "evalue": "UndefVarError: x_restored not defined", "output_type": "error", "traceback": [ "UndefVarError: x_restored not defined", "" ] } ], "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](L\"$x^\\sharp$\")\n", "subfig2[:set_title](L\"$x_\\mathrm{restored}$\")" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGgCAYAAABMn6ZGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAH7RJREFUeJzt3X9s1dXh//HXbWtvO8e9hQD3Url8KMJEFOkAqSUy57yxKCpsLANkKggyycDxw0FZFIguARnbAKkjKbq6xA1xGYxYV4LFxh+7FgW7gfyISo1VuBdZ03tLwQLt+f7hlzuvLbWwvvvj9PlI7h89Pe97zzlh9pn7ay5jjBEAAIClkjp6AQAAAE4idgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWI3YAAIDViB0AAGC1lI5eQEdobGzUsWPH1KNHD7lcro5eDgAAaAVjjGpra5WZmamkpNY/X9MtY+fYsWMKBAIdvQwAAHAZqqqq1L9//1bP75ax06NHD0lfHpbH4+ng1QAAgNaIxWIKBALxv+Ot1S1j58JLVx6Ph9gBAKCLudS3oPAGZQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWI3YAAIDViB0AAGA1YgcAAFiN2AEAAFYjdgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWI3YAAIDViB0AAGA1YgcAAFiN2AEAAFYjdgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWI3YAAIDViB0AAGA1YgcAAFiN2AEAAFYjdgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgNWIHAABYjdgBAABWa5fYKSgo0MCBA5WWlqacnBzt2bOnxfllZWUaOXKk3G63Bg8erKKioovO3bJli1wulyZNmtTGqwYAADZwPHZefPFFLVq0SCtWrNC+ffs0YsQI5eXl6cSJE83Or6ys1IQJE3TrrbeqoqJCCxYs0OzZs7Vz584mcz/++GM9+uijGjdunNPbAAAAXZTLGGOcfICcnBzdeOON2rhxoySpsbFRgUBA8+fPV35+fpP5S5cuVXFxsQ4cOBAfmzp1qmpqalRSUhIfa2ho0Pe+9z09+OCDeuONN1RTU6Pt27e3ak2xWExer1fRaFQej+d/3CEAAGgPl/v329Fnds6ePau9e/cqGAz+9wGTkhQMBhUKhZq9JhQKJcyXpLy8vCbzn3jiCfXt21ezZs36xnXU19crFosl3AAAQPfgaOycPHlSDQ0N8vl8CeM+n0/hcLjZa8LhcLPzY7GYzpw5I0l688039eyzz6qwsLBV61i1apW8Xm/8FggELmM3AACgK+pyn8aqra3Vfffdp8LCQvXu3btV1yxbtkzRaDR+q6qqcniVAACgs0hx8s579+6t5ORkRSKRhPFIJCK/39/sNX6/v9n5Ho9H6enpqqio0Mcff6y77747/vvGxkZJUkpKio4cOaKrr7464Xq32y23290WWwIAAF2Mo8/spKamatSoUSotLY2PNTY2qrS0VLm5uc1ek5ubmzBfknbt2hWfP3ToUO3fv18VFRXx2z333BP/9BYvUQEAgK9y9JkdSVq0aJEeeOABjR49WmPGjNG6detUV1enmTNnSvryJabPPvtMf/rTnyRJDz/8sDZu3KglS5bowQcf1O7du7V161YVFxdLktLS0nT99dcnPEZGRoYkNRkHAABwPHamTJmizz//XMuXL1c4HFZ2drZKSkrib0I+fvy4Pvnkk/j8rKwsFRcXa+HChVq/fr369++vzZs3Ky8vz+mlAgAACzn+PTudEd+zAwBA19Mpv2cHAACgoxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKu1S+wUFBRo4MCBSktLU05Ojvbs2dPi/LKyMo0cOVJut1uDBw9WUVFRwu8LCws1btw49ezZUz179lQwGPzG+wQAAN2T47Hz4osvatGiRVqxYoX27dunESNGKC8vTydOnGh2fmVlpSZMmKBbb71VFRUVWrBggWbPnq2dO3fG55SVlWnatGl67bXXFAqFFAgEdPvtt+uzzz5zejsAAKCLcRljjJMPkJOToxtvvFEbN26UJDU2NioQCGj+/PnKz89vMn/p0qUqLi7WgQMH4mNTp05VTU2NSkpKmn2MhoYG9ezZUxs3btT999/f5Pf19fWqr6+P/xyLxRQIBBSNRuXxeP7XLQIAgHYQi8Xk9Xov+e+3o8/snD17Vnv37lUwGPzvAyYlKRgMKhQKNXtNKBRKmC9JeXl5F50vSadPn9a5c+fUq1evZn+/atUqeb3e+C0QCFzGbgAAQFfkaOycPHlSDQ0N8vl8CeM+n0/hcLjZa8LhcLPzY7GYzpw50+w1S5cuVWZmZpNIumDZsmWKRqPxW1VV1WXsBgAAdEUpHb2A/9Xq1au1ZcsWlZWVKS0trdk5brdbbre7nVcGAAA6A0djp3fv3kpOTlYkEkkYj0Qi8vv9zV7j9/ubne/xeJSenp4wvnbtWq1evVqvvvqqbrjhhrZdPAAAsIKjL2OlpqZq1KhRKi0tjY81NjaqtLRUubm5zV6Tm5ubMF+Sdu3a1WT+mjVr9OSTT6qkpESjR49u+8UDAAArOP7R80WLFqmwsFDPP/+8Dh06pLlz56qurk4zZ86U9OX7ab76CaqHH35YR48e1ZIlS3T48GE988wz2rp1qxYuXBif89RTT+nxxx/Xc889p4EDByocDiscDuvUqVNObwcAAHQxjr9nZ8qUKfr888+1fPlyhcNhZWdnq6SkJP4m5OPHj+uTTz6Jz8/KylJxcbEWLlyo9evXq3///tq8ebPy8vLic/7whz/o7Nmz+vGPf5zwWCtWrNDKlSud3hIAAOhCHP+enc7ocj+nDwAAOk6n/J4dAACAjkbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALBau8ROQUGBBg4cqLS0NOXk5GjPnj0tzi8rK9PIkSPldrs1ePBgFRUVNZnz0ksvaejQoUpLS9Pw4cP1yiuvOLR6AADQlTkeOy+++KIWLVqkFStWaN++fRoxYoTy8vJ04sSJZudXVlZqwoQJuvXWW1VRUaEFCxZo9uzZ2rlzZ3zOP//5T02bNk2zZs3Se++9p0mTJmnSpEk6cOCA09sBAABdjMsYY5x8gJycHN14443auHGjJKmxsVGBQEDz589Xfn5+k/lLly5VcXFxQrhMnTpVNTU1KikpkSRNmTJFdXV1evnll+NzbrrpJmVnZ2vTpk1N7rO+vl719fXxn2OxmAKBgKLRqDweT5vtFQAAOCcWi8nr9V7y329Hn9k5e/as9u7dq2Aw+N8HTEpSMBhUKBRq9ppQKJQwX5Ly8vIS5rdmzletWrVKXq83fgsEApe7JQAA0MU4GjsnT55UQ0ODfD5fwrjP51M4HG72mnA43Oz8WCymM2fOtDjnYve5bNkyRaPR+K2qqupytwQAALqYlI5eQHtwu91yu90dvQwAANABHH1mp3fv3kpOTlYkEkkYj0Qi8vv9zV7j9/ubne/xeJSent7inIvdJwAA6L4cjZ3U1FSNGjVKpaWl8bHGxkaVlpYqNze32Wtyc3MT5kvSrl27Eua3Zg4AAIDUDh89X7RokQoLC/X888/r0KFDmjt3rurq6jRz5kxJX76f5v7774/Pf/jhh3X06FEtWbJEhw8f1jPPPKOtW7dq4cKF8Tm/+MUvVFJSot/+9rc6fPiwVq5cqXfffVfz5s1zejsAAKCLcfw9O1OmTNHnn3+u5cuXKxwOKzs7WyUlJfE3GB8/flyffPJJfH5WVpaKi4u1cOFCrV+/Xv3799fmzZuVl5cXnzN27Fj9+c9/1mOPPaZf/epXGjJkiLZv367rr7/e6e0AAIAuxvHv2emMLvdz+gAAoON0yu/ZAQAA6GjEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrORY71dXVmj59ujwejzIyMjRr1iydOnWqxWuMMVq+fLn69eun9PR0BYNBffDBBwn3OX/+fF1zzTVKT0/XgAED9MgjjygajTq1DQAA0MU5FjvTp0/X+++/r127dunll1/W66+/rjlz5rR4zZo1a7RhwwZt2rRJ5eXluvLKK5WXl6cvvvhCknTs2DEdO3ZMa9eu1YEDB1RUVKSSkhLNmjXLqW0AAIAuzmWMMW19p4cOHdKwYcP0zjvvaPTo0ZKkkpIS3Xnnnfr000+VmZnZ5BpjjDIzM7V48WI9+uijkqRoNCqfz6eioiJNnTq12cd66aWX9NOf/lR1dXVKSUlp1fpisZi8Xq+i0ag8Hs9l7hIAALSny/377cgzO6FQSBkZGfHQkaRgMKikpCSVl5c3e01lZaXC4bCCwWB8zOv1KicnR6FQ6KKPdWHDLYVOfX29YrFYwg0AAHQPjsROOBxW3759E8ZSUlLUq1cvhcPhi14jST6fL2Hc5/Nd9JqTJ0/qySef/MaXx1atWiWv1xu/BQKB1m4FAAB0cZcUO/n5+XK5XC3eDh8+7NRaE8RiMU2YMEHDhg3TypUrW5y7bNkyRaPR+K2qqqpd1ggAADpe697k8v8tXrxYM2bMaHHOoEGD5Pf7deLEiYTx8+fPq7q6Wn6/v9nrLoxHIhH169cvPh6JRJSdnZ0wt7a2VuPHj1ePHj20bds2XXHFFS2uye12y+12tzgHAADY6ZJip0+fPurTp883zsvNzVVNTY327t2rUaNGSZJ2796txsZG5eTkNHtNVlaW/H6/SktL43ETi8VUXl6uuXPnxufFYjHl5eXJ7XZrx44dSktLu5QtAACAbsaR9+xce+21Gj9+vB566CHt2bNHb731lubNm6epU6cmfBJr6NCh2rZtmyTJ5XJpwYIF+vWvf60dO3Zo//79uv/++5WZmalJkyZJ+jJ0br/9dtXV1enZZ59VLBZTOBxWOBxWQ0ODE1sBAABd3CU9s3MpXnjhBc2bN0+33XabkpKSNHnyZG3YsCFhzpEjRxK+EHDJkiWqq6vTnDlzVFNTo5tvvlklJSXxZ2/27dsX/zTX4MGDE+6rsrJSAwcOdGo7AACgi3Lke3Y6O75nBwCArqdTfc8OAABAZ0HsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALCaY7FTXV2t6dOny+PxKCMjQ7NmzdKpU6davMYYo+XLl6tfv35KT09XMBjUBx98cNG5d9xxh1wul7Zv3+7EFgAAgAUci53p06fr/fff165du/Tyyy/r9ddf15w5c1q8Zs2aNdqwYYM2bdqk8vJyXXnllcrLy9MXX3zRZO66devkcrmcWj4AALCEyxhj2vpODx06pGHDhumdd97R6NGjJUklJSW688479emnnyozM7PJNcYYZWZmavHixXr00UclSdFoVD6fT0VFRZo6dWp8bkVFhe666y69++676tevn7Zt26ZJkya1en2xWExer1fRaFQej+d/3C0AAGgPl/v325FndkKhkDIyMuKhI0nBYFBJSUkqLy9v9prKykqFw2EFg8H4mNfrVU5OjkKhUHzs9OnTuvfee1VQUCC/39+q9dTX1ysWiyXcAABA9+BI7ITDYfXt2zdhLCUlRb169VI4HL7oNZLk8/kSxn0+X8I1Cxcu1NixYzVx4sRWr2fVqlXyer3xWyAQaPW1AACga7uk2MnPz5fL5WrxdvjwYafWqh07dmj37t1at27dJV23bNkyRaPR+K2qqsqhFQIAgM4m5VImL168WDNmzGhxzqBBg+T3+3XixImE8fPnz6u6uvqiLz1dGI9EIurXr198PBKJKDs7W5K0e/duffTRR8rIyEi4dvLkyRo3bpzKysqavW+32y23293iugEAgJ0uKXb69OmjPn36fOO83Nxc1dTUaO/evRo1apSkL0OlsbFROTk5zV6TlZUlv9+v0tLSeNzEYjGVl5dr7ty5kr58Zmn27NkJ1w0fPly///3vdffdd1/KVgAAQDdxSbHTWtdee63Gjx+vhx56SJs2bdK5c+c0b948TZ06NeGTWEOHDtWqVav0wx/+UC6XSwsWLNCvf/1rDRkyRFlZWXr88ceVmZkZ/6SV3+9v9pmhAQMGKCsry4mtAACALs6R2JGkF154QfPmzdNtt92mpKQkTZ48WRs2bEiYc+TIEUWj0fjPS5YsUV1dnebMmaOamhrdfPPNKikpUVpamlPLBAAAlnPke3Y6O75nBwCArqdTfc8OAABAZ0HsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrETsAAMBqxA4AALAasQMAAKxG7AAAAKsROwAAwGrEDgAAsBqxAwAArEbsAAAAqxE7AADAasQOAACwGrEDAACsRuwAAACrpXT0AjqCMUaSFIvFOnglAACgtS783b7wd7y1umXs1NbWSpICgUAHrwQAAFyq2tpaeb3eVs93mUvNIws0Njbq2LFj6tGjh1wuV0cvp8PFYjEFAgFVVVXJ4/F09HKsxTm3D865/XDW7YNz/i9jjGpra5WZmamkpNa/E6dbPrOTlJSk/v37d/QyOh2Px9Pt/4fUHjjn9sE5tx/Oun1wzl+6lGd0LuANygAAwGrEDgAAsFryypUrV3b0ItDxkpOT9f3vf18pKd3ylc12wzm3D865/XDW7YNz/t90yzcoAwCA7oOXsQAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YidbqC6ulrTp0+Xx+NRRkaGZs2apVOnTrV4jTFGy5cvV79+/ZSenq5gMKgPPvjgonPvuOMOuVwubd++3YktdAlOnHN1dbXmz5+va665Runp6RowYIAeeeQRRaNRp7fTqRQUFGjgwIFKS0tTTk6O9uzZ0+L8srIyjRw5Um63W4MHD1ZRUVGTOS+99JKGDh2qtLQ0DR8+XK+88opDq+862vqcCwsLNW7cOPXs2VM9e/ZUMBj8xvvsDpz493zBli1b5HK5NGnSpDZedRdnYL3x48ebESNGmLffftu88cYbZvDgwWbatGktXrN69Wrj9XrN9u3bzb/+9S9zzz33mKysLHPmzJkmc3/3u9+ZO+64w0gy27Ztc2obnZ4T57x//37zox/9yOzYscN8+OGHprS01AwZMsRMnjy5PbbUKWzZssWkpqaa5557zrz//vvmoYceMhkZGSYSiTQ7/+jRo+Zb3/qWWbRokTl48KB5+umnTXJysikpKYnPeeutt0xycrJZs2aNOXjwoHnsscfMFVdcYfbv399e2+p0nDjne++91xQUFJj33nvPHDp0yMyYMcN4vV7z6aeftte2Oh0nzvmCyspKc9VVV5lx48aZiRMnOr2VLoXYsdzBgweNJPPOO+/Ex/7xj38Yl8tlPvvss2avaWxsNH6/3/zmN7+Jj9XU1Bi3223+8pe/JMx97733zFVXXWWOHz/erWPH6XP+qq1bt5rU1FRz7ty5tttAJzZmzBjz85//PP5zQ0ODyczMNKtWrWp2/pIlS8x1112XMDZlyhSTl5cX//knP/mJmTBhQsKcnJwc87Of/awNV961OHHOX3f+/HnTo0cP8/zzz7fNorsgp875/PnzZuzYsWbz5s3mgQceIHa+hpexLBcKhZSRkaHRo0fHx4LBoJKSklReXt7sNZWVlQqHwwoGg/Exr9ernJwchUKh+Njp06d17733qqCgQH6/37lNdAFOnvPXRaNReTyebvFNqmfPntXevXsTzigpKUnBYPCiZxQKhRLmS1JeXl7C/NbM6U6cOuevO336tM6dO6devXq1zcK7GCfP+YknnlDfvn01a9astl+4BYgdy4XDYfXt2zdhLCUlRb169VI4HL7oNZLk8/kSxn0+X8I1Cxcu1NixYzVx4sQ2XnXX4+Q5f9XJkyf15JNPas6cOW2w6s7v5MmTamhouKQzCofDzc6PxWI6c+ZMi3Mudp+2c+qcv27p0qXKzMxs8se7u3DqnN988009++yzKiwsdGbhFiB2uqj8/Hy5XK4Wb4cPH3bs8Xfs2KHdu3dr3bp1jj1GZ9DR5/xVsVhMEyZM0LBhw8T/pR26mtWrV2vLli3atm2b0tLSOno51qitrdV9992nwsJC9e7du6OX02nZ/zy4pRYvXqwZM2a0OGfQoEHy+/06ceJEwvj58+dVXV190ZeeLoxHIhH169cvPh6JRJSdnS1J2r17tz766CNlZGQkXDt58mSNGzdOZWVll7ijzqmjz/mC2tpajR8/Xj169NC2bdt0xRVXXMZuup7evXsrOTlZkUgkYTwSibR4rs3N93g8Sk9Pb3FOd3051qlzvmDt2rVavXq1Xn31Vd1www1tu/guxIlzrqio0Mcff6y77747/vvGxkZJXz67fOTIEV199dVtvJMuqKPfNARnXXjj7Lvvvhsf27lzZ6veOLt27dr4WDQaTXjj7PHjx83+/fsTbpLM+vXrzdGjR53dVCfk1DlfGLvpppvMLbfcYurq6pzbRCc1ZswYM2/evPjPDQ0N5qqrrmrxDZ3XX399wti0adOavEH5rrvuSpiTm5vb7d+g3NbnbIwxTz31lPF4PCYUCrX9orugtj7nM2fONPlv8cSJE80PfvADs3//flNfX+/cZroQYqcbGD9+vPnud79rysvLzZtvvmmGDBnS5CPR11xzjfnb3/4W/3n16tUmIyPD/P3vfzf//ve/zcSJEy/60fML1I0/jWWMM+ccjUZNTk6OGT58uPnwww/N8ePH47fz58+36/46ypYtW4zb7TZFRUXm4MGDZs6cOSYjI8OEw2FjjDH5+fnmvvvui8+/8FHdX/7yl+bQoUOmoKCg2Y+ep6SkmLVr15pDhw6ZFStW8NFzB8559erVJjU11fz1r39N+LdbW1vb7vvrLJw456/j01hNETvdwH/+8x8zbdo08+1vf9t4PB4zc+bMJv+xkWT++Mc/xn9ubGw0jz/+uPH5fMbtdpvbbrvNHDlypMXH6e6x48Q5v/baa0ZSs7fKysp22lnHe/rpp82AAQNMamqqGTNmjHn77bfjv3vggQfMLbfckjD/tddeM9nZ2SY1NdUMGjQo4cwv2Lp1q/nOd75jUlNTzXXXXWeKi4sd3kXn19bn/H//93/N/ttdsWKF85vpxJz49/xVxE5TLmOMae+XzgAAANoLn8YCAABWI3YAAIDViB0AAGA1YgcAAFiN2AEAAFYjdgAAgNWIHQAAYDViBwAAWI3YAQAAViN2AACA1YgdAABgtf8HyT6sQPPu5OoAAAAASUVORK5CYII=", "text/plain": [ "PyPlot.Figure(PyObject )" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "1-element Array{Any,1}:\n", " PyObject " ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(log10(En_array - minimum(En_array)))" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Convergence is much slower in this setting." ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.5.0", "language": "julia", "name": "julia-0.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.5.0" } }, "nbformat": 4, "nbformat_minor": 0 }