{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Homotopy Continuation\n", "=====================\n", "\n*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_matlab/) for details about how to install the toolboxes.\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"], "metadata": {}, "cell_type": "markdown"}, {"source": ["This numerical tour details the homotopy algorithm to solve $\\ell^1$\n", "regularization."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 2, "cell_type": "code", "language": "python", "metadata": {}, "input": ["addpath('toolbox_signal')\n", "addpath('toolbox_general')\n", "addpath('solutions/optim_8_homotopy')"]}, {"source": ["$\\ell^1$ Regularization\n", "We consider the following linear inverse problem\n", "$$ y = \\Phi x_0 + w  \\in \\RR^P$$\n", "where $x_0 \\in \\RR^N$ is the (unknown) signal to recover, $w \\in\n", "\\RR^P$ is a noise vector, and $\\Phi \\in \\RR^{P \\times N}$ models the\n", "acquisition device.\n", "\n", "\n", "Size of the problem."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["N = 200;\n", "P = round(N/4);"]}, {"source": ["We consider here a simple compressed sensing scenario where the operator\n", "$\\Phi$ is a realization of a Gaussian random matrix."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Phi = randn(P,N);"]}, {"source": ["Sparsity of the input data."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["s = round(P/5);"]}, {"source": ["Sparse vector."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["I = randperm(N); I = I(1:s);\n", "x0 = zeros(N,1); x0(I) = sign(randn(s,1));"]}, {"source": ["Noise level."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sigma = 0.05 * norm(Phi*x0)/sqrt(P);"]}, {"source": ["Measurements $y=\\Phi x_0 + w$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["y = Phi*x0 + sigma*randn(P,1);"]}, {"source": ["To recover an approximation of the signal $x_0$, we use the Basis\n", "Pursuit denoising method, that use the $\\ell^1$ norm as a sparsity\n", "enforcing penalty\n", "$$ x_\\la \\in \\uargmin{x \\in \\RR^N} \\frac{1}{2} \\norm{y-\\Phi x}^2 + \\la \\norm{x}_1\n", "      \\quad (\\Pp_\\la), $$\n", "where the $\\ell^1$ norm is defined as\n", "$$ \\norm{x}_1 = \\sum_i \\abs{x_i}. $$\n", "\n", "\n", "The parameter $\\la$ should be set in accordance to the noise level\n", "$\\norm{w}$.\n", "\n", "\n", "In the case where there is no noise, $w=0$, we let $\\la \\rightarrow\n", "0^+$ and solve the Basis Pursuit constrained optimization\n", "$$ x_{0^+} \\in \\uargmin{\\Phi x= y} \\norm{x}_1. $$\n", "\n", "\n", "In order to avoid technical difficulties, we assume that $\\Phi$ is such\n", "that $x_\\la$ is uniquely defined. This is the case with probability 1\n", "when $\\Phi$ is a realization of a Gaussian matrix field.\n", "\n", "\n", "In the following, for some index set $I \\subset \\{1,\\ldots,N\\}$, we\n", "denote\n", "$$ \\Phi_I = (\\phi_i)_{i \\in I} \\in \\RR^{P \\times \\abs{I}} $$\n", "the sub-matrix obtained by extracting the columns\n", "$\\phi_i \\in \\RR^P$ of $\\Phi$ indexed by $I$.\n", "The support of a vector is\n", "$$ \\text{support}(x) =  \\enscond{i \\in \\{1,\\ldots,N\\}}{ x_i \\neq 0 }. $$\n", "\n", "\n", "Standard results from convex analysis shows that $x_\\la$ is a solution\n", "of $\\Pp_\\la$ if and only if\n", "$$  \\choice{\n", "      \\text{(C1)}  \\qquad \\Phi_I^*( y-\\Phi_I x_{\\la,I} ) = \\la \\text{sign}(x_{\\la,I}), \\\\\n", "      \\text{(C2)}  \\qquad \\norm{ \\Phi_{J}^*( y-\\Phi_I x_{\\la,I} ) }_\\infty \\leq \\la\n", "  }$\n", "$$\n", "where $I=\\text{support}(x_\\la)$ and $J = I^c$ is the complementary.\n", "\n", "Homotopy Algorithm\n", "------------------\n", "The homotopy algorithm computes the solution $x_\\la$ for all values of\n", "$\\la \\in \\RR^+$. More precisely, it can be shown that\n", "$$ \\la \\in \\RR^+ \\mapsto x_\\la \\in \\RR^N $$\n", "is a piecewise affine map, and the homotopy algorihtm computes the\n", "solution vectors\n", "$$ x_{\\la_1},\\ldots,x_{\\la_K} \\qwhereq 0^+ = \\la_K < \\ldots < \\la_{1} $$\n", "at the breakpoints $(x_{\\la_k})_k$ of this piecewise affine path.\n", "\n", "\n", "This algorithm was introduced in the over-determined setting in:\n", "\n", "\n", "Michael R. Osborne, Brett Presnell, and Berwin A. Turlach.\n", "_On the lasso and its dual_.\n", "Journal of Computational and Graphical Statistics, 9:319?337, 2000.\n", "\n", "\n", "It is extended to an arbitrary matrix $\\Phi$ in:\n", "\n", "\n", "Dmitry M. Malioutov, Mujdat Cetin, and Alan S. Willsky.\n", "_Homotopy continuation for sparse signal representation_.\n", "In IEEE Int. Conf. Acoustics, Speech and Signal Processing, Philadelphia, PA, volume 5, pages 733?736, March 2005.\n", "\n", "\n", "This algorithm is closely related to the LARS algorithm, that is an\n", "approximate method:\n", "\n", "\n", "Bradley Efron, Trevor Hastie, Iain M. Johnstone, and Robert Tibshirani.\n", "_Least angle regression_.\n", "The Annals of Statistics, 32(2):407?499, 2004.\n", "\n", "\n", "A theoretical analysis of the homotopy algorithm when the dictionary\n", "$\\Phi$ is incoherent is performed in:\n", "\n", "\n", "David L. Donoho, Yaakov Tsaig,\n", "_Fast Solution of l1-Norm Minimization Problems When the Solution May Be Sparse_.\n", "IEEE Transactions on Information Theory 54(11): 4789-4812, 2008.\n", "\n", "\n", "For large values of $\\la>\\la_1$, the solution is $x_\\la = 0$.\n", "The solution $0$ satisfies the constraints (C2) until\n", "$$ \\la = \\la_1 = \\norm{ \\Phi^* y }_\\infty $$\n", "The support is then a single Dirac\n", "$$ I(x_\\la) = \\{i_0\\} \\qwhereq \\abs{\\dotp{\\phi_{i_0}}{y}} = \\norm{ \\Phi^* y }_\\infty. $$\n", "\n", "\n", "We thus initialize the algorithm with this value of $\\la=\\la_1$ and $I$, together with the solution $x_\\la=0$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["C = Phi'*y;\n", "[lambda,I] = max(abs(C));\n", "x = zeros(N,1);"]}, {"source": ["One Iteration of the Algorithm\n", "------------------------------\n", "The algorithm proceeds by computing iteratively the value of $x_\\la$ at\n", "the breakpoints $(\\la_k)_k$ for increasing values of $k>1$.\n", "\n", "\n", "We suppose we are at iteration $k$ and have compute the solution $x_\\la$\n", "at $\\la=\\la_k$. For small values of $\\ga>0$, the solution $x_{\\la-\\ga}$\n", "is supported on the fixed support $I$ that has already been computed.\n", "\n", "\n", "Compute $J=I^c$ the complementary support."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["J = setdiff(1:N, I);"]}, {"source": ["Compute the correlations $c=\\Phi^*(y-\\Phi x_\\la)$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["c = Phi'*(y-Phi*x);"]}, {"source": ["If $c=0$, the algorithm should stop.\n", "In practice, one can test if |norm(c)<1e-6|.\n", "\n", "\n", "We can display the correlations $c$ to check that conditions C1 and C2 are enforced\n", "for $x_\\la$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAR/klEQVR4nO3d25ar\nKgIFUOmx//+X7Qe70qlomZsKC+Z8OKNOzK5CBZYgMWWe5wkA0vyndgEA4BMCDIBIAgyASAIMgEgC\nDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIM\ngEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEidB1gppXYR\nADhF5wEGQK8EGACRBBgAkQQYAJE6DDALNwBG8K92AY4kugDG0dUIbJ7neZ5rlwKAK3Q1Atu0HpY9\nhNzmuM17vMd7vGfw97Q/HijtF/Fdpfx/p+5/BmjQEhk6qg90NYUIEEd0fUyAARCpwwAzZwgwgg4D\nDIARCDAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkw\nACIJMIBOlVK7BOcSYAA96j29JgEGQCgBBtCjea5dgtMJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJ\nMAAi/atdgLeVn4+Xz1ufcli2bm4CoCdhAVZKuYXT/c8Pr6w3AdAZU4gARAobgT21nkIsv59oaWQG\n0Id+Auyv2UWJBdAlU4gARAobgc3zvF6FuIy3NjcB0KuwAJu2wslsIcCATCECEEmAARBJgL2plOn3\nunwAqhBgAEQSYO+zVASgAQIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgC\nDIBIAgyASAIMGIMvQuqOAIN++fq6m+U4OBp9EWDQKZ31veVbkHwXUl8EGHRKl03vBBgAkQQYAJEE\nGACRBBgAkQQYAJH+1S7A28rP4uB5tbxqZxMAnQkLsFLKLZzuf55+0mt55WETAP0JC7Cn7mPs/pWb\nsYJt2fehdhkYRlcBtjk4GyuxAIZhEceFrn8wnfAG+iXAAIgUNoU4z/N6qeEyYbi5qS3z7PmqAEcJ\nC7BpK5zc7gIYkClEACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCbA0nuUBME2TAAuzpJcM\nAxBgkTwxC0CAhRFdXOz67wCClwkwACIJMOBvBv00TIABEEmAARBJgAEQSYABEEmAARBJgAEQSYAB\nEEmAARBJgAEQSYABEEmAcT4PhAVO8K92Ad5WfrrC+Y+ntJVS/trEMW5p5DgD9YQF2H04bQZVcaV/\nmdfTa56NwIDDdTWFaOx1EQcZaEDYCOwDD2MyCQfQh34CbAmq239vQSWxALrUT4Dt3xsLtowge9oj\ngCOEBdg8z+tViP0klqwCeFlYgE1bU4IPr3QSZvf62yPOVopqQ/e6WoUYT4/DIZahvI8u0DsBBt1Z\nroRcDzXF82hOIMC+ttRLVZOnVJJczl2TBBhcwrRerkPOnQHxCQTY10zX8Ao1JJdz1yoBRgPGmYPV\nFeb64NwNUqvrEWCXGKeDBhYmjc8nwGjAx+MSvQPNMto+nwC7itp8OOnFsc6YKdHwzyTAgIZddpni\neiiQACPWzrXtK5fSOqwzHHtUr7yNZDlxIAHGkNxgP8+BR1WcsEuAfcHawmY97fg6vtzur052eZo4\nggD7Wn/9xXmOOlYuHf7SwmGRN1xFgD2z86jDjq/iz9BC3wp0RIC9Rko1xen4iyPDSATYM2f0CGOO\nRfStwKEE2OXGTC+AowkwACIJsMuZSQM4ggADIJIAAyCSAANINvC6MAHGazz8Aho0dqsUYDB6LwCh\nBBhM0yTDyDT2quZ/tQvwtvLT0cyrM7eziQPUOqqlDN5KadrS7aiiNYQFWCnlFk73P98sr2xuogm3\ngc6LJ+iygZEKwwcM3KvqagpRaAGX8pUUVYWNwF7xMPwqvy+RhFxl8/zelODy/lAml+BMXQXYklUP\nESWxALrU1RTi9H1c+bQTB3LxdCzN8ySxBzZsBDbP83qp4TJnuLz+1ULEzFP4NvNa9GSzPqvkYwgL\nsGkrmZZXDpgqfPf2zCte/IXa25Wsyw/14g3RQa5Ej5XZInqbQmzL0pAabE6ZlfUYX56OBs8mD6wM\nHIYAO9PrDUlj27fM0VcPj2uuSFrY0yAjtJ1blVAxfhNgDOaQqWa4ngxbEWCZRrtIb2RS6LMCjHay\n0jV4vu7r/1+VsLUyX0KABRqypsZ7OGs7vWT1nCbOqIMzARao2Q5uvPbzkkaGj7wu7mSNWscEWKzW\nKqv0ektrpy9IIyt6aIAAy6cxA0MSYOcYMFQMKbjGldNlAzbkKAIsn+SAU8mwVg0ZYBdUR6ECHRh1\ncUSKvGchfuu23lSlnIZ/BqMra0bTV5MfbwTmkmpt2H5cNYBk4wUY98T5NPzu5xr2wusbfdV2AQYE\nGvXZExcJ+bCdADtO3fP9UOESKt8THewC52lzJKHSXqvzAJun167Rvq92B9ba6tc+tf569R0nTpsx\n1oGQmwudB9h7vuk6q5/vdQHqVr7vo6j5xrNNAO/r+xrlaaXte/cvN3CA3dek6vFzr6nCfOazJtrN\njuuh2KeGHGSAAIvuEKs46isfBzzyA+7yuwY/RMM2jXMMEGA79quRwT6f0T21b/DW3cvujx1gcJJe\nOghomQD7m+to6NU4rXtzcW8vuy/A4By99BHQLAEGNMO86+G6XjbS1dPoy0/tnzs9W1DHwwdOXvwn\n7zbD9tOrr0e5d6CfACul3HLr/mfgau1HEV0whQg8c800VMRFZ0QhvxG1gLafkcr2COxnVnGeprL1\n37c2/e+XH/cLX9l00q/tdVNThWlkU1OFOWnTrW2GHq7pjzfUPZLtp3XnAXafZNM8b8/Lr1/cfGVx\ne31niv/F2f+dIm3+nlfKuf/izi/cecNn/2qnMA+bXtyvF/duZxdeOS/r92+Wc3PTW+f948K/WMLN\nv/j03z49XE//1dPj+e6mv47b+g9N00u/6pW/tf8XN//Wu63g3TO+WbDPNu28+MHNzqpMIb6m65U8\nEE/bPMo8Bx3MfhZxzPP8ZBXilxO7OSc1xgcL1TjP0kBy7n88uqwuqbTN6CfApr9y64FOswW5vWT3\nlskrSDDSFKJpwHYs0xQ9nYujhi8Vw+P6M5I+5qO2kQJskl6cTAWDC3U1hQjVHBJdo+Xf/Tq9oRh6\nHmSwERhPaVRwtvvbGVrcFwQYwOVGG22fwxQidxIb1eGzMZu/0GXysRzPRWKLa4kRGF14+oyGd3vM\n6tM7I3TxI+xjy/JvxRmBEe7pNey7H2za/IWulI/leHIEAcYAPhtO1R2E6eKPkjzCOFd+HTOF+BFN\nAqA2IzCga/njDP5iBPaRzh6DBNcwdcGhBBhAA6qne/UCvE+AAZcYed4iMBsiuAcGDXv4pI5+sGPD\npvsXjMCgtp3Pk7b/HUAytRuBQ2QB1gBdADvu+5SLu5jXa6Y6vCMwGFIIsKo0e6ZLhlnn1bT+vpuU\nHO6BVaXZs2iwJjRYJPjNCIwXGClGM0Ji0V1DFmD8If9J1e1yVOEIAoxnTr1+15V/rKkrjLO/RO3U\n3WzqSJ6nx4G4AGNXdzW+CYcc1fZX2C+6DwbqsYiDPzztGQ+5br24/y2lfo9/YAGq78vN/peofZ9h\np+7pUb+8hdo1GCMwPtVsW3XJ34FBpvX4jhEYX6j7lY/vajZx35VywL+xfI920CkLKmpH8gKs/LTe\neVVjdjYxCqc+1DquPjiVt3Hb99VghKuEfGEBVkq5hdP9zzfLK5uboBPL6KQWPTvNCAuwfUKLUXxc\n1b+5t3Tq86ia+j3H/ipO0+EijofhVyllmVq8/QB8ImXhPsNoegT2kDdPB1jL+x/eNv8sNDA+q6Pi\nrfgD74h048ul7Y4kLWk6wD6IHClVx5UhcX0giUBoUtMBtjbP83qp4TJneJsnfNhKTXXPgjpAOnc9\ndoUF2LSVTMsrEquaK4/89WdZvYJW5QUYwChcP+3qcBUiXMcMTyJnrRdjBJj6CtCdMQLsAymZt1nO\ndwvf4INTmyrMX77/gqWI3exMj1+LNSwB1pcOOsT2d6HBvIchjbGI490LrpQLtHU5P3tKXlP7G/cY\n8m8MsptwjjECbCgd9ImN70LjxYNhmEIEIJIAgxDuusFvAgyaJ7pgiwCD5vkeE9giwCCB9IIVATYM\nn14C+iLAKrk+TlzCA33xObCRyDCgIwKsElkC8B1TiBzKPTbgKgIMIJCLRVOIHMzUKHAVAQYQyMWi\nKUSARpkkfEaAQb5vejq9JLEEGLAi1VpgkvAZ98Ag3LHdnOhqhPR6gREYDGyeHztKT74nR16AlR87\nb7iyPOs/X/Ovw/fOSC/tghOETSGWUuaf1nX/8/0bLi8U8DdNktPkjcB2bEba1dZzMjAyc5KcJmwE\n9oGHMVn9hIPRaHSco+kAeyt7ljff/nt7s8Ta1uXETin6ShhH0wH2Vvbs3xujc4fkcZehzklu30mr\nt6mn6QBbm+f5NiyTWF/p7Ih1tjukUPGqCguwaWtY9vCKMONDag6vU1sa0NUqRIAhmO6epkmAxVBf\nAX7Lm0IEGJrZyx/DB1jKyEaVBfjNFCIAkYYPMCMb4Fi3j4hxsuEDbBogwzQnuF6tjmWkli7ABnDg\n01RHahs9e7imcVqP5YneVxFgY9CcuKc+9Gqw7Bx+FSJvGaltdO7+VDqtZDICAyCSAAMgkgADIJIA\nAyCSAAPoyzCfixBgQEt87p6XWUYPtMfK/m8Mc/QEGNCSYTpfvmcKEYBIAgyASAIMgEgCDBiVtY7h\nBFirLCYG2CXAWmUtFpxNKwtnGX3DtC44m1aWzAgMuvB0tnnA6egBd3kweSOw8lMp59Wl084mADoT\nFmCllFs43f88/aTX8srDJsBcGf0JC7Cn7mMMBqLOrzkmvesqwDYHZ+X3PLhsA+hD0wF2SPZILIAu\nNR1gsgeAvzQdYGvzPK+XGi4ThpubAOhVWIBNW+F0e0VuAYzDB5kBiCTAAIgkwACIJMAAiCTAAIgk\nwACIJMAAmuGr2N8hwACIlPdBZoBueRrDO4zAgI+Y5qI2AQZAJFOIwEdMdlGbERgAkQQYAJEEGACR\nBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYfMczbaESAQZAJAEGQCRPo4fveCg7VJIXYOXn\nlsO81XEsWzc3AdCTsAArpdzC6f7nh1fWmwDoTFiAfaD8XiQm2AD60FuAracQJRYcrxQ3/6iu6QB7\na/C0P7sIQGeaDjAhBC3SMGlD0wG2Ns/zehXiMt7a3ARArzqfajOXCNArT+IAIJIAAyCSAAMgkgAD\nIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMg\nkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIFKHAVZKqV0EAE73r3YBjiS6AMbR1Qhsnud5nmuX\nAoArdDUC27Qelj2E3Oa4zXu8x3u8Z/D3tD8eKO0XcdPDmbjfi1L+v1P3PwPQk9QRmFgCGFxX98AA\nGEeHAdbU4KyRhZGKcU8xmirDpBi/KcbrOgwwAEYgwACIJMAAiNT5KvOIaVyABrWfDp0HGAC9MoUI\nQCQBBkAkAQZAJAEGQKTUZyG27Lb08f6Zwg+vXFmY5Y/WKsPD3617KOoWY+dcXFme1oqx/t/LirHz\nR+seiouL8bQmVGy2+wTYKW5nfZ7nik/Hv1W7WmVYN8WKh6JiMe4/zrE+F5ednb++Wqh6MfaPzwVl\nqFVRNw9FlYq6UxOuLMa7TCEer5ET3E5VK6W0U5gqGvmq1Ydi1CrSQzGq1I3NM3J9RW2nYtQuwoeM\nwM4yeJd9s56UqOV2LVm9JE1RURcqamJNEGDHexh0VyzD/X9HtnmrgwYralwHeqxaFbWFmvAZAXaK\n6lVhfZelbnloU2sVtW5hRhZ68AXYwR4GPcscd/U1PLXKsP67Ixfjr8JUOTttVtSp0tFopIZUKcYr\nNaGFirHJVQ8AkaxCBCCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgAD\nIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMg\nkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyDSfwEa9M0sSOoEgAAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "stem( I, c(I)/lambda, 'b.' );\n", "stem( J, c(J)/lambda, 'r.' );\n", "plot([1 N], [1 1], 'k--');\n", "plot([1 N],-[1 1], 'k--');\n", "axis([1 N -1.05 1.05]);"]}, {"source": ["The mapping $\\ga \\mapsto x_{\\la-\\ga} \\in \\RR^N$ is\n", "affine as long as long as $\\text{sign}(x_{\\la-\\ga}) $ is constant equal\n", "to $s$ defined as\n", "$$ s_I = \\text{sign}(c_I) \\qandq s_J=0 $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 13, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sI = sign(c(I));"]}, {"source": ["Since the sign of the solution is constant equal to $s_I$,\n", "under the condition that $\\Phi_I^* \\Phi_I$ is invertible,\n", "the equation C1 shows that the solution can be expressed as\n", "$$ x_{\\la-\\ga,I} = \\Phi_I^+y - (\\la-\\ga)(\\Phi_I^* \\Phi_I)^{-1}s_I $$\n", "and hence one can compute explicitly $x_{\\la-\\ga}$ from the knowledge\n", "of $x_{\\la}$\n", "$$ x_{\\la-\\ga} = x_{\\la} + \\ga d \\qwhereq\n", "  d_I = (\\Phi_I^* \\Phi_I)^{-1}s_I, \\qandq d_J=0. $$\n", "\n", "\n", "Compute the update direction $d$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 14, "cell_type": "code", "language": "python", "metadata": {}, "input": ["d = zeros(N,1);\n", "d(I) = (Phi(:,I)'*Phi(:,I)) \\ sI;"]}, {"source": ["A useful vector $v = \\Phi_I d_I$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 15, "cell_type": "code", "language": "python", "metadata": {}, "input": ["v = Phi(:,I)*d(I);"]}, {"source": ["This affine expression of $x_{\\la-\\ga}$ is valid for $\\la \\leq \\la-\\ga\n", "\\leq \\la-\\ga^\\star$, where $\\ga^\\star>0$ is the smallest $\\ga$ such\n", "that one of the following three even happens:\n", "\n", "\n", "*Event 1.* condition C2 is violated because of saturation at +1:\n", "$$ \\exists j \\in J, \\quad \\dotp{\\phi_j}{ y-\\Phi_I x_{\\la-\\ga,I} } = +1 $$\n", "This condition is equivalent to\n", "$$ \\exists j \\in J, \\quad \\ga = \\frac{\\la - c_j}{1-\\dotp{\\phi_j}{v}}. $$\n", "\n", "\n", "Compute $\\ga_1$, the minimum $\\ga>0$ so that event 1 happens."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 16, "cell_type": "code", "language": "python", "metadata": {}, "input": ["w = ( lambda-c(J) ) ./ ( 1 - Phi(:,J)'*v );\n", "gamma1 = min(w(w>0));\n", "if not(isempty(gamma1))\n", "    i1 = J( w==gamma1 );\n", "end"]}, {"source": ["*Event 2.* condition C2 is violated because of saturation at -1:\n", "$$ \\exists j \\in J, \\quad \\dotp{\\phi_j}{ y-\\Phi_I x_{\\la-\\ga,I} } = -1 $$\n", "This condition is equivalent to\n", "$$ \\exists j \\in J, \\quad \\ga = \\frac{\\la + c_j}{1+\\dotp{\\phi_j}{v}}. $$\n", "\n", "\n", "\n", "Compute $\\ga_2$, the minimum $\\ga>0$ so that event 2 happens."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 17, "cell_type": "code", "language": "python", "metadata": {}, "input": ["w = ( lambda+c(J) ) ./ ( 1 + Phi(:,J)'*v );\n", "gamma2 = min(w(w>0));\n", "if not(isempty(gamma2))\n", "    i2 = J( w==gamma2 );\n", "end"]}, {"source": ["*Event 3.* $ \\text{sign}(x_{\\la-\\ga,I})$ changes because one coordinate\n", "vanishes:\n", "$$ \\exists i \\in I, \\quad x_{\\la-\\ga,I} = 0 $$\n", "This condition is equivalent to\n", "$$ \\exists i \\in I, \\quad \\ga = -\\frac{x_{\\la,i}}{d_i} $$\n", "\n", "\n", "\n", "Compute $\\ga_3$, the minimum $\\ga>0$ so that event 3 happens."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 18, "cell_type": "code", "language": "python", "metadata": {}, "input": ["w = -x(I)./d(I);\n", "gamma3 = min(w(w>0));\n", "if not(isempty(gamma3))\n", "    i3 = I( w==gamma3 );\n", "end"]}, {"source": ["Compute the minimum\n", "$$ \\ga^\\star = \\min(\\ga_1,\\ga_2,\\ga_3) >0$$ so that any of these three events happens:"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 19, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gamma = min([gamma1 gamma2 gamma3]);"]}, {"source": ["Update the solution\n", "$$ x_{\\la-\\ga} = x_{\\la} + \\ga d. $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 20, "cell_type": "code", "language": "python", "metadata": {}, "input": ["x = x + gamma*d;"]}, {"source": ["Update the value of $\\la$, to define the next breakpoint location\n", "$$ \\la_{k+1} = \\la_{k} - \\ga $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 21, "cell_type": "code", "language": "python", "metadata": {}, "input": ["lambda = lambda - gamma;"]}, {"source": ["Update the support $I$ by either adding or removing an index."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 22, "cell_type": "code", "language": "python", "metadata": {}, "input": ["if gamma==gamma1\n", "    I = [I i1];\n", "elseif gamma==gamma2\n", "    I = [I i2];\n", "elseif gamma==gamma3\n", "    I(I==i3) = [];\n", "end"]}, {"source": ["The Full Algorithm\n", "------------------\n", "We are now ready to implement the full algorithm."], "metadata": {}, "cell_type": "markdown"}, {"source": ["__Exercise 1__\n", "\n", "Implement the full homotopy algorithm.\n", "Record in |X(:,i)|, |Lambda(i)|, |Sparsity(i)| respectively\n", "the solution $x_\\lambda$, the value of $\\lambda$\n", "and the size of the support $I$ at step |i| of the algorithm.\n", "or bookeeping\n", "nitialization"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3d65Ks\nqBKAUTwx7//Knh92u20viJBAXr4VExO9q6stNVMSEa1lXdcEAIA1/5u9AgAA1KCAAQBMooABAEyi\ngAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooAB\nAEyigAEATKKAAQBM+m/2CoSzLMv2w/5d2MuyVHwvdt1fYYxTlCuCdc0TaEOUp6OADXVM8eoKtCc9\ndBKJcmpoFjEAUdaAAjbNnrjpb006vn58JfMGqLWu6x7i/edUHGWYQJRnoYANdU3u23Tff7W9Tu/M\nllNMt3/uL36KMqFXSyrK9ERbUMBGI5UjaI/ysR2ETiLHMgMqLZiFOFQmTddfI9cHPUhFmWTQrD3K\nFK12nIGN9jTv6HZ24u07od/TXNPCKJ+ujJIAOjVG+fad+IQRdqW4+BEBUY6AKPfDECIAwCS6BgAA\nkzgDAwCYRAEDAJhEAQMAmEQBAwCY5LyAcatgBEQ5AqKMK4cFjER3jxBHQJTxytWTOMh49whxBEQZ\nhVwVsNvHYl4PhpLHbvKe3zeklEqXM+CewqcnnxLllvc8RTml7Unqf5ZDlI2+R1uURTi8kVnqi+aQ\nfpI+pUN+a3AKK1FutEX5ugufXh+DKMvSGeVGDq+BQdCW1kaTG4We4kvcPXEZZQoYAMAkChgAwCSH\nBYyBcvcIcQREGa8cFjAAUTDhPjYKGACbtupFDQuMAoZaNByYizHG8ChgqELnF0pQxgKjgKEBbQeA\neShgqELpAjAbBQwAzFqWyCP5FDA4FfvABiKggAHAcFK9q9iD+RQwAJiBEYJmFDA4FbtninvUDF9C\nFjCujgABabt5kT5Ws5AFDAhCT2OtAQXDnZAFjDxGBNpOOJTg8HckZAEDNJMd4qa9hl8UMMApShe8\no4ABKSWm9rhAEIOhgEXFca4WZ06CKGmuUcBMkToaubx/i8ph3VMESXWnKGCIYav9NGTRbCWNrolT\nFDBTpI5DjmpAGzpY31HAEAM1G4UoJHZQwAqQ0EAQHOmmUMC4OgLg19xzdEYIPqKAFSQNWeUVvRZs\nTpnAIW8EBewXKRsN9xJgQyaYRQFDVIa6LAxxF6rbS0zwMeu/2Svw2fKbo+sl4ZZD+l5/C0Oeoiwf\nYvJkHuEoU+PjMVbAlmXZE/r4806sbm0HA63bDPkoR+yarKu/1lk+yh73koxl8dqUeRtCXJZlCZXE\n8QaXwoU4JJkoqz06Rq6V6yt83grYuq7ruh5Tv/5IcNpnse4a4jSsqqltEN15ivKs9THMdTtmbAgx\n73bYwfmIU7Bhk6dodowyg8nDCUR5j1rL0XEMvWwajD9snSawnzMwemfuEeIIiHKrSDvQ2BnYcVTh\ndAX49lewKBPlOSHeP2tAxzlM63ONpuSBLBWp6hXYP31wW7Rf8YrRBhorYOkup/dXZtYtBppEPUWZ\nroknTzMPNUZZyVj9a2U6raeGde7JzxDiBFzSN2SP1P7oS52xG99h17kfPlG4Ffmbo+vWtnw+ocIe\nQB8UMCFhMgYClDS4GtahnbmtGLnC3tslCthfX/vm3vPDvGs0pz83aGKDe9wbsy7SiOsR0GuMBO+m\narmw1/LnHlHADvR35Qavof4dMkvLKZSSZmj6Cthy3V1KTqOfKF89IYEL2FPffK5M2nVNx+vnms7+\n48prCKsenvbG+BTdpkiYPjR8CVzAnrQf4Z3y21PT09X4Z+doaNQ0rMNI17HQ668Kl8OTeswKXMD2\n/JM97K/nMRabFYvrvOv9VI7TVObrrOWWoJf84dcm22gS5s2qH2rr1inKatdTVOACNlimw3ginnn5\nYcntP3/pfrtF/drx/f6b8kDXuT2/LLmu5q+Gba5b/SmZ/WV+JIELWKdZWJl2cy8k/Zq2r33tzLZz\nYF/d1onrPwdM06hoowloD3rObkXiq2dzysQrYF0j1HUeRMvSXu/ex5Vgb0AQwRovnwadImKqkMwS\nr4A96Zouexe4a1+YjrYqNEA+jJ8TNP7jjkw1ICELWL8I5QeXZBeO3tTOR22ndsVQiAimlIIWMCCy\nuomOnlSMDA/uPpbMx2mMkeYnghajgNVqCbzppLGf9GJemwDTJ81eozxlSFCWSE01nZy/KGBVSh6W\n8XprjkIUpwjqWi5VudG+Mr2bbz376uq494yXMQpYlUzUj/27r0msJOmnzLmySOeukKo0OreuXeF2\nVRy8+bstpajqScwWqYCNCfzx7P7TXVaDH3V4q25oQtUR1bgyJX/efu2hk+PK93gUiKqS9mll9ORn\nb+WHuQv2vpFZ2NdBv8bY51uH3odZnMM443qvXqeHip0+cfuUkSGo+Kyb8QODjV35LBXBpjxz/HY6\ntBs3YY+y5XoW6QxsfJxK+vJPXX7ZOUIiM8Kf1kfVAdDyBCmp+WblJzGdTl5vN8R996Xius6p/yq1\niyZ2U1QdjP1FKmC3jvF+PeDHJEfXEbn2ps3BJK6Tln5Giaea2q90Sf25odbw9cJtpqdY/s/yTxym\n8esxjQtfwPIKxyIam7bbTNL2WA09a3KrscbMaox6PLgZV5mRgz0ET0dcSdm7XXKS6KZkltD7Wxcs\nCF/Ael+xr9AjLwU3U1tlffX1sva1yei0yY6f9NHJlNY80780pDxVTA2xxChgUuM/Ji4C3er3vFEl\nif56LitLw0E+ZR30RPyTuqtigz+6x58fPW/amg4DGKZ6qDEKGD6pO4ZtNW2vlzrym6Pn8B5wGrfv\nionxXZb18HPlQvRETcTXsc2nzTd02F44L2BrOkw2LYzTtWk7zW4qaamn58SYFdDWItyGpmJXtExl\nHGZI3/zmE7UFfdMvOl+X3Kkn13tiV9Ia2SznBeyfrq3PMWU1nL6I9Jend7qL/eumCC9X+mbzKSep\nT594fF3/18Vda2fhBCsfyg/G8k62C2EKmKzXR2zMPebLH5YjMkpWcf+NoBgHas5ra+VsF2koqBk9\nVk/5Js/jqoAtv/68Oib2ei7Vyi6z5GsdBroJ8chjW2GUO/u5vD/W/YFcweAO/6O9Kfu0Dw12dPwU\nsGVZ1l8CqW/LcXtPPxfesHL923KjmolciK03VY3ycSyZSatGLsq61zwlud1bceHtJMYR4aeAfaD/\nMOgn8rZPJ3JVUmRRn/R+ViSuKq5ESEXHVJSX1Uuh3jpu559NBUOnpfx5rp1z6T7EiSiPRZQjMFIX\nnJ+B/Un565Nj9leeHifzqR90+1fXKarHjzst/PSrzEc8vZifXXL7zqf1/P15vd2E/B+OdRPl4z+f\nopweNuR2aV//NhPKMVF+Ws/MFt2+ri3Kt5tz/GdJpF7Dusvnxm0o99crovz6hmscr+v5+tHXhauJ\n8lfOC9g9O+GRl2/7Tuz2eUuOYUHlt5TKrlXkTBYx/sFUEOXn+8COl3zP46Jr/6/aav+g267c6xSM\nAdulpozlQjx2PUrvyOm966bMF+hMb5RvA6ps7zUxuC1+Cliam+5PXm9wziRN75uvX5U0wWOTPhdi\n2YLxaWkrF1wljTuQh3VtWzRml+vkDDmEuBnzRDWDw8p/DB6L86fr3hsWHafNX5P9GXUbkUDMPdYM\nHuwhC9jgozFzMT/z2/yvWj73q9dhTGtJ74dIJhO+Ouw3BUIWsJI5Qu6VbP6nxx/r7KSPXyuLeVUy\nac3idvkg1U3xGMGQBSxDZyucd3t/68hk1XlgWAlly95r3/PlS9DcTdFP5zFyZS2+FDDd8nnf9fkx\n1zuKMuugsH+XOc9WdZT6ntiGNPBaQLyHcVDADgpb4a4T8z6dQl0LTMWcxuq1tZPlZ9Nrbf5hxD1W\n7zVYhaOI03edLLs5nNGyUdaiTAETNWBOfO/c+jp33E6u/yhc4TEjZk/7MH9T+chm1+jIoYa07LHT\nSjoiT9tuLogFKGDfZTopJf2XfpP6Tp9+6uY/3Q7ScpuIhmZivHybXniD3eBdp2RGqzYTt6v86S3V\nSzsxdWpViALWpjwh9r62eA41jg5dHzfgLstr5LspaWzbVzjMmFmla9EtbCuV58PEs4qWkyFBx2w8\nfZzmwAmhgLV5bREyT+L4Olj31e3x8/RKgFxv8vr1Y8Na0lNYy+uQYxoGxwrPgW5/vn2R5zQWoIDZ\n97WjVz1rQEMzMV1Js1JRYF6VTLI4uh0r3t+j/LyqnHj36/XZb6eP/rTY127Qq05RM3toU8DaHOea\n57tU1zfMbUTMpuxMXeflC55ntzCXGCKb/3Wr63qBr+9pr3DBOC9gP/GvnsF1228qnEa/LJWT8tXO\n+wp7OJXP7CrvpnSNspsTrAGe4nWdD5UJ1mt6FE6VeurvDpgNazNhnBcwk8ZkktmUvZq5GSXzTm3x\ntC1Hr7N/X1XvGa5TduPq61SulpTWY9+n+lpRxf3Fnz7laehA4XmYWsducuHp7/XPj//8NHvzq2Gt\nFSl09HVvXNuN9fv37Jz+WRGRlntdXOMMzAISN+tn7/TuIDd6neOOAQqHgjPUjvCH5PwM7JuKbntv\nnIoVOl48KPd0iaLrgNJIjnrurbs7H6+n86rbyRrP+1NfC+IcZ2DNXLQOQY1sbcgTQT1Og8qTwcFV\nTy/ZSAH79ZqUXa+I2GJ6YxU+eUTDOqh3k3PstPAoYM0CPsnC6MZOrLuUKBGyN2Jbj0jLPHvr2/6L\nAiahdzaIfxNKy2rYTf2Srkb7RX7MVRKpbtFc704U0U/sAmalVaqbmCC+DnZL187oJlhJVD3G77GJ\nhdPHsVkldgEz5Ok5VZwxjBGtjbC7sdfz7OuMjx5bt65pXZepd9UHFLuAfWqVqic+tTyHvkLjsyGo\nfLNkHhN8jWa/vrwDt7cP6980Dr3vYhewCoKHgdF8Nbra+ulvYe3SsG/1nMQrWQ0J3Mh8p/CWxkK3\n90Lq5PJWgcxqC966LvVA+vLn5hkNB56obRYUi3QGxgFfKP+VV74PM0NJ4jsQgiq+gEJqaeV/PiWa\n9lPI3hnY8psB62XvL4fkuP72g9Pf3j4fVkNLN/7ZV6MeTfQU5ZoQP71NdtcNu0wlQsdT0ySj7JLm\nFNLBWAFblmVP6OPPu8d0b6xnoTJp9vbmo9y3RXPfXBZGtn+Fk4/y7TprC6ijp1NqYKyAvdr6bvJt\nnM4Do9861D35dIheIa7+HpYnGjLkaotdZt10tK29ovy79PvXXz9uyoAHnnkrYFvGnzp3yd9AhIbN\nmTQMdQ3xzz9F1krDju0tfxKgYw88RnnMgDCMUF3Alr/H2GsRun2Dt9L1ifhZRQefovz02zXfKJe4\nXbKG0xHd4SskGWW8CvOlZaoL2Kd8zfXOwuo61ie0t8uj9h5iEqDc674SzRzJKLew/jzPCi46QE9U\nF7CrdV2vM5e2jL/9VSVPIe90taPn/slEWSzEJhRuo82O9jWaRQeyhUGF+fadYzM3yhkrYOkup/dX\nQjRq4lS2BU9RHh1ifXvmkXgct9P3nunxNPMwF+XXSSiDKa8QenZUH5FuZDZtP07aR7f3v41wY3IE\nneKotmlunKRDwjti7wysO835fao9jYsKezD7G4bqsSFudk5Xw/aSbH/CyywPCpgRgscJDdPsO7UB\niGAIESEpP/sUrK9hS7Wtk4zX9RS/wKn8EChDAQM8stJw92PrEm/JelrZloEoYIAagtc4dy462pWs\nbPjE9bSyix5QwFCGHv1gjS2L8YZptLDpbWug9YIC5ojZLMQP2bMlapgJGsKkYR2qMAvRhd6ly2x+\nA+9mpbfIYRv76dWcgblg63o1AEjgDAwAZhDpccbutnIGBuAXl1FtiV29EgUMQEqULuOiho8hRAWi\nJh8UCd+Xh0WcgQGm0N2xrkcEo96uzhmYAiEzDwAaUcBcoFceBH0d4IACBgAD0QuRQwFzgUMCwCt3\nQzVM4gCGcNd2wCpHqUgBA4AY3D1zjgIGDCEy0dlR3xlzOKpeiQIGADCKSRyAEb76zkA7zsAAACZR\nwAAAJlHAAAAmOSxgCzO1vCPEEUSMcsBNbuNqEkfEjA+GEEcQN8rrSg37xNUZ2LquKzO1XCPEEZiO\ncmsBivrFKHVcnYHduvbmTsfGbXeP91S8Z2KjQ5SHvYcoZ96Tkq71sRjlTxYrK3pyisRxK5bl30Yd\nf4Y5T1E+hZUo21V4IF//qdOycPo0lNUhxPWv8SvQdZi+38JtLdlxlMmfzfQQJ9GNOm0BUe4t4hCi\n/iV3Xbih7CxncXcZzZ+JLMaCKHdl4KwcAIArq0OIAIDgKGAAAJMoYAAAkyhgAACT/M9CFHSc87NN\nftlfaZwLs9/jcl1g40ec7orbX1/XtWXJp78VX+2JiPLTKjmOsniIU8ExUr1wwRR9janyEFPAvhG/\nY/qYi9cFtnxE/q779pU/5rfgamtAlI/LeVpJN1EWD/F1mU8fV7fwJJqimSg3rvYADCF+syyL7L0X\n/e7fvF2y1PorTGVBRDl5D3GSjnLXG7GvC+dA3nAG9s2xtzJ7XWpchyBa7L02kaXpQZR3dnfCK9lj\nYTDZFLUbZQrYB0ZjvBNc/9PwgifWN0pq/R2HOBnfLg7kHUOIpYz21Hbi62836TOI8pHLECfjUeZA\nPrJ65jhFv1lYg+enSc1c2v/c5fy0yFG+hrjHak/Uaa7d4FmIjYEoibLyEFPAAAAmMYQIADCJAgYA\nMIkCBgAwiQIGADCJAgYAMIkCBgAwiQIGADCJAgYAMIkCBgAwiQIGADCJAgYAMIkCBgAwiQIGADCJ\nAgYAMIkCBgAwiQIGADCJAgYAMIkCBgAwiQIGADCJAgYAMIkCBgAwiQIGADCJAgYAMIkCBgAwiQIG\nADCJAgYAMIkCBgAwiQIGADCJAgYAMOm/2SsQzrIs2w/ruu6v7D9/Wk7FX2GMU5QrgrUvIR1SBaoI\nRpkQ16GADXVM8ZYKdGzdoI1UlGnUNGuP8nYUVxc/JArYRHviprvudqYDvuU6NcyEPVKnqL1GOf1t\n4KAZUZ6FAjbUNblv033/1fY6vTNbTjHd/rm/WB5l+uaaiUSZY7wRBWy0TMqWnFQdz9hIerUao5zo\nklvQHmU0YhbiUJm0Xn/ll3B8G22cTu1RpvnTrz3KaMcZ2GhP845uZyfevhP6Pc01LYzy7TuhDVGe\njjEopRgejIAoR0CU+2EIEQBgEl0DAIBJnIEBAEyigAEATKKAAQBMooABAEyigAEATHJewHiiQQRE\nOQKijCuHBYxEd48QR0CU8crVo6TIePcIcQREGYVcFbDrt+9c/5meH0LIe67v2Z56f31PSuvvm/8t\nZ8BN8bchvn1F817V9p58lPc3Kozy7RreLi27kBDv0RZlEQ6fxCH1fbhIP0mf0iG/j6/P2rWnsN5H\nee4qmvK0q/RHmSCX0xnlRg6vgb1blnTXW8HVltbX5Daa7rj1FE0TUTaxkhqYjvKTkAUMAGBfyAJm\nusuBQkQZ8M5hAeOil3uEOAKijFcOCxgAIAIKGADAJAoYAMAkChgAwCQKWEjbnXDcDOcbIb7yt0Ni\nR5kCBiCGraFX0tzHLjxSKGCmSCX90wM2ItPToMhGGTuXae9scz6igBXQ07QBcCN27RFBATOFjO9K\nye5VshqAehSwAjQoAKAPBQwAYBIFDIBfXMB2jQIGwDtqmFMUsAdkvHvczR1BxdR5ssIOCtgdVTc8\nAgDuUMBSSpdaxbRDl26jTN/cE+4BD4YC9ny+RdPmiemzarIrAqL8HQXMeG+LpC/k8jFCOCG+wVDA\nfpH6eEWSKFfdmVPSCyTBPqKACZmVeWR8BES5RGP1UlLD8AUFDMhikNYKynw8FLC2jhtNmxV1kSK+\nEXB91CwKGPCMpg1QjALWhqYNiIDhFpUoYBQhAK75rb4UMDzwm/TCTO8o0yv/SnDr7HZzHceXAtaX\n79bBOkKDvH4ZMjL3XF/H/W/2Cny2/MZ+vYRkOaTF9bflH1D5h5kFSmVPmDb3KcoyIU6ie3JdJZc2\nPsTzmrbuUW6x3xxWtwL7JjitHEoYK2DLsuwJffx5Ny3dd1vi7r0eu9VLtl3+Ih/l+SFOf6OMKiOi\n3JLDXfN/3sHljLchxGVZlr+ZcX0l//fy6/RJZgVoLlNKDwH9EGUN47qZdSDKKaXnKJf//eTH0nca\nuCvZKA0ZPoq3Arau67qupyGI0t6cki9iCJN8da4hTp+iLLQS4z5LkJ2m7SnKs9ZHhZJHXhmJrxRX\nBaw1v11f7XxhpGkTaMKURLlH39xCBEuIRdmZktQ9Xrzwkg8ZfgrYh+EF5TQ0bSp3pp8Qi3P0oCz/\nUR5ZWryXMWOTOI6jCqcrwLe/cqVrIl6vKs/L+0yUnYc41e728ulCaqYPXKNZfyAzp+ZKTaC7MlbA\n0l1O7698a9RCJf3Iq4ASnqL8uW6FivInChq4p5mHNQeyISP3vIIod+VnCBHCaPS9chbZHhc1nxr9\n4+veR+dMCFzAehzGtwl9nDu0JX31MFFvSiY4yBLcnMw0sGNkadoceBpRJ7iaBC5gI5XMf4V+JbWQ\nKNvyVJBOsb6G3lk/z6bABayxoSlJ+qc3kPqzfAp6YdN2enFdie9Q1xh1fYIGwdXE3iSOVpnR7fLU\n/HqEyD5TCq+mnAZ5v2Cuy7KklASOqJKoZd7ALKGpAp+BtWjMV3/D6M42Z1M4YJjpEnnaLWo3Z/rg\nnqfqpTbKD0IWsDEJZy0V6uncTE/NCq62W8ckzsHKP04RnQfdcCELmAZKjoeWWZF4pSrK7ZRszmAK\nj46nSWHxokwB60Nh0qdua2Uq4yWpjXLmdg44cDsXLGR8KWDBnLKcWZHOZKZNEuU67DfFKGB9kPQR\nKIwytSqskEGngBUzca3odQ1DZvkHPqKMjDEhJkZDUMA0Ebm3uveRY6KJrzNmu4iyBvofjjOmymre\nAwXi3cjcQvnpS9fbpY0nepHjDlR+43nXKGve8E1LNu5RHlbDKvZqfq16PEXIpkhnYFLdDdluy3FR\nxxQXn+AutSjlrdvXnfb6/GVxjVHO7H/BDM9/0FwV25iZca52M5/IHsjmNv+vSAXsq2tzIH7GPbgP\n2MJFuhcRD8dr2ujZq8coq+2n60zFTJQrVrXr1p26UJZFKmCCafT0el02ZB4Oq+0oTeovHgjuMdko\nf/oIJdRG+av885dvudl21yIVsFciU/ikhoNkSR2ryq8MlSjZWNltVLLHPkV5o2TNK9yOGYqMPT5R\nsq8qomwZkziySi4GnN6zXRzWnBxSz5tRvInfvO6Qa4Ou/DQ0detqKL88tql7xtLTe6JF2VT3JfwZ\n2LHYDBgyel0Zux+huWw3rtiYr6sX/4geF/Ome0qz/IFcHsHpz7Yvdxvl8kbs8M41FTSDKoUpYIJX\np0reVlELZSccPh3nprKz0JqymywyMnx6W/4K5cgod17OTdM2JYWW5SfK/XzdLiVR7vQthkYaCucF\n7F/TloRm3z69TWqaX2/lQ0BvO+1P0/bz0pym7eUNUmulM8qZVRL8LA3nXruKGRnl7yk/5M0xUpM+\ncV7A/lAev/2QaG8oRXrKhX8+/UhuD+vITThGuZHg+ZDyQyPpW8P8+vSYCvspYeqya/qx/F2AAqYt\n9W9NT53vt7CM+y7BEmOi3Pgp19sKRWi++qjE4D6KYGd0X04aUsOsCVDAdl0HBwS/m0d25IekvyXY\nrNhS0aTq6QIe77D2NLXqVacLXdWfokakAqaZwdRRqm4K9RinKA+45zqO/PydYfsnO9mv1adFxUiJ\nkAVMT4u2k522tIuRxGZ0ivInpIQJ1WeZweLrqoAtv06vnt83IMbB0miY+xAnlZ2SE/1rqMZLlFtu\n6upN6hN1TujVx8+TOJZlWX+jfvz59q2D1smor/tn1P78EOKkL8r793fMVTENcuxqf4rympakajKR\nLFszbGdwdQZ2L+zJEJPTIhgQ4rqZ3C69PgTEAVPb8taNteO+42YqGDo9dXFvXu+cS4998zBRVnG6\nQZQjMFIXnJ+B/aT88ZbPdX384fTO6z/zbv/q9OLp435fX9L68/r1eUW3H/G6Dre/um5putst57+7\n/5P1aQOHq4zy9U+unv4k/1fpEuWf9byLcuZD73738oanOL7uhOun59N4uJco7/8sifJrWHf53LgN\n5doe5ec3XONY/tvMwvP/VMx5AVMhnw2a+5Uleax5/YcpaTJ0KlztaFHO92ZssZucBfxM4ljXdZ+2\nNG1cdP1yof71zSVP+cs/Cq/xrsbBd4y+URHiJB3lxjVJzQFSE9/N0CgLHj79Ai0SZaf8FLA0t1Gr\nsM8JNp2aY9ffWIiTgnZHJEBj21B7UfbBYKVkCPG7YUfX9MM4f5nn9DZnXG7ULddDTC/aN1xJcx/z\nIKWAzfSaUhqODad5Pw470DdD8TW0qsUoYKP0e4Kqx7yMgthFMCzK7V1eZZe9X1HARqlI4jF5f0rZ\nity1lvS6FEa5MRmuMap73gpRrjOyp6Lwa1e7oYCF156sdtJ9BJ1NVeHlzK5LwKvG3dt+0561C6Ku\nZiEKKGkRTHRCP2Xh8c11M5FMJb2fKH8SLcolukZZ8AtLuaDwgDOwv6aHv3oErx/fTfkUtVFe098/\nnL4hkOXvWOuMAnYw5mrELE/fMhPtmBkcPtHdu1Q/CpEoF6reUb1vV5+7DlpRwKqINIK32VYxBn38\nkx4ZbLRgN+q61fkoPyTGktbzGVhGvLashsoo55b5FNao4aaAySnMoa6p9nTAXI+lliu9kb2Gb7mM\n8lX0SOpcr3KVL5M02JTH6/TV53N3YNRnDlDA5HwdgaxLeruPBvZhTBPw8Ck/3XOi7EDdcwwyPdGQ\nQhawHne0HKd1Ze62iZ1tE/SO8jXiPTvjH4YQ8ck1T65HsXi/Id9Nuf8dCfBHyALW7wFoU+rTvjKZ\n8cN+n6v2iBJZsdepLmMfsvDzjWLYjUy/0wht10cTTB+TNCLqfWCfkuPTdy4cO+ktKVhy88d1xY7r\no7a0jCF7/B937GnKTN0H7eVf5BEb1atRuPz29eyq99yl8fmWHmoAAAerSURBVHMa090pYMualH+W\nKSHPwMpl+t1fM6k8RRofrzBsdsZa8LWzJig/ldy0d/9bPldtlMuT8Da+vQtD5qMb34mUEgXsRd3R\n+6lBvN6gKlgY1LY7qjxF+dM06IJw30e5JE8yqzGgmDnOorqZFJ+W/ylAT92U3lE2G2IK2JvqWX9f\nJiU+XrY1m1jGtO/nglvx7qNMiMd4qgFj9n9Jp1bDRHxrKGAdiHeXbOaWK7eXGz/FJX8WlX8DZnmK\n8kMpqn9UCqo4L2A/2fQ047lO3RJ8XC5SqUuUK1eFKPfSd7d+vfD8FX2UPpwXsJ+s7Jc9oRos5VeY\nZaOs6jkLIymP8qZxJRvPp1vyIU4iDeG8gKU0b/qWcibaqTJLIsr+naOcXOXwaF52XYACVq5wPphs\n4zg4k1o+Lk5V6FECZzUZLr/0q3Gmbo9tvP2+m2G8FKSvKGDFTBzYyNA/E0xQ1BbtAx+xrr4k7wIF\nbLZZ410Bx9kCTlMOGOVZs3jWNb1+GUqnzw0s6qOkrizmQcVRutY+ztzi/rmatRXtz4uqmLIfNsrR\nfA23oyhTwHQz8ZQjNCK+guy2zsofOKkSQ4h3tD0PV1taK9kt7dxsCDK8RllbszCDvTOw5Tcd10v8\nlkOmXn87Tb+Jf7Muq6TujcJTlLuEeHrPV2GUhxgaZQ0ymeZpMwcyVsCWZdkT+vjzTjLdCxfltX83\nTz7Kc1q01yjPesJLJ/138qAoi/e3Grs72qKsan2+8zaEuCzLUh6SKd+WUqd8fYxn5KtvIU5vUdYz\nT0/JaujwOcrZZf37f+Gbe9CTaY54K2Druq7rekx9ySPh4SNfUnNK4votY9cQpzFRfn2D4CPtX/mN\n7+Ypyrm/yX9rtsgx2B7l8a2B61RRPYR4ytfXUYXbN6yZMYRPmWQlD7ZhgeqDZPhmfory029foiyy\nUc66z2MDLRnlCoV/qC3EVtqceVQXsE/5entJTN70C/4l9K/hQXnU6kPs78ymsYMyfDO7R9llZ7Ra\nmNtvVBewq+OowukK8O2vZD/bRvVqN/vSbibKfUMcwR7c2TvwGs1BB7IhLWecMRorYwUs3eX0/kr3\ndHefEPsGzt7SpygTYgFqtvFp5mH0uiUixj70Nokjp/qsIsz5uAeNUUYETAj0IlIBqyae67SVPYjs\nVUKjHL1JHNgbQqz09Rrv6f09+mvDBqljjIbX3Lp3/BO6KRGIB4UoT8UZ2AyCN6bkcXRFEOqkpPze\nc2e7xc2GiKKAXYwZHxf5iNec7lEpfRxIA6I8+EJLhJPsT26Tn6tfvoQZQgxL8HD1UbrMeR0B7tFB\n8THsbKUnquqDTOEMzLLBOb2P3lDJHKOhnGLwMeVlfJUzMLNoaCKYEmXxJ7gDfVDA8B0nYYCg8d2U\n/RONH8gMIaKW8dTHC07xoR4FzJq9bPDICcc0fD2ml8skcIwCVmv6gd3YQW5cf6YjjzExzYhvbzw7\nphnXwKwZ+a2JmNU67NMo5gaLVAlieqbVooDVshlvWNKSY+TnJ1Na8PbJUOG7sxQw4JnxwxvwjQIW\nWOzRc+DHxG4KPaQ2TOIIjzLmG/GFX5yBBUbvD/AhajeFAga4RjcFfjGECAAwiQJmBI9FAPAk6oMF\nKGCmhMzRQOigAF9QwIyY3sOibQWgDAUMUGNiN4Ux6mHYyXIoYCgw/fwPvRHfCNx1U5hGDyClRA0b\nZfp+nr4CcihgABCDo9K1cTiEuDg6QcYtQhwBUcYrV2dg0TPe7Jf6lIse4hiIMgq5OgNb13U12IIL\nHK1hDnijIU5yIYoQartRTjECpIerM7Bb197c6di47e6NfE9KQp/195/jt2tioxMnymLZQpQ7vEcq\nOsvy512hovzJYmVFT06ROG7Fcgj+sljdQKTnKJ/CSpTtKjyQr/8Ekt0hxPWv8SvQdZi+38JtLdlx\nlMmfzfQQJ1O7a8zCDV2DjDiEqH/JXRduKDvLWdxdRvNnIouxIMpdcVYOADDJ6hAiACA4ChgAwCQK\nGADAJAoYAMAk/7MQBR3n/GyTX/ZXGufC7Pe4XBfY+BGnu+L219d1bVny6W/FV3siovy0So6jLB7i\nVHCMVC9cMEVfY6o8xBSwb8TvmD7m4nWBLR+Rv+u+feWP+S242hoQ5eNynlbSTZTFQ3xd5tPH1S08\niaZoJsqNqz0AQ4jfLMsie+9Fv/s3b5cstf4KU1kQUU7eQ5yko9z1RuzrwjmQN5yBfXPsrcxelxrX\nIYgWe69NZGl6EOWd3Z3wSvZYGEw2Re1GmQL2gdEY7wTX/zS84In1jZJaf8chTsa3iwN5xxBiKaM9\ntZ34+ttN+gyifOQyxMl4lDmQj6yeOU7RbxbW4PlpUjOX9j93OT8tcpSvIe6x2hN1mms3eBZiYyBK\noqw8xBQwAIBJDCECAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooAB\nAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBM\nooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKA\nAQBMooABAEz6P3qqAUEgyLLRAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 23, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo1()"]}, {"collapsed": false, "outputs": [], "prompt_number": 24, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Display the signal to recover and the signal recovered with $\\la=0^+$.\n", "Note how the noise impact the recovered vector. This is because choosing\n", "$\\la=0^+$ does not produce denoising."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAANo0lEQVR4nO3dYXOr\nNgKGUWsn//8vaz/QsqwhTuwA0iudM51OmjixAoIHYd/eUmt9AECa/7QeAAB8QsAAiCRgAEQSMAAi\nCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAA\nRBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiDRiwUkrrIQBwua/WAziTdAHMY6gV\nWK211tp6FB0p5aHp0LmLjtMZDv+hVmB706/J6sNGgN5ddJz+9cf2vx4YPGCP1vuglNJ6AI9aa/Nh\n/DsYw+hrGD2MwTCWxKxPfeIwlsP/0+8NuPAd6hYiQJwO2p1KwACINGDAergXAcDVurj1fJ1O7q03\nVIobFNC7i47Tv/zYiJPngCswAGYgYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIw\nACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgk\nYABEEjAAIgkYAJG+Wg/gbaWU5YNa63dfOvwqACMJC1gpZS3T9uOVbgFMYrRbiKWU7ToMgFGFrcB+\ntKzAnhZq+wcA8CTu6n+ogB3GSbEAfmN7toyI2Ti3ECM2NwBnCVuB1Vr370JcbhgefgmAUYUF7HEU\np/UzugUwj3FuIQIwFQEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRg\nAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJ\nwACIJGAARPpqPYAzlVKWD2qtbUcCwNXGCVgpZe3W9mMAhuQWIgCRxlmBfWe9rzit6TcABLjoOB37\n8B8/YO4lMolSHq8n+4sH/Pi9f/muz374RZ4G09XYuhJx6e8WIgCRxlmB1Vq9CxFgHuME7KFbADNx\nCxGASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAA\nIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiPTV\negBvK6UsH9Rav/vS4VcBGElYwEopa5m2H690C2ASo91CLKVs12EAjCpsBfajZQX2tFDbPwCAJ3FX\n/10H7N32HD5AsQB+Y3u2jIhZ1wF7qz2HL4kBMKquA7ZXa92/C3FJ1+GXABhVWMAeR3FaP6NbAPMY\n7V2IAExCwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgC\nBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBKwa5VSWg/h\n8TCM/2cYXY3hYRj/r5NhRBAwACINGDDXLwAz+Go9gDNJF8A8Sq219RhOVsr/filJYxrLnH8x4V88\n4Mfv/ewZtw/74OdfZDvsz37xWfRfh8EDBpNYLtVeTPwXD/jxez97xu3DPvj5F9kO+7NfnH6knuuf\nllbb30LAAGaQ+hqYRAFMbsB3IQIwA3fbAIiUeguxZ+vrc/s3Q95/ubC+IthqDE/P23ZTtB3Gi31x\n53h6G8b+P28bxosnbbspbh7GjzOh4WH7moBdYt3rtdang+TmM+b+ee8cw/5QbLgpGg5j+56j/b64\nbe/s/1RJk4m6H8br7XPDGFpN1MNN0WSivpgJdw7jXV4DO18nO7ifqVZK6WcwTdRae/j1n4bRakhP\nw2gyNw73yP0TtZ+J0XoIH7ICu8rkp+xVP3+ofL2WbD6SrpioCxM1cSYI2PmeFt0Nx7D998z8z1kO\ndThR406g52o1UXuYCZ8RsEs0nwr7V1najoc+9TZR2w5mZqEbX8BO9rToWe5xN38PT6sx7J935mF8\nN5gme6fPifpotDU6mSFNhvGbmdDDxDjkqgeASN6FCEAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEE\nDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAi\nCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAA\nRBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAABlTKo5TWg7jYV+sBXKsMvwMBjtXH386BtdbzBnOJ\nwQP2aL0PSik9TALDMIxux2AY1w2jlM9PgBFX/24hAhBJwACIJGAAROrinu+5tneQO7mpDXC/Uh4f\nn/8iTp5DvYkj4lVHAE4x1C3EWmv/lwwAnGKoFdihp2WZwgEciruJNX7AFAvgN7Zny4iYDXULEYB5\nCBgAkQYMmHuGADMYMGAAzEDAAIgkYABEEjCAlmb4mycvImAjc2AAAxMwgJa8b/pjAjYyBwYwMAED\nIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgC\nBkAkAQMgkoABEEnAAIgkYABEEjAAIn21HsDbSinLB7XW7750+FUARhIWsFLKWqbtxyvdApjEaLcQ\nSynbdRgAowpbgf1oWYE9LdT2DwDgSdzV/1ABO4yTYgH8xvZsGRGzcW4hRmxuAM4StgKrte7fhbjc\nMDz8EgCjCgvY4yhO62d0C2Ae49xCBGAqAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIG\nQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIn21HgBM\noZR/Pqh12GeEm1mBwZh0i+EJ2HvWq1p4y5KTq6NSyvMU/eMz7n/gnGyEPgnY20xlmMpyyDvwOyRg\nMI7TV3jN70P2sARsvhH4zlBv4ij/zvR65Yz74Gcv4/r7oKJflj9rIzCPJukq5XiW9jN1HUqrcQJW\nSlm7tf1484DHw15v4S+nIXttZrX21TB6c3CiD3UYsNL87gNApv7rMM4K7Duv98FnF/h/XxaccjPw\nrGE0maUfXORuN9r22/c/av297rmUvvRZlh/+mz31epscPuCzDbU+eP0h333vjw9YH/a4Zh4ezo3t\nZ67bdx8f46dsjb//XhFX/+MH7ArdX5f8Vttf5N1jbD3hJhxZJ3trQ61b6cfv6mQmdzKMc312/TTh\n3P6L2QPW6sg5ZXGQe9h/fJSuF/5PP+ppY+ZumThOuC98MA9N3beME7Ba6z3vQjxLwhiv8scX559a\n5SX3J/dvDbuAJsYJ2COkWyxO3Fd2+/3WYr24gLA4a2WeLe8PMjOsAQ7jiNf8fvMGDbiCgDE4J9DT\n1fqrVe89//vHF6bd9b/cQQMQMIbV/AT6d8uZ6PRf4baFXfTGp38CxsicQA9NslnmWYhMa6g3ccD9\nQu9TObMzACswACJZgcGfWMpAK1ZgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAA\nIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASF+tB/C2UsryQa31uy8dfhWA\nkYQFrJSylmn78Uq3ACYx2i3EUsp2HQbAqMJWYD9aVmBPC7X9AwB4Enf133XA3m3P4QMUC+A3tmfL\niJh1HbC32nP4khgAo+o6YHu11v27EJd0HX4JgFGFBexxFKf1M7oFMI/R3oUIwCQEDIBIAgZAJAED\nIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgMHsEv7iJzggYABEyvvrVIBz+WuICGUF\nBkAkAQPm4jW/YQgYMAvpGozXwIBPJMbAq32DsQKDKST2Bl4TMAAiuYUIUzj97tkvf+Cy8ivF7TvO\nZwUGTMpt1XQCdq3SxyFiGFuGcecYav3nn7bD+A3DiOMWIjApdzXTWYEBEEnAAIhU6tCraHeTmcZy\nIJvwnKb/OgweMJjEcqnmaGYq3sQBI5AuJuQ1MAAiCRgAkQQMgEheAzvf+tbH9Q0y+8/cOZjlSVuN\n4el5226KtsN4sS/uHE9vw9j/523DePGkbTfFzcP4cSY0PGxfE7BLrHu91vp0kNx8xtw/751j2B+K\nDTdFw2Fs/zjHfl/ctnf2f6qkyUTdD+P19rlhDK0m6uGmaDJRX8yEO4fxLrcQz9fJDu5nqpVS+hlM\nE7XWHn79p2G0GtLTMJrMjcM9cv9E7WditB7Ch6zArjL5KXu1vynRynot2XwkXTFRFyZq4kwQsPM9\nLbobjmH775kdvtRBhxM17gR6rlYTtYeZ8BkBu0TzqbB/laXteOhTbxO17WBmFrrxBexkT4ue5R53\n8/fwtBrD/nlnHsZ3g2myd/qcqI9GW6OTGdJkGL+ZCT1MjEOuegCI5F2IAEQSMAAiCRgAkQQMgEgC\nBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACR\nBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAA\nIv0XFwHpn5Cj9UYAAAAASUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 25, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "subplot(2,1,1);\n", "x1 = x0; x1(x1==0) = NaN; stem(x1, 'b.'); axis tight;\n", "subplot(2,1,2);\n", "x1 = x; x1(x1==0) = NaN; stem(x1, 'b.'); axis tight;"]}, {"source": ["Display the evolution of each entry of $x_\\la$ (different colors) with\n", "$\\la$. We display the evolution of the correct components (those indexed\n", "by the support $I$ of $x_0$) in color, and the wrong components (indexed by\n", "$J=I^c$) in black."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3dXXqj\nOhKA4eI8s6/YK+vuldlZGXOBjQUILIQEVaXvvZhxfo4hjsPXAgFd3/cCAIA1/129AgAA5CBgAACT\nCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACT\nCBgAwCQCBgAwiYABAEwiYAAAkwgYAMAkAgYAMImAAQBMImAAAJMIGADAJAIGADCJgAEATCJgAACT\nCBgAwCQCBgAwiYABAEwiYAAAk7wHrOuk665eCQBAed4DNqBhAOBOGwEDALjjPGCfkReDMADwxXnA\nAABeNRCwvn89YBAGAI40ELAQDQMAL9oI2DgIEybWA4ATbQRMpg0ThmIAYF4zARORvmcoBgButBSw\nAUMxAHChvYAJQzEA8OB/V6/AdYaGjekKGzYbpQEA9Gk4YIO+jwy/Zp+hZwCgT/MBk4QznWv0bLms\nvU+bstuT9ALwi4AFEud3VDpgVuNpx+ekZADcIWDrPM1XHFaejAFwxGHAuq7ra2ypjW79w+4yIAPg\niKuAdaYHSZXMJlsOGJABsM9VwIaBFxmLiE5UYUAGwDJXAYua9azK3kVDGJAB8MJ/wFovVtRGxoSS\nAbDBf8Cwau0EOAZkACwgYOAIGQCTHF7Ml32G+WaXOR5wsWMAKjECwwIDMgAWOByBoRgGZAAUYwSG\nbxiQAVCp9YB1f+ODif4v2+UFziEDoEmdywaqsXFdxLV0LRGziLW9iK7fTgBUaTRgs3qNifpaNWI2\nR8kAXKTFgIWV2gjSdswo2RwlA3CutgK2NvD6/jzrMaNkc9GSuX6bAbhEQwHLrtf8OZn3kYiSAaip\nlYAl7jbc9+SxkpGxOTIGoI7mAlY8MAzIUi1L5vq9B6C2JgJWr16TZVGyFGQMQCGtn8hc0Npc/HPy\nacbybGhOhQaQxf8ITP68Hp/fDw6SfcFoDMAB3gM2JuRfLxdtHtm1+AUZA5ClrYANrvqJKdkWMgZg\npzYC9m/1Z2RMpgsZA5DMe8CmN/3YuI+VnpKRMTIGIEUTAZv9iNou2seALI4zoAFsch+wTmTrR5xt\nJK99MShZBBkDsMJ/wFJ+wHAjqeH1oGQR7FcEMEXAwm9+PdDzklCyOTIG4I2Ahd/8eaztVWG6x4Sq\nPb8ALkLAZt//eazwhWFANkHGgLYRsOV/8nms9rWhZB9kDGgVAYv+V5/Hyl8eSibCgTGgUfYC1r23\nVtE1H776uQtzVsDEVMMGlIyMAa0xFrAwSMs4fW5fuXiQtazPY0Mv0lrJpJGYsUcRaAYB217c57Gp\n10lks2Qjt0kjY0ADvAVseBB+z+wZXM7p2JZSsoGrnrFHEfDOT8CiXzo4Ans/2+exqVcrLrFnTmLG\nUAzwi4AlLvfz2NQLlsT/zkbfvz+gVcYCJiv7CcdcrX2pxHI/j629Zvuc37MzlshQDHDHXsB2KRgw\naalho/TjZ9dKzVuDv0LALwK29wk/j12/chFWYjaKV42hGOAFAct4ztcD16+cAdlB7f/NPuYXCZhE\nwPKe9vXA9Ytn0o5zBoKMdX+Mz1IBmkTAsp/59cD16+fE6nW2pkOx7k/wJXoGqEfAsp/589j1S+hc\n97ebDcWi6BmgEAE78uSfx65fxQZMZ3asZWxEzwANCNjB5389cP0qtiE2O7Gtq5YA1hCw44t4PXD9\nQjZjc5I9PQNUIWDHF/F57Pq1bEbyiWL0DLgWASuylNcD169lY/af7+z/kpKAMgSs1IJeD1y/nI05\nds2O7Z5RMuA4AlZwWa8Hrl/R9pS49FTrt8kG6iBgBZf1eez6RW1SoSsoMiwDCiJgZRf3euD6RW1V\n6asAMywDDiJgxZf4euD6dW1YhYvZUzIgDwErvsTPY9cvbcOq3ZCFHYzALgSsxkJfD1y/tM2rfF8x\nhmXAVwSs0nJfD1y/us3rFo2pce8eSgasIGD1Fv164PoFxnm3eKZkwAwBq7foz2PXrzFE5LyMvZa2\ndoczSoaWELCqS389cP0a4+3chr2WybAMDSNgtVfg9cD1y4zAFRkTSoYmEbDaK/B57PqVRuCihr0W\nTsnQDAJ2wjq8Hly9IjjXpRkTSoYGXL99r0pDwISGtezqjAklg18qtu/1aAuY0LAGKWjYgJLBGRXb\n93qUBEwYhEFNxoSSwQst2/d03XtDsFzz5Zf0BExoGERXxoTzyWCcou17ijBIszgN9Ro+M35JVcCE\nhkFOugDVXtGSkTEop2v7/tXXgA3C75k9A7PqoYKyodhoWTIyBrVcBWz5JW0jMGEQhhAZAw5Qt33f\n5iBgQsMwQ8aALBq37xucBUwUbaxwKZUHxgZkDGpp3L5v25hqqHwWYohBGCK0DsWEjEElpdv3UtQG\nTGgY1pAxII3e7XsRmgMmNAxrFO9RFDIGNVRv34+zEjDRtYGCDmQM2KR6+36c8oAJgzB8pXiPopAx\nXEr79v0g/QETGoYUZAxYMLB9P8JWwETddgma6N6jKGQMpzOwfT/CRMCEQRjS6R6KCRnDiWxs37NZ\nCZjQMOxiLWM0DDWY2b7nMRQwoWHYi4yhbZa27xmMBkw0bougEgfG0DBL2/cMtgImDMKQh4yhSca2\n73uZC5jQMGQjY2iMve37LqYDJuq2P7CAA2Nohr3t+y4WAyYMwnAcGUMDTG7f0xkNmNAwHKd+j6KQ\nMRxjdfueyG7AhIahCPUZ48AYshnevqfwETBRt82BNWQMHhnevqcwHTBhEIayyBh8sb19/8p6wISG\noTjmd8AL89v3bZ4CJho3NTCLjMG+/65eAXyhb8MCF/p+8t7qusgOxkvNirXcuwiYH6BsczACG7Aj\nEbVYOzDGUAwjJ9v3NW4CJjQMVZExGMQuRHuU7emBC7M9iqJupyJ7FLHkZ4AS5WkEJgzCcA7dozGG\nYhgxArNk3Ixo+pcx3ImOxtRgKIaRqwHKkrMRmDCrHidTPNueoRjsbd+791/U2pqH0fIXMGFHIk7G\nHkVoZWwX4hCkQRfbrRH9pDPsSMSpdM/vYI9iy4wNULZHV8Nn3I/ABozDcAHFozGGYg0yNgLL0E1d\nvTrlefyZoJXi+R0MxRpkbICyMbqaxWn4kuMRmDAIw7W0zu9gKNYOY9v3lN2DjexCHNAwXIyM4TrG\ndiEOh7gGYaWuXSsNeA1wDfYo4jruByjOf0BhEAYlGIrhdMZGYFhiVj1U0Hp/FoZijjkfoLQwAhsw\nDoMWDMVwFufb9wYDJoq2GGgYGUN9zrfv7QRMThmEFZkv085vpHU0DJU53743FTCp07BKkzyb+r00\njYyhGufb92YDJgc2FInF2vvCrj1tU7+gRhlpmJAxa5xv31sLmGQ1rFKx9i6utd9Uc4xkjIYZ4nz7\n3mDAZLKhOLT3r/ZLR8laRMZQjvPte5sBk9dWYne9rnqtKFlbaBgKcb59bzNgG7sElb8alKwhZAyH\nOd++NxiwaQM+P7utlyFastZ+lf4ZaZiQMa2cb9+bCtjyhjJatw+pGJA1QevblKGYfs637+0ELHo7\ntPeXZPqlc9aopO15ko38ip3T+jYlY5o53763ELCNdAXfI9PvqbpGFX2d8e/+1+2Z1rcpDVPL+fbd\nfcBS6hV88+RD0y9Myrlrvn/1bml9m5Ixhfxv373+gLvSFfxXkw/dvDbsY3RF69uUyR3auN2+D7wG\nLK9ewX8++dDZK8TUDye0vk0Ziunhc/s+chmwcAOd/dNp3TgUtoyZv/eDc1pvFETGNHC4fQ85C9jB\ngVfsCScfOnqpJsiYbVrfpjTscq6270ueAlZk4LXyzJ/HXl6tCDJmm9a3KRm7kJ/te5SPgBUfeMUW\nMfnQ/mu2iowZpvVtyuSOq3jYvm9wELAT6vVe0ORD4y/bF1yqyjCGYngzv33fZjpgp6VrutDJh2Zf\nvCRkzCrFb1MydibD2/cUdgN2Sb3ei558aPP124f9iiYxFGue1e17IqMBqzdfY886fB4bfAlzkDF7\nFP9ri4ydwOT2PZ25gF048FpSvHGoSNWvAEm0/muLyR21Gdu+S7B9Wa758ku2AqZh4LWkdeNQFxkz\nRvG/thiK1WNp+y7TIC3jNH5m+UA55ZtLxRuHitijaI/if22RsRr+u3oFSjK6fVFeLxHp+8nWIOFC\n8B70fb/8F1LKVfBxmdnbVNMva1as5d5FZLAxQBltj8Ci3zP7qrafV+duwzVtDsXEwhsJEwzF2uAq\nYMNWZnozYr0/oP6BV1SzDRMyZovidyqTO0rRu32PSjwGtvEZJWwNvJYU/wO3OqP/8miU4ncqQ7Hj\nlG7fN6xNNYxuVhQGzM3mT/GW4Qxufo/+KR6KCRk7Rt32vSxtAXO21aNhs89Y/4V6pvjNSsOy6dq+\nF6cqYNZ3G0bp/tftGciYGbrfrGQsg6Ltew1KAuZs4LWk+F+3JyFjZih+s9KwvVRs3+vRELCr6rV2\nxlKlFVC8WTgPGbOBoZgX12/fq7o8YOfvNsw707bIutGwgfsBtxO6369kLAUBq7jo8MPaq1H2ChHZ\na6t7m3AqMmYAQzHjCFit5YYfVttrtxqtxCWmZG/vytOwETfMtEH3W5aMbSBgVRY6Pq6x9OPdynjm\nXc+v+9+1ZyNjBuh+y9KwNQSs8OLCD8suul639i4ucaG6/117AUqmne63LBlbImAllzU+LrXQk6OV\nsRrba6J7g3ANMqaa7qGYkLEpAlZsQePj40tU0q01uybo07AoMqaa7nctDRsRsAKLCD88sjjl3VpK\nLNnyu1T+NNegZEoxFLOAgB19/vFx3oLMRSsqpWRkbAMZU4qhmG7+AzY8KP5jZg+8Ss300+lrycjY\nNkqmDkMxxQjYoadNfPIap1ulL1HDmXChvu/J2DYypo7uoZi0mrFWAjYoPr1icwLe2SMtnXsjN1+H\nyFq5fj/uQ8Z0YSimT1sBk6JHqpZ3f/76PKedGXbOCuyyq2Su35W7UTJFGIpp4j9gy5s1y+FDVuHN\noL8+Q8FXOLFY0R95+T0l1ijT+upN1sr1e3M3MqYFQzE1/Ads7UsZB66G/+TkNpS9Su+S1pLJGDPX\n79AclEwFhmIKtBuwwfLHv/B2JNlLV6LSOXBDyVy/T3OQseupH4qJ94y1HrBRym635X9SaulFHJlR\nUknGu+vbqvp+w+5Gxq7HUOw6BGyHjMgVXHSR56lz/5Qdr8muJ8++BmODuBn0lRiKXYSAKaLn+rxS\ncX9g5iJ2XYCxWQzIrsRQ7HQE7EqXv/inTUgpuKCuExE201vI2GUYip2LgKVKeaHOPxWslDNPgi4S\ns7WMJf7njWC/4jXUD8XES8YIWHWGXuFLruVxpGfv/5SSbWFAdgGGYqcgYCc54ZDSQdMrx1+ThLyY\nbWds7b9qDRm7AEOxylwFbHnpXj0BM2djyuX59wj9uhrT76VkW9iveCqGYjX5CVh466/xMQGTb5un\n7JdIz5zJjJJJ81ttBmSnYihWBwEDAETor8N/V68AAAA5/nf1ClhV/DBbqasyejU7Jre8l03By1YV\n+YdnvWMf9Xb1bLyMF54yqH8cgKsQsEwnHDqqcUmnvu9vt9v44fP5PPJsX5clIs/n53lut9RIL/u0\nfY/pUleoKrcVDp+zyFOK1D9KsfcWd7VvvEC6sM3PMTBRMwtR4UV1q5r8tNOfPVqvvf3YTle6c85y\nqzTwOucA+4Vn4he/eTpa4CpgSy6Dca342yX2LprVK6MfRTZq515kJHzyQs/5dzEuKVqvlL+Rhzxu\n/a3gQjdWwPcWCWWxC9GzcVsw7Db8/f1N/P74HOvFd288VViv+10yrl54fLt25nlsNQZeldKVWKzx\n8WnpEuqFnRiBVXfyK9y9D2v172NdiQe9JKVbaT/LUK/7Pf7VvVeVLHUssObVsGYLOvyERdOV+FcQ\nRmtAuqAcI7Avjh9OO3MPSU691i+L+1rRc/uR/XJdd/mrcFmHn61EunbdzubZPZdfIl0wgYB9sZwY\nIt82ENv3vTyzZ6v12rx6YNlDTdkFygtexnKzlR14HUlX3g3YSBesYxdijsSYZYzeDs2ymw6/xnp9\n0vXtqrcF03Vwv9/X//zy+1sWHHhlpOvILUNP7taAmRqogRHYjqHV8tv6vt+o1PaJt9vPvPYfJprU\n6+vV2ouec/Z4RM73Sn8qW+mSY/VKTFepfwCRLjjTdMC2z4dN2Wqkb1lqH0JbHv1KWXB2LZYe8xkA\nSdK3bpd3670as6XnPs96una9VVJ+/Eu6JaQL9bW7C7HUtRvMyb5W0MazbV9uI2Vxe9fqkvdtkd2G\n83T92/Hf7v2plaRLqBfqaHEElvG3tDYy256voUrGntLtJxlVqpeq+30UGXh1f7uquRpd1S0hXThX\ncyOwurPYC10/8IQFFVkN2V8vc+mSAwOvvb+m4z8j6UJT2gqYnnswpkiN2XD0a+204UJLiSpeL20b\nwb0Dr+KHrxJd2C3R91tDOxrahXjyH1XZM6C/r/zj0d9u4+TD6FWj6r0CG/VKHE4p3Ah+rde+3+yf\nKlfgvTZdwkwNXKqVgF34d1Xk4NNazMLJh+Fpy72Uu5rsinH4NavX17O8Z58xlK4dv74/02cona7L\nuyWkCwq0EjAlysZsbZPx/HbR3uOmVzvMOalW24GuUbBe3eIzqz4D7pr3PYl2S0gXWtXEMTDlP2Op\nGRk/Pz8y1mvnBXMrMTFHY3D88FXV+54oHHINLv/FoWX+Azb8gF9HLUqUikp6OQrKOC36umnxJWdb\n1EuXziHXQPmfElrQRMAs7vG46vSyry/O2qGvNUrORK40P/DkdJ3cLdH0bw5gyX/Avn6P8ldAw4nS\n40v0tV56XvA9r9tnfdJXrVK61HZL1P+loEENBWz7qhkmXoeN9f/5+Xn+/u5t3Wz/am0aziIPTm+Y\nfT55QRXSpWRXoZAumNJKwFJmb0e/TaHu+dw4bXljouPJl39UMtJaORw4+560ZZ2VrvO7JewthEH+\nA7Z3ZoH+F2R27lf0nOXXNxT6WcY9h2vdVHVAa3tllKRLT7eEdMEszgOb70ZLad6Flmcuj7Pnlxv1\n4kek1E4aTFkxDenSs6twQLpgGgF7CY+QKW/Y0vP3N9zw7Nr793UTtnfa4XE1bmBdpF7F03VVt4R0\nwQUC9qG/YdELRy2vu5FyzYutpSy+P+9mlUeWuCbjl3JtuvQPuYR0wSwCNqG/YXGb63n8Fi3RS0Zl\nvzhVczVd0OzZ0v6rEvsMtXVLGHLBIzvb6Cx5Edp3GfizbA2/6qzk+aeglXq1lyue8sTH06WwW0K6\n4Je9EVjKzPiDf5/GxmF1Vu/57IY9h+Ghr+JJK/7aXpIund0S0gXvjAUszEk0LeOUQmcNSzz6VUp4\np8rQtxnqV16G4/x0qe2WkC60wVjAtpX9E9XWsLgKa7X3PsvBulw1yX7+marp0twtIV1oiauADWax\nOXKespKGrQ6/Kq/PafPms52ZLuXdEtKF9qgO2N72RO/+5WlfYnjP5UrOP+srz2np0t8tIV1oleqA\nZfwR1vi7vbZhk+HX+5OVhl8m6pWXLtlzatdatERZt4R0oW2qA7YUtmQ2m2P4fKUbV17VsLFeEgy/\nKs3dWJu4oUrVs5INdUtIF2AuYBL7Kx0+U/uv99pxWD/be1hz6TqHX/XOSrbeLSFdaJW9gF3o5Iad\nOXVe887DSumy1S0hXcACAdtHxZyO0gt1Vq+NdG1ES7R2S9hbCKzQenpTIZUac8K1pk6bOp991ldt\npdJlNFoD0gVsYASWQ8U4rDQ99cqYZ7hM1+PfQ0Se/57R71feLSFdQAIClqlqw84ffqmt1650DdFa\noz9aA9IFJCJg+cJbOUfPoT6u6twNbfXala6xWz6iNSBdwC4E7KjiQ7HomcvDkg4+c0hVvXbtMxyO\naT1ktVu2ojUgXUAGAlZApd2J9YZfqs5Z3hh4bc+/mLHYLSFdwAEErIxSDTth+HXttMOwSXe5hV96\nyFNE0ttqtFgj0gUcRMCKKTsOqzR34/x6RUdRs3TJu14b7n/uw4O9d0lWiHQBRRCwkg42bHX4Vchp\n9dre9RcdeM2GU8tp8WI/XVxKAyiLgBVWZBxWY/h1Qr2+Xpxp5XDXbfyMy24J6QLqIGDl5TVsHH79\n/P1bZ71eytYr/YqCG5M1ot0S++miW0BVBKyKIuOwGsOvgvVKv9Pj3nS57JaQLqA0AlbLroYth18F\np86XrdeuOxSvpcvfkGstWgPSBdRAwCpKvFRH1anzpU752tWtwbJeTXWLaAG1EbDq0odixc9cLjJx\nI5qu7XOwounq/s6/zWi3GGwBShCwM2w0rN7w62C9MoZcsnJRqJTbIuvHYAvQhoCd5Os4rN5V5/fW\nK2PIJevpCgde5rrFYAvQjICdZ9mwemcuZ0zcyOvW4Os+Q1vpYrAFmEDATjVrmDxel1QvO/zaNXEj\nb1fhKDrPcLbP0ES9GGwB5ji5m/AanbdLnmwrH4/51PljK5x46OtgtyTtcJfydBEtwDRGYBcIx2Fy\nv8vPT/i1UkuJ1ut4t8T44a7taAndAuzQOEApSOcITIbJh/f7+OHPz0/Z4deyXvXSJer3GRItwCWl\n2/dSVAdsEGTs4KpGdx4W6ZYYTBfRAtxTun0vRWfAZheO+g1OW85e22W9jswqDCWmSxTUi2NaQFPs\nHQMbN1Jr2yOd0Rotp87//PyMDStxJ8xIvfJuXryWLtE08CJaQLOMBSzcvke39V93HOkxmzp/5Or1\nG4e+XKaLaAEQcwHb9jo7WHHDNq46n30HluVZX+Pwq+AOQ7l0nyEHtAAsuQpY1Gzbp3FL916ljIZt\nHPoqNU1DLkpXyj9ENP42AZxFdcB2tSe8a0m49dezjUu56Vf2OGxj4kbSuiXfc7JeuigWgF1UB2zX\n1mr72Jhei1VNb9js0FdYr+wLapycLqIFIJvqgC2FG3dbxYoPv1ZWO6VhGxc8zLuMoZy1z5ADWgCK\nMBYwiW3dZp/RvPkb6/XVdsOOH/o6eeBFtAAUZy9gFn2uu/GWctX5lHFYRr1OSxeT3QFURcDOszF3\nY020YRuHvr4Km1IpXdxMC8A5CFh1y+GXyI6L9s4a9r6D2Ev6xI21gVepw110C8DJCFhdu+ZurAkb\ndr+/7oI5u+hGkXplpGutW0QLQG0EzIZZw177EhMOfZ2cLroF4DQErKIiw6/R4/G59UrXdQ957Uzc\nW6/j+wyX6aJbAM5HwKrLmLuxNEzcCBt2l/vYsKjlfI3i6aJbAC5EwGo5OHcjFJ71NdmXKPf1e8rM\nl1l8nyH1AnAtAlZXkeHXaJi48ZDHXV4DsZV7ynwe9z3pAuDTf1evgE81hl/hWV/hzsNZXSYf/ekO\n1qvruvD5+76nXgCUYARW0fG5G7MLHobTDnuJnOM8q1f43zJTA4AzBKy8YfiVftnDFLdbv7zixuwc\nZ5EgMEG92GcIwCUCVt3x4desXuG8+ekdqN8NO1Av0gXACgJW2Gz4lT13Y3boa7Q862vesD/Bl/bU\ni3QBsIVJHCWVmruxcegrstBOJjsP/70Xm1yv2UwNoV4ALGAEVl7BqfPhzsP1eg16kfcH/3YUiHQB\nMIqAFVN8+BWduDFZ4nLO4b/xS1/uUs0kQwDWEbDCDk6dT5m4MYjOmN++ByZ3PAHgCQEro8jcjdmh\nr1Fqvf72EmsY3QLgEgGr5kAbtg99bdTrveTZ+WHLVaNbAMwjYAWUHX4drNfrM4uBF9EC4AwBq2Nn\nLcKdhxsTNxLr9V4FigXAM84DO+r48Gty6Ov+uUrvbPi1q14A4B4BOyR+2cPcoU940Q3qBQDbCFhJ\nGVPnUw59US8AWCJg+cruPKReALALkziK2jP8Wjv0FfrU69jtUQDAH3sBG2eHL2fZbXyp/GoUuuq8\nyOrEDeoFABuM7UIcLi0xmJ3nNHwY/dJJsoZfaxM3qBcAbLM3Ats2Zmz2mVGRwdnB4dfXQ1/UCwC+\nchWwMU7hdWyL706M1yt5KdFDX9QLAPZSHTAzF0PKW7HYxA3qBQCJVAdMYbGK7TyMTdygXgCQTnXA\nlsIJGrMdhtEvnbNOid8YvVvKRr1IFwBsMBYwicWp3uGumWJT5xeHvqgXAOxlbBr9hYrN3aBeAFCC\nvRHYJYZ6zR2o1+tpuUwUAORiBLZD3s7DtUNf1AsAjiBg3xW7Z0qw85B6AcBBBOyL5c7DXfdMiR76\nol4AcBwBS3J05yH1AoDSCNiWIzsPl4e+qBcAFETAVh3cefjxHn5RLwAoiIB9UWrn4V1uny9TLwA4\njIDFldl5SL0AoBoCFnFk5+Hy0Bf1AoAaCNiqo9c8vD+EegFANQRsruDOQ+oFAPVwLcSJMjsPV+pF\nugCgIEZgERk7D2eHvqgXANRGwD7KXPPw/qBeAHACAvZSauch9QKAc3AMTCSo18Gdh596MWUDACoj\nYJF6vSQMv8Z63e/BN1MvAKjP+y7ExyN+M+W36Ff3XvOQegHA+bwHTERWKtU9n+HnJzsP9xz6ol4A\ncInmdiFGYzbfeZiAegHAtbyPwO734f+Hbi3r1d9u/e02frhr5yH1AoALNTcCG4TRur0f79p5SL0A\n4Fpdn3GHRju6ruv7fjbwWq2XpM48pF4AcDnvuxA33YKSiVAvALDE3gis697NiK358NXxS8MILPo8\nt6xDX114yUMutAEA1zF2DCwM0jJO42c2ujWgXgBgnbGAHTfbbUi9AMAoY7sQv47Ahgfj58dc/f7+\nisjPz0/4/dQLAOxSPQLrJt2IH/QKv3mjbZF0SeqtUqgXACikehZiP5XxDL+/v8/ofVIkp14PeQr1\nAgAdjO1ClNh+wnDuxtqXxv9Y3t+RtqzP48ejuz+pFwBooXoXYtSyuONnUmOc8G2zdIkI9QIAVeyN\nwHZ5ds/hwa2/pQ+/qBcA6Kf6GFhJ0/kgid841Es47gUA+jgP2E1eV6N/yuP1qS9TGT+Px3rdbtQL\nANRxHrA56gUAXjgPWDcc/RKRcBAW/U7qBQCmOA/Y4NOw9xeWLwUAAALkSURBVJyOUNfF6wUA0Mx5\nwB4PeT6757Nba9hsbsdwqvKA4RcAaGbvPLDjnt3z1t/i6Xq8Jn1QLwBQzvkI7H7vh/tPhoOwu3yp\nFwBAP+cnMn89++uzzzCoF8MvANCvxV2Ig/BwV4h6AYAJ7kdgncjkB5z9uN3fTkQet9eH1AsArPAf\nsOEHfD4/OxPHSnV/uzFdQr0AwJRWdiHebn3YsNnACwBgTisjsEHYsBBjLwAwx/k0+ploqKgXAFjU\nVsC6v939OfmMknp1yXd7uRArWZCJ9WQlCzKxniZWMtTQLsThuNdA2/29Zrs6dWIlCzKxnqxkQSbW\n08RKhloZgWmuFwAgg/9ZiGG6hHoBgBfGBox7zeol/y5aDwCwRn8d/I/ABq+B19+LVwMAUIrzERgA\nwKtWJnEAAJwhYAAAkwgYAMCkViZxqDWeOTieA6/qqORyrRSup4mVHFn5jYtI3/fKV1Lzb3x2VQu1\nL6a8V1Xz23INAbtM+P6eXDFE2cnw49t6+AvUuZ4mVlKCX7rmlYyumJ6VnG1t1f7Gw9Xouk7nSsri\nX1Q6V3INuxAv0/e9/veH/jUUIyspRrYI8t7aXr0WX8ySoJmV9bSIERi+0/8XqH+Da8hsV5JCtkYJ\n+o1jL82/9CgChi2zvQpq6d/mDusW/q9O+n/XtihP7GzH5rUrk4FdiPhC85+f2Pmr699E8Utq5cUE\nBqr/ddACzXPSltOoRPd6ap6TNtL8G5fFWplYyehnNJgNv9Su5PBA8298DQEDAJjELkQAgEkEDABg\nEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABg\nEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGASAQMAmETAAAAmETAAgEkEDABg\nEgEDAJhEwAAAJhEwAIBJBAwAYBIBAwCYRMAAACYRMACASQQMAGDS/wEU/IdVzVFDzwAAAABJRU5E\nrkJggg==\n", "output_type": "display_data"}], "prompt_number": 26, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "I = find(x0~=0); J = setdiff(1:N,I);\n", "clf; hold on;\n", "h = plot(Lambda, X(I,:)', '-'); set(h, 'LineWidth', 2);\n", "h = plot(Lambda, X(J,:)', 'k-'); set(h, 'LineWidth', 2);\n", "axis tight;"]}, {"source": ["Note how the noise enters into the signal (wrong selected indices) for\n", "small $\\la$. Selecting a $\\la$ proportional to the noise level $\\la$\n", "ensures that $x_\\la$ has the correct support and is close to $x_0$.\n", "\n", "\n", "\n", "Display the evolution of the sparsity $\\norm{x_\\la}_0$ with $\\la$.\n", "Note that this sparsity is not necessarily increasing, because when event\n", "3 happens, one element leaves the support."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAPK0lEQVR4nO3dbW/i\nuhqGUXtr/v9f9v6QIZOGl9ISiG97LW2NaE+lE0HLhZ2HUFtrBQDS/Hf2AQDAbwgYAJEEDIBIAgZA\nJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBp8IDVWmo9+yAAeIM/\nZx/AJ1QRA/ih1trZh/CNKQLW/8NQa3WQh4g4yBJynA7yQBHHuT3IiNf9g28hllK6/50B4DfGDxgA\nQwreQtyucFtr65f9r9MBeF3Atuw9u+3aO7dtIQL8WMRJu+wtxFrrt2caE85EAvBjwVuI5bJb+ETD\nvmw2vveYADJFTB5uBQfs+RSJFsC3tk+VETFL3UJ8/s4VL4Ahpa7ArscOTSECTCU1YOUqVLoFMJXU\nLUQAJjd+wFyQHmBI4wcMgCEJGACRhg9YLSbpAUY0fMAAGJOAARBJwACINHzAlot0nH0UABxt+IAB\nMCYBAyDS8AGrZugBhjR8wAAYk4ABEEnAAIg0S8BM0gMMZvyASRfAkMYPmClEgCGNHzAAhiRgAEQS\nMAAiCRgAkQQMgEgCBkAkAQMg0hQB81YwgPFMETAAxvPn7AN4Va21tbbcWL/ZrLkARpcdsPr1Qoe6\nBTCP4C3Ede21/U69f+1eV/UFGEn2Cmxn3Uvchq3WWkpbb1ulAdz0YAHQp9SALXf0+m9r7V6ZWmvr\ng6JeAPdcvfTvXeoWYrsofxP16L6WLYDxpK7AdrYNs8wCmEF8wNZc6RbAVFK3EAGYnIABEGmKgK0T\nHgljNQA8ZYqAOTsGMJ4pArZSMoBhzBUwAIYhYABEEjAAIgkYAJGmCJgxeoDxTBEww4cA45kiYOXS\nMCUDGMYsASultGYLEWAcEwUMgJEIGACRBAyASBMFbDkB5jQYwBgmChgAI5koYCbpAUYyUcAAGImA\nARBJwACIJGAARJooYMboAUYyUcAAGMlEATNGDzCSiQIGwEgEDIBI8QGrl6mMevHgh30kGMAw/px9\nAC/Z1qtdzm5tbwMwquCALaF6vOQqm8iVUkpppZRajXIA7H37dNqb4IA9absaS3t0AD7n67NlwNNl\n6jmw5c7d/vsMk/QAw0gNWLsoX181ADCJQbYQtyfD9AxgBvEBW3OlWwBTSd1CBGBycwXMBekBhjFX\nwAAYxlwBM0YPMIy5AgbAMAQMgEgCBkAkAQMg0owBM0YPMIAZAwbAAKYLmBl6gDFMFzAAxiBgAEQS\nMAAiCRgAkSYNmEl6gHTTBUy6AMYwXcBckB5gDNMFDIAxCBgAkWYMWGvOhAHEmzFgAAxgxoBZfgEM\nYMaAATCAGQNmhh5gADMGDIABCBgAkf6cfQC/Vy/DGK217ZfrdwAYWHDAyiZdyw3dAphH8Bbida5q\nrfXpGXnD9ADRsldgu1ztFmQPfgaAnecXAJ3IDtj27Ne9Mt38vqtJAew8eOnfp9QtxN2dG3FfA3Cg\n1BVYa+3eFKJNQoAZpAasXIVKtwCmkrqFCMDkpg6YE2cAuSYNmHQBpJs0YMv5MmfNAHJNGjAA0s0b\nMO9lBog2b8AAiCZgAESaN2D2DwGizRuwhYwBhJo3YCbpAaLNG7BiEBEg2dQBAyDX7AGzCAMINXvA\nAAglYH8XYdZhAFkETLoAIgmYeXqASAJWilEOgEACBkAkAQMgkoABEEnA/nEaDCCIgH2hYQApBOyv\ndZhewwAiCNg/24bJGEDn/px9AD1aSrY0zBucAfoUHLB6WSW11q6/fJ2MAfQsOGDla7rWbtVaj2pY\nkTGAXgUH7MBKPfH/VcrXGUUxAzhXcMDKZtvw+Z95JXvXGQMYxjPPqF3JDtju7NeDnwHgse2zZUTM\nUsfoT7xzW/N2MYDzpa7AWmvvnkL87gBKrc6EAZwmNWDlKlS2CgGmkrqF2AMbiQAnErCXuO4UwFkE\n7BgaBvBhAvaq3SWAlQzgM4KHOPqxNmyxNsxYCcD7CNjxrktWxAzgaAL2RttoWZYBHEvAPsQGI8Cx\nBOzTbpbMx7UA/JSAnebmqTIAnmSM/nxKBvALVmBd2DbMRiLAM6zA+uLaVABPErDu+LwxgGcIWKc0\nDOAx58D6tW3Y7sSYsXsAK7CuSRTAPQIWzLgHMDNbiL27Phm2G7W3nQjMyQoszPUbxUwtAnMSsCQP\n3ubsDWTAbGwhxvj2Ih3L/2pHEZiEgPVuCdLzl5jaZqwoGTAuW4gZftohFwgGhidgw7L2AsZmC7F3\nL3ZIxoBRWYGNzHg9MLDgFVi9PDe31rZfrt8BYGDBASubdC03dAtgHsFbiNe5qrVWW2YAc8hegZXN\n8qtcLcjWH9j+/FSrtJ++hwyYWdwCIDhgy329q9e1qYq1Y4gDeN6Dl/59Ct5CLJu7O+K+PsvEBQdG\nlroCW4q1HUTcDSUCMLbUgF1XSrduWnYR3TfAeLK3EAGYloBNzalDIJeADW6dpH/8AwBxBAyASKlD\nHDxp+8FgRjmAkViBzUvSgGgCNgtX5QAGI2AT2TbM8gtI5xzYXKzDgGFYgSFpQCQBm46dQ2AMAjaj\n1v5lTM+AUAIGQCQBw2QHEEnAAIgkYPz1zCLMQg3oh4Dxz+M+fXthe4BPEjD+MZEIBBEwSvkuXbX+\nvfSUwgH9cCkpHlk2DLfdWkYWlQw4nYBx23W6ALpiC5EbHm8Yet8Y0AMrML6w8AJSCBilbEbkn0+X\n1AHnsoXIP2oEBBEwSik/HpFffljwgBMJGACRBAyASMFDHPUyyt1au/6Sd/OOZuBcwQErX9O1dqvW\nqmEAwwsO2JOVql/fc6ttR/nF5D3Qs5p2hYLggC2W9daD+12x3sq9C8PYPltGxCw4YLudQz7MBaWA\nc2VPIarXuTQMOFHqCmxZfm0nD00hAkwlNWDXldItgKlkbyFyOruIwFkEjJesw/QAHyZgAERKPQdG\nJ9Yzj97RDHyYFRjHcDIM+DAB4zAaBnySgHEkDQM+RsA4mIYBnyFgHG9pmIwBbyVgAEQSMN6iNXuJ\nwHsJGG+kYcD7CBgAkQSM97IIA95EwHg7DQPeQcD4BIP1wOEEDIBIAsbnuFw9cCABAyCSgAEQScAA\niCRgAEQSMD7HGD1wIAHjozQMOMqfsw+AWawz9EvDjNQDLxIwPm1J17oUUzLgdwSMc1iQAS+KD1it\ntbW23Fi/2Twd5tguyDxuwPOCA1av5gF0K5d9ReCnggN2vfBabstYrt2+YlEy4L7ggF1bk7Zt2G6h\nJm8RfIQYfN71tlbnxgnYvTIpVigNgw978NK/T4O8kTnivuanNAx4YJAVWGttbZglF8AM4gO25kq3\nhrQswjy2wLVBthABmI2AEcCZMOCagJFBw4AdAaN3rTkHBtwgYGQwUg/sCBgAkQSMGBZhwJaAkUTD\ngJWAkUfDgCJghNIwQMAIY6oeWAgYkZwMAwSMVBoGkxMwgmkYzEzAyKZhMC0BI97SMBmD2QgY8aQL\n5hT/icywTtX77GaYihUY43A+DKYiYAxFw2AeAsZozHTAJASMAS2Xm9IwGJshDoa1bZjhDhiPFRgA\nkQSMkW0n7IHB2EJkcN4lBqOyAmMWphNhMPEBq5cnpHpx7vHQM9OJMJLgLcRtq2qt7bI9tL0N15aG\n+R2BdMEBWyplycUvmLCHAQQH7Em7wlmcsbCXCDtx64HxA6ZY3LOOdfgdgfL12TIiZvFDHPAKYx2Q\na5AVWGttfb1gycVPGeuARPEBW3OlWwBTsYUIpZjpgEDxKzA4itl6yGIFBjdYjUH/rMDgn+00x7Zh\nFmTQIQGDL9ZWbaMlZtAhAYPviRl0SMDgZ27GTMng8wQMfu/mJz6LGXyGgMEBHuwxLl+qGhxOwOBg\nN2PmUlVwOAGDN9quwHbvLdMzeJGAwXvdnMsvTpvBywQMzmE0H14kYHC+ezEregb3CRj0xU4jPEnA\noGuPdxoNNzIzAYMYBvRhS8Ag0nYFZpuROQkYBFtyZaCROQkYDEXMmIeAwbBcOJ+xCRhMwYXzGY+A\nwVweXzhf0ggiYDCv65hpGEH+O/sAgC5s9xh3l7OCPlmBAX9dL8isxuiZgAE3rG+ULjJGr4YKWN1s\nfDR/c/AyV1ykZ0MFrOgWvMH2c6WhH6MNcdRaqxPQcLTW/mZsHVaE0425Aqu1rkuxXc8s0eDXtifG\n7CuOJ+7V/1ABuxknxYJj7a5/zzC2z5YRMRtnCzHi7oYxeFlID8ZZgbXW1oZZdcG7ufQUpxsnYEW3\n4LM0jHONs4UIwFQEDPi97Ww9fJiAAQfQMD5vqHNgwFl2s/VOjPEBAga8ZG3Vg4/KhHcQMOB4YsYH\nCBjwXvdiVvSM1wgY8Dm7Ylmc8QoBA07z7U6jN0rzgIABXbgXMx9Fxj0CBnRn+7kt29uPf57ZCBjQ\nqW2WHidK3uYkYEA8eZuTgAGDk7dRCRgwNXnLJWAAd8lbzwQM4Jfk7VwCBvAWNxO1vjfgXt6E7XkC\nBvA5Ny/evyVsz/OBll2oCZ8G6CAPFHGcDvJAzx9na7f/Wz75evffWQfZCSswgAAPNiSf/OHxCBhA\nqsn3IQUMYDSThK220AN/TtyWLsAZboeg8z4MvgIbO88A75Dyyn/wgAHwUymv/I3RAxBJwACIJGAA\nRHIO7GS1/h0EXQcmuxo8uT6qDo8z4iBXKY94KaW11vlB9vyI70agu70zy+VQe/61vEfATrP9/V6f\n1Ha3e7D+Wi9/gX0eZ8RBls2D3vNB3jywfg5y92zb7SO+PYxaa58HWa5eUfV5kPfYQjxNa63/34/+\nj7CEHGQJeUYol2fbs4/iG7sk9CzlOBNZgfG9/v8C+3/CDbLbSupQ1iqhf+vaq+cH/SYB45HdrkK3\n+n/OXY5t+2+f+n+ss3Se2N3G5rkH8wu2EPlGz39+Jeevrl2Uju/SlDsTFl2/OphBzzNp12NUpe/j\n7HkmbdXzI16ujiriIG9+pwe75Ve3B7nc6PkRv0fAAIhkCxGASAIGQCQBAyCSgAEQScAAiCRgAEQS\nMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACI\nJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkCk/wE5IqP0\nqxLs5QAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 27, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "plot(Lambda, Sparsity, '.-');\n", "axis tight;"]}, {"source": ["_Remarks:_ the LARS (Least Angle Regression) algorithm is an approximate algorithm\n", "that only allows the support to increase, thus forgetting about the\n", "events of type 3. Its complexity is the same as OMP (Orthogonal Matching\n", "Pursuit)."], "metadata": {}, "cell_type": "markdown"}, {"source": ["__Exercise 2__\n", "\n", "Test the algorithm for $x_0$ of varying sparsity $s$.\n", "What can you conclude in term of noise robustness ?"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3dbXaj\nyLKF4eCuM692jayrRlbuken+wMaIhCQhv2In77POOsvtklFIRGorE4Sm1+tlAACo+b/eBQAAcAcB\nBgCQRIABACQRYAAASQQYAEASAQYAkESAAQAkEWAAAEkEGABAEgEGAJBEgAEAJBFgAABJBBgAQBIB\nBgCQRIABACQRYAAASQQYAEDS/3oXgDfTNM0/HH1T9nKDyM0iG1n/efxvl9tE7vG02pRilt+nbA2D\nqd3w9t5mN/o/UkNm/1+tDSECzJ2lm4+aON7c6z/c3cjp2Ij8yXqMnd5R/DbLf84/pGwNQ6ra8Js2\nO93a0Q3CIZDZ//dqwwYB1k7KO6wGHbyZ99zeyHoLDDyEnDT8huf+L1LboxBgTe2uv+3+U+KiRMoC\ny24NkbsIX3c293JaWGIxu8M18b4gIaXhT/d4ZsNv2uxG/6cUGd6mRm3YIMDaWQ+/1+t1uiYe2c7u\nz4k13LhNyjC+NOrCpZXlyQlfpCAqpeEt4U3S7YbfbbOUsiO/PKrzaurcqw0bBFg7m5bd/OvmXVj6\nNnc3clpApt2jVjfegdaoDU6cNvyNne6kx+h/J3jWmoofFTg6Ayp9dXGznfAA8uafdg8sb2qIn4hx\n9FiObnNUSeLWoCX9rL9Ij5Vt+NO7W285HCa7qZze/5dqwykCzLuqb8143wdvGvckQ0AaOw8AIIkr\ncQAAJBFgAABJBBgAQBIBBgCQ9JTPgS0nx3LOCgCMYfAACz81P01mxoUecIHcmbpcygRF+O/8wQPs\nwHzh52Kba/BRkjafVuGB7G6qyHZa2lwVQusJf/jd+XloEp3/iGNgr9fP/xYKewe4Y/PSI/FKBNww\n+AeZd99fbIbz0E8AcoleqWE3tBQfCHqR6PwnLiHOO2UZ4NNEhuERuNoeBvPEAJu9XmQYRhC5LHL6\nnwCKHnEM7AiHxKBu+SqpzZdvpf/5rFqBQEUCq5w5pt/vX47we/cLHVY3GPnJwB3OjwTsftNHTiB5\nfrBoyXnnz561hLjJMzN7/X6xlojB5HyrNbMxCHn0EqJ9RxpriRgJIYSHEJgk5jiaBe9OxVhLRMj5\nQkqpJUTPj/EJ1q9IP0c62r4kBV3juvNnAiXmiL/6hEfIyDBsOA8wu3UW4przR/cEDtPr9RLofHt4\ngH3d5r17yDCsSQzjUJhh4YExxcc1Hm/p9fopQaDzBUrMkbgPtiuKf37+ZOinB+ckhnHodBKm+KDG\n0z29Ipclkuh8gRJzXNoHTMUQkhjGoUiAKT6cIblKr/CuJDpfoMQcV/cBUzFsSAzj0Fw2a4Zu9U2v\nlOvBSnT+00+j39h+0vnfn/3MmcmQ83rXuxx88Z9eKgQyNsftNxFHU7Ghny3skHgfGhIt+wl20svN\nsuH7LQVaiBnYvqOpGPMwALd1TK9pGvCyeQIZmyP/TcTuVGzo5wxvJN6HhkTLHlvf9Lp6PxItxAzs\nxO5UjHkYgEt6pdeQE6+FQMbmKPgm4m0qxjzsMSTeh4ZEyx5Vx/Rau3Q/Ei0kUGKOsvsgXE4c+smD\nmcgwDomWPSQP6XXjTiRaSKDEHDX2wWYqNvTzB41hHBItezxd0qvIifISLSRQ4kb8a2fnfw2vz124\nhoTvycQYJIZxSLTswXRPr5x7kGghsZM4dr9APfzX2t+HtEms8MtZADxc+/Qa+3yNXQIZu7b77Ufh\nv4Y/1KqHqdjoJN6HhkTLHkaX9FrLvweJFhKbgZ2apmnzvM+/qXR3TMUGVrVzMLBteq1nRq9X7fSq\ncw9OCWTsWmQGdvTVtG0e4P7FzaBP4n1oSLTsAeyk189/C0y8VlsWaCGBEtfcBpixnDgoiWEcEi1b\nXcf0Kr55iRYSKHFj9wvUl7g6+qdWtb1dwJ4MG4DEMA6Jli2tZXo1uKK8RAsJlJij/T7YZJgRY+Ik\nhnFItGxdg6WXibSQQIk5uuyDr/ZiKjYEiWEcEi1bVLP0avlVXhItJFBijl77IMwwI8Y0+RnGRx/h\n3/29n7KHN2R6mUgLjXYavRNf+/3P+wsNJ9njrvhH+Nt8fh+hLun1qBPl4wQyNkffNxE/PcdUTJaT\n96Hxj/Cn3wYFtUmvxhOv1f0KtBAzsIp+9j5TMVS2e2Ea5mT1vKXXn1qntHeZeAl1jkDG5vDwJuKt\nEzizQ02n84Dep+yvV/wTkLZ3YKx75w9sm14//yE/8VoVINBCAiXm8LMPWE4U5aSFEj/CH/8linhC\neplICwmUmMPVPmAqpshPC23ONpwLC+dqy42dlD2YJb0qHfTyEF0ziRYSKDGHt31wlGFGjHnlrYUS\niZbt3HPSy0RaiJM4mnrrB87sAHS0TC9OlE8kkLE5fL6J2HQqVwH2zGcLnRIt262q6eVq4rWQaCFm\nYB2su2Ka+FIxwLUHppcKgYzN4flNRDgQmIo55LmFIkTL9qb2R5Urf91KFokWYgbWzWYeZny/M+BJ\n1fRaf0tzie09lEDG5vD/JmK3ifl+Zz/8t9Au0bL9qPphL4llQ4kWEigxh8Y+OMswI8b6kWihkGjZ\nTjRLL8+7SKKFBErMIbEP7LinmYp1p9JCG6Jle1AvvSQmXguJFhIoMYfEPpilZJgRY80JtdCaaNnd\nkV4LiRbiJA4v1h9dfAszzuwAmmiTXnxIuSCBjM0h8SZig6mYK4otZLJld1QpveQmXguJFtKbgU3f\nLv2TkPD0+q/fMxUD6iC9RAlk7FrityIt/yTxJmJX5FQlpmItibaQaNldNEgvxV0h0UJ6M7C4eQbm\n/3k/dTQPM6ZiQDk10osPKTcj9lqf+LV+6xnY/Buth7mIDwNOsq9KunnGeBtXW6X0WtPdCRItJFDi\n2o0A03qAofQMM2KsAtEWEi27pZ9L9LJsuEeihQRKXHtggFnCqGAqVo9oC4mW3Uzx9Bpm4rWQaCGB\nEjfChZ3IgqHEPkhxKcOMGCtHtIVEy26D9Eoh0UJ6J3G8vq1/c/RPw4ic0/F1A87sABJUTS8+pNyY\nQMbmkHgTkS5lhZ2pWFmiLSRadm1l02vIiddCooX0ZmBPdjoPM6Zi44p/Tl/98/sNkF7jEcjYHBJv\nIq5KPNOJMzuKcNJCkdOX7P1T/Ee3ebh66TXq0yzRQgIl5pDYB/ewnNiGkxY6Pf82clIuCqbXcyZe\nEi0kUGIOiX1wG1OxBpy00I0Am3/wUHxfldJr4OdVqHlcDM56nLz61HMjw4wYu6JLC20OaG3yafPz\n5pbhbZ5s7nyWDW+QaCGBEnNI7INM6eOKqdgNTloofgws/KWTsvsqlV7PWTZck2ghgRJzSOyDfPcy\nzIixBH5aaLOwE08sP2X3UiO9HvWMSrQQp9GPIOX0+q9bcpK9rM3n9DcvLv5fa1oivR5CIGNzSLyJ\nKOXSYGMqlki0hUTLLqJIej1z2XBNooWYgY0jfR5mTMUwqOLpxdWhPBPI2BwSbyLKurrowZkdcaIt\nJFp2phrp9VgSLSRQYg6JfVBcToYZMfZOtIVEy86Rn14sG65JtJBAiTkk9kENN95FMhXbJdpComXf\nVja9nvTMHZJoIYESc0jsg0oyM8yIMTOTbSHRsu8hvWqQaCGBEnNI7IOqmIplEm0h0bJvmH5PLBvW\nINFCAiXmkNgHtTEVyyHaQqJlX1UwvR7wbF0j0UKcRj++S6fXf/0JJ9nDPdILAhmbQ+JNRBv3hitT\nMdEWEi07XU56sWyYQqKFmIE9xY15mDEVg0ul0osPKasTyNgcEm8iWrq9bPLYMztEW0i07BQF0wsR\nEi0kUOLG6ZetnX7xxMMVyTB7TIyJtpBo2adupxfLhldJtJBAiWsp34pkfK3fmZz3oU+biom2kGjZ\ncUXSa7hnpRaJFhIocS0eYOF3q0vsgy5KZZiNHmOiLSRa9pGcjyqTXvdItJBAiWs3Amz+QethtkGG\nxUk3j8SrTyKWDbuQaCGBEtciATa9d+vut9ZiI/PN6UNiTLGFRMsOsWzYi0QLjXMa/fr7av0/707c\nO7f+5885yR41kV6IE8jYjXBhJ5yKcQzsKqZiR0RbSLTstXvpxbJhKRItpDcDW8+0lt9sbtC8KHlM\nxeBKfnrxIeUnEMjYHBJvIvzIX3gZ7yR7Py10dFLJ7u/9lH1DkfRCJokWEigxh8Q+cKVshpl+jDlp\noaOF8fUHH8dYPL+RXiwb1iDRQgIl5pDYB94UeSc7zFTMSQvFA2y2e1Kuh+LTZaaX1GP1S6h5XAzO\nepy8+sgpnmEmG2NOWigSYKc/qyC9XJFoIb2TONBA5jkdXxvhzI67pne9y2nhanpNE+kFZmA4VuoF\nQnoq5qSFBp+BXWw1oqsBiRZiBoZDReZhxlSshPkaabPllI3d3+shvXCXbNOnER7VbhR8vVA8s0O0\nhWTKvtJenG3YkkQLCZSYQ2If+Fcpw0whxkRbSKPsu+nl/5ENQKKFBErMIbEPJJR9+RCaiom2kEDZ\nyS3FxKsLgRYiwHDJA6dioi3kvWzSyz3vLWRmBBiuetpUTLSFXJd9K73cPppRuW6hbwIl5pDYB3Lq\nZZj5izHRFvJbdlr3MPHqzm8LrXAaPS4rdXr919Y4yf45SC8UJZCxOSTeRIgqvrzjcyom2kLuyk4O\nJZYNnXDXQnuYgeGmsvMwYyo2sFV/TP8e5hJXh8JVAhmbQ+JNhLQarziuzuwQbSFHZb+n19EOZdnQ\nG0ctdEygxBwS+0Bd7QyzrjEm2kJeyr6eXh6qhvlpoSiBEnNI7IMBVHoB8jAVE20hF2UnpBcTL7dc\ntNAZgRJzSOyDMTTIMOsRY6It1L9s0ktc/xZKIFBiDol9MIx6C0Edp2KiLdS57IvpJfgEj0+i8wVK\n3Ih83XX4TxL7YCRtMswaxphoC/Us+yy9mHhJkOh8sdPo5+d0tvtNtZF/QgOv18/r0eas6Nwtc5K9\nBNILDQlk7Frid86uv/RP6wEOY5jlRNEW6lP2lfQSfFKfRaLzBUpcSwmwzW3mH7Qe5hjUlxOlm6f1\nq8/mo8rhjJmJlxoCrLx4gC1fsr57e7RX9R13m6mYaAs1LZv0GpFE54sdAzvl/xl/lOKXm3rb+OpV\nkkNi3VxJr/UhUiCfWIDNJ2jMNuuEy//PelaJFTKslEhvd+v5aHpxbUPUJjBJzCExC36CZmuJVno5\n0UkLRRbPd09ZalF29JQNlg3VOen8OIESc0jsg4eo/Wa80iExJy3kLsCS08vBk4c7nHR+nNgSInRV\nXUu0hy0nbmwW1de/rHR/Pz/+y7LhaISOwhBgaIcMSzS9i9/y6KP98y9rFPfzY5Be7wWUv3M0UKtz\nKiDA0BQZluL1rnc5K2npxdmGaENglTOHxDLuA7U8HmZ5h8T8tNDmg9XrQ1/r3y83Llz2+9uNdXox\n8RqSn86PECgxh8Q+eKzlha/SLipyWodoCxUu+/h0edJrVBKdzxIiuql02d+f7Q+xnNhfWnqxbIj2\nCDD0xCEx7w7Si7MN4QEBhs7IML+O02uN9EIvAqucOSSWcWGOT+sQbaECZSekl+ATg1QSnS9QYg6J\nfYBZg1fGG6d1iLZQbtl76cXE61EkOp8lRHhRey3RWE5MRHpBBAEGR8iw/s7Si7MN4YfAJDGHxCwY\nG43XEi26nCjaQjfLDtLL/rxtRPCZwE0SnS9QYg6JfYBQmzMFUg6JibbQ5bLD6DLS69EkOp8lRHjU\nYC3RWE5cnKUXy4bwafwAU/leAGyQYY0kpBfg0/gBZlJfb4O19Rt/MqyKaHox8YJzAqucOda5NfYj\nHVvH0zokjgSEkso+Sy88mUTnC5SYI5x4jf14B9brtA6JYRw6L5v0QpRE5w+1hJj4DbYsJyrikFhJ\nx+nFsiGECGRsovX7hfC7/kLDPPBHaT8Psz+SrXL49jn4XkozJl7YITEDEygx0dUAWwzzDDyEn4+I\nebb/6sOyIZIRYE3dDjAgTnGM7Lz6RJcNgQ2JAPtf7wIAVLa7bGhmf/y/QAExjw4whi9OaU/ig+LX\nB71of6gTmCSmW15r1muJu7cc6VGjKomFFAvqnKZpU/Qy8VI8pIf2JDp/qBlYytNdY5c02NNtmokH\nomj3Xdr6wf+sGVZLr8ZP+MB3N/BDq2GoADslvauAXTvnK33/3CC6gI604/eU9gEM+CAxRrZvpafJ\n3s7X6FAS1Pnv/KcE2NgPEw+xeUO2Pei1/s/fkzHrwugGDzDgIdQPZgA3DHUtRADAc/CuDQAgiRkY\nAEDSs06jLy686GLBGW3kc9n17qXeaS+bi1UWv5fNl5dy/k4bVYfA5o7W2x/s7qz+AFnf0fJzs0dX\nCQF207oJdq8jXMTSXvMrcu17sfewrJQutR9I1bvAos0+XWvQqL3urvGTeXT+qtxgYQnxptfrVXtP\nt+mkka6LwbeVttRgCGzurtl9Nb679rExzEhhBuZdm8s7Vd1+M3yHzvCW1Yg299XgXrpYL+30riUL\nMzC/5ndJDTqs6lvp+VVg/f+VqA9FxDUbDos2U8xmA2Qx0kghwFwbYO71+maVD4BX2jL8aHlF3TZ3\nZK0GyGKwkSI/heyr3ilYuxcN4izE+PbX29Q9sUpLm7MQw+HAWYg5d7TZuO5gIcAAAJJYQgQASCLA\nAACSCDAAgCQCDAAgiQADAEgiwAAAkggwAIAkAgwAIIkAAwBIIsAAAJIIMACAJAIMACCJAAMASCLA\nAACSCDAAgCQCDAAgiQADAEgiwAAAkggwAIAkAgwAIIkAAwBIIsAAAJIIMACAJAIMACCJAAMASCLA\nAACSCDAAgCQCDAAgiQADAEj6X+8C8GaapvmH1+sVv0HkZpGNbP7p3t0d1XC6tfVmT2vDAy2NkdIM\nKWNhs9nwb5t1/tFtaPscBJgjcysvQ+WooeONvv7DzUY2/2TvI/DS3YVDN3K/m/Iivz997BjVujES\ne8nS3irF/7ZZ5+/eI22fiQBrZz2WTt+pabXyZuyljOGj7cT/HIpSOj/eGLcdbbZUj13t/KONFCzp\nUQiwpjbzofCfLk1iXq9X+irK7hZ2/zNy4/DtZ/q9nNaW+DYWik47/56csXDawzU6P+xt2j4HAdbO\neoy9Xq/bzZr4ji9SwPzD5j9T7i78ZWQuxWjErFTn72559+f0v23W+UyzaiDA2gnXvtfudfald52R\nQ2KZlkRcb7ngG21Iq9H5u/eSuNlenc8QKI4Aa2fT2WE3b96ozj9cPYgd2Uj4Tyl3l3joLtx4+nrm\nbm0YxmnnR26fvsKW0jaRgRD+ecHOn2+wWfCk7fOxzuNd47U4lv7gU9XOpO1FsdsAAJK4EgcAQBIB\nBgCQRIABACQRYAAASQQYMILiF2EC/Bv8c2Cf06eZfbw+OtcBVBOJrrn/jSGAQQ0YYPGP4gJXOf+o\nye4lkZboWv/nL/vVrizoc975NliAHQXVX/tb/B1opU8+stmqW76x2ZHe/fy1v9ZwNubh08Hda+he\nwO0aJDp/qACLXJ2TtRQ80NLt6wkZ6+oYRv93B8Vtrq05v+sMMYARp3WRupQrNW/WFY1RgGMe5o6n\nBEq8avcLUsOhu8YwxhGJYWxXvmpgMxZofuyS6HyBEq96D7Cf3y8PNBJmDGZsSAzj0GnZ61FA2yMk\n0fkCJV51FGC2yrAZSYZTEsM4lFI2UzFESHS+QIk5wvM5dh8uSYYjEsN4x9L6V2KMbsdCovMFSsyx\n7IP8JDOG9yNJDOMdm4+F/bWPj8NHQYYhJNH5AiXmiCwnziKPnmkZTGQYh6ZpCouOxBjLidiQ6HyB\nEnNMv99S6/V7fzZmJBkOSAzj0OfnV5d/bC6+cfxYyDCsSXS+QIk5NgE2i8SYkWR4JzGMQ3OA/cy3\n4qczrf+Q5USYmUjnC5SY4+cYWJBkc4zZ9QnZjDB7CIlhHNove/djJQEyDCbS+QIl5th+urPohGxG\nko1NYhiHDssmw5BGovMFSsyxuw8iMWZ3J2RGkg1KYhiHYmWnLSeSYQ8n0fkCJeaI74PIuqKRZDAz\nkWEcil1Kaj48tj65gwxDQKLzBUrMkbIPKsWYkWRDkBjGoXjZOxlm+21Nhj2WROcLlJgjfR/EY8yS\nPwq9iyTTJTGMQ+fXQkzOMCPGHkmi8wVKzHF1H+x+bux9gzt/dekpPAozXhd8khjGocQAmzEVQ0ii\n8wVKzHHzq0jrx5iRZCIkhnEo6WK+y4edP14pZ3aQYY8i0fkCJea4vQ9OM+x7+9vf3Lg3kswziWEc\nSiz7LcPs/CR7Muw5JDpfoMQcmfvgdowZSTYKiWEcuhpgtpthttPHXHHqISQ6X6DEHEX2QeMYM076\n8ERiGIfSy95OwowMg5lI5wuUmKPUPkjMsO873f7mdgkkWXcSwzh0qeydDDNi7OkkOl+gxBxl90Fm\njBlJJkhiGIcKBJiRYY8m0fkCJeYovg8uZZiVjjEjyZqTGMahq2UfZphdO7PD6MNRSHS+QIk5Ku2D\nqzFmRdcVZyRZGxLDOFQywIwMeyKJzhcoMUfVfVAkxowk801iGIdulJ2ZYUaMjUWi8wVKzFF7H9zI\nMKsTYzNOxC9OYhiH7pV9IcOMqdjgJDpfoMQcbfbBvRgzkkyBxDAOVQmw703v/PL9voixAUh0vkCJ\nOZrtg9sZZjVjzEiybBLDOHS77PsZtvZ6kWHqJDpfoMQcjfeB2xizgyTjZeWUxDAO1Q2w1d2cb9D+\nLj/Tb0IkOl+gxBzt90FOhn1toXmS8bISITGMQzllX8uw93s937j9pd8kSHS+QIk5eu0DYmwYEsM4\n1CfAgiJi/yr4rD6KROcLlJij7z5Yx9iNDDNizAE/w3j67oZNPdOqS5Z/yiy7WIYtCDM1fjo/QqDE\nHN33Qf5U7Gs7zU/0IMlm3VsoLGNT0m6F7gJsjTBT4KTz4wRKzOFkHxBjury0UDTA5h9OU+2Sqhk2\nN9uH/YrdyMHT/mROOj9OoMSNo4UU21tL8bMPSmXY19ZKX5hqjXXFNSctlDID29xm/sHvPGzz9ZiE\nmRtFmqcNF4MzXWQYJ/6mr9oxZjUnZM+MMSctFO/89NtcUjvAvu4l/MQYa4wOOOn8OIES104DbP6h\n7DAuq2yGGeuKlTlpoaPOT/k5R5sMs6OLd5Bk/Tjp/DiBEtcSZ2DrH+Z/8vYwi8eY1VxXfOb1gr01\nz6aeeJPLBZidTvoJs7YIsPIS32PuHhXwpkaGWacLU9m4Sea5hSIKlt0ywyzxOoqRMBPcWT5JdL5A\niWuJJxNLBNhMMcbsSUnmv4V26QbY150mXkqRaVk1Ep0vUOJG5ECXxDGwUKUM+9p4zfMVZ2MfJ5No\noVDZsj1kmJ02FWFWlETnC5SYQ2IfzBrHmDVJsgFiTKiF1gYIsK+7vndVe8Ism0TnC5SYQ2IfLKpm\nmPWLMVNOMq0WWhQvu2OGWc4XjJFkd0l0vkCJOST2wUbtGDPWFa9QbCEbLsC+asj8jjHC7AqJzhco\nMYfEPgg1yDBjXTGNagtVKHuEDFsQZmckOl+gxBwS++AIMeaBaAuNGmBflRT8umeS7IBE5wuUmENi\nH8TlfydL6h31WFf0H2OiLVSpbD8ZZmVjbEaYrUh0vkCJOST2wak2U7Gv+3ofwjWePK0YE22hJwSY\n1ciwxeM/Ky3R+QIl5pDYB4mIsS5EW6he2d4yzKrGmD03ySQ6X6DEHBL7IF3LDDNizMxkW6hq2Q4z\nzGrHmB0nmWB7pJDofIESc0jsg6s6TsWsVYyZmyQTbaE2AWZPy7BFm3OfupLofIESc0jsgxv6TsXs\nSTEm2kK1y/Y5CZt1jjEbJMkkOl+gxBwS++C2IWPMnK0rirZQg7LJsB8jTsgkOl+gxBwS+yBH4wyz\n58WYaAu1Kdtzhln7GLOhkkyi8wVKzCGxD/IRY/WItlDjADMybG2IGJPofIESc0jsgyLaZ5g9I8ZE\nW6hZ2f4zzHrFmGknmUTnC5SYQ2IfFESMFSfaQi3LJsNOaJ7rIdH5AiXmkNgHZXXJMBs3xkRbaP5y\nVzJso2eMmdiETKLzBUrMIbEPaiDGShFtofDbyWtzfkLHWucYM40kk+h8gRJzSOyDeppdCHh7v/1i\nzEq/Hom20LTaB2RYqPvJQWbelxYlOl+gxBwS+6CqXlMxGyXGRFtoen/2OaFjV/+p2MzlhEyi8wVK\nzCGxDxp4QoxZnbfVoi1UKcBOJ3afn9OvX4d/7vOZdB1j1i3JJDpfoMQcEvugjY4ZZsoxJtpC096r\nYcoD2f3DBvo+yS5WFBc+JmQSnS9QYg6JfdASMXaVaAv1yqHiWj75vmLMOieZROcLlJhDYh80tskw\ncxBj5vhkRdEW6j6ROj2hI6fCutcpdrKiuOi0tCjR+QIl5pDYB130nYqZToyJtlDVANs8IUfLlfdO\nSrxXefF95C7GrPWETKLzBUrMIbEPOuoeY+b+o2OiLeR5CfHe83npERXZZe5WFGetJmQSnS9QYo7G\n1yNQ1H1F8asMrzHmfxiHn1n2nF67qkZa5u5zGmNWfULmv/PtIQFmZNgZYuzoVcn5MF6Xt/m5X1G5\nbj/hKY/63sY9rigu6iSZ886fCZSYo+MwVnxiPawomrNPQDsfxkMGGJzw3Pkz14Mz3xjDuPE+IsZs\nFWMEGB7Lc+fPXA/OqyodDAifIulXh8Q97jbGrO264vAB9nq91oeKwz+cb/Cz5e/G2LREZLn+3qHo\n+DM/fX5u6/z4SNnsx8fHf//9988//5jZZ7CR5b7fN70t43P69LWKaGaba5/+senf7//487K0UXP2\nuN1xPTgv2R3JzZKmwYlVzYQvVV+/75Rh5iDGqtxTCTkBdno2/G5XHwVYS5v0SowuM/tY3fJ2enkT\nHsD+8uen8tMHsdcvAukgUGKivgGGgTkfI7tTn0jn78dS2rUTu6fX7eiyVXqNEV0puWW3omv+E+dr\nD7P/9S4AQJb4q0z8XxNnXR7cXjOcnaSXTnQl5pZlRJcQAmzL7QC+bf0lF/Z9WYTMl4MfOiP/niEn\n8fcu9dtr+lUqukw5vQrm1tcGBR70OQLMbMTQWoRf0UR0PfJTmtYAAAsCSURBVMp6gfEojN32f36v\nnqSX+x4+yq3X79ft+ZP7B32BwCpnusSzEEd6yHGn6XUzumysQRAlcSQgVPwjve2nX/m9Gksv3w0c\nyS3LWPq79IcSnT/UDCznYMB4NulFdD3N0ZTL/0Ao0quK6XWaW9Of4J8qRJeQoQIswv+gLatKenkd\n9jhS7IuYW02/Sr3NOkwvlz0czy3L/hiJywddxiMC7FHptTll45f9tdUAZuIFn4odmj1KL38NfJpb\nlj1z8vegCxs/wB6bXky8UETt6VfB6LLd9PK0fHZ4MuHOlU2CG9ydcl36Wy3jB9hzRNKLiRd8KnZo\n1szC9HLzQp6eW1Y6ui79uRwCbBDr9CqzbEh0YaX49KtsdNkmvf77z0MDZ+aWEV1nCLARbNPrGxMv\nZIq8BN/fZunosjC93u6gdfdeyi3Ly55n5taCAJNXOL2ILlRT9nDX4jC9un4P0drR/JXVwkwEmLCd\nEw6/3XldcHPAAE6UPX2jxsTLjtKrYeveyC0jugohwFQx8YKQMqcUvft4385XerVq3VK5ZZxemIEA\nk1QyvYgu7Ck4/VrSq1R02e7Eq0nr3sst40BXHQSYnt30Ys0QBRU8d6NFetXv25QPHR/+LdFVDQEm\npkx6EV1Ikzn9Kp5e22XDo6+mLCQnt4wPI9dHgCkJ04voQnGlFg8Lp9c0ffzzz/oXVdNrN7rSn5CC\nF4xP/9sHIsA07J5weO11gehCQ1XTq3F05eSWEV01CXzjSw6Jr7Q5VWDZkDM17hJtodtlF5l+lUyv\nhhMvP9HloeMkOp8ZmHe5y4ZEF9oKP618d0OTmbVJr7LRxZSrGQLMtaz0IrpwUf70q9jnvYL0IroQ\nIsD82qQXEy84V/YLEGqnl5PoYmjmIMA8Ck/ZYOKF2jKnXwXSay+6rHR6dTm9kClXJQSYO+HE60Kf\nk17oITe9Vn1bb+JFdI2HAPPl/rIh0YUMOdMv5+mV+WFkKxddDMriCDBHbqYX0YV+stLrILqsRHrl\n55bdCiGmXC0RYF6s0+vjI7nhSS9kuz39uv95r/e+LTjxqpRbRnS5RID1N31+/rVfy3+mHvMiulDC\n7ev23kyv4GW+SHoVyS0rF10MxzYIsM4+P6e/q/8kvXBk+t7pm+sjTKtmyLl0wrU5Son0yl82LJVb\ndiu6mHJ1R4B1w8QL6dbX9Qmv8XM7t+4tHt5Jr6Bv374SpVB03fwAAAe6ZBFgfdxJL6ILB+ZJWJsr\n1/VNr4K5FdZlRJcaAqy1efyTXihojq7NzOw01W5Mvy5f53A3uq6nV8Glwt26jNXClWn30bpEgDV1\nJ72IrufZvILEp1ZH/xr/q8zvXE6afmVPvIrnVliUEV2B5f1Q70LO6QXY0aFss2JHsyvZLBvaaXpx\nbtNTpXdv/nde3Jh+nafX3ruu9PSqkVsHRaXeOOVP0J5YgMUPZZvL3DImXsj2er02b93m/g9/n6LU\ndy4fbD028bJoepU9xHVUkR2Pp6NZB+PPJ7EAOxWu+7c8vr1f0tX0Irp88LaEcnTmoaNzN+4e8ao0\n5QorsoPxRG6JGi3AwqPZ3aPLSC9NQkcCLqly7sZeUJxOvOrl1m5R4WCK7FtGngTXAXbpUHbKDVq6\nnF5EF+orf+7GUXQdT7xq51ZYVHp0Mey0uA6wS4GUfzS7lPXb1aT04mQNNFfg3I2L0dUgtw6Kiv1r\neBsIcR1god1D1jlHs4vbnXhZenoxklBNsXM3EqLLVulV6dSMhKJi/2qMNn1iAWZ74eThcNfs2rIh\n0QX3dqZfZ6fIW3TiVeW8x+PBRG6NTS/AfMpdNmRIobKcczf204voQm8EWAHX0ovogqjk9ArXDGuk\n1+6aIbn1KARYls3pxSfpxcka6KTA9Os9vfpOvI6iixH2NATYfUcTL0tJLwYWhKSlV5c1w91fMrwe\nwsup55VUOrc+a+I19BM+Hj8fz7jk7aJrmdOvVQMffQVlx+jaENxXTkl0PjOwy0gvqLj3seW3Dp8m\nC3JrNqdXr8Nd2ztlYD0SAXbNemy/Pj4+P1dvTjfpRXTBkxu58s/v32F0dT9T4+0eGVXPJjBJzFFw\nFryZeMXSi+gaiMRCSmiaJvv36+f0XPmIfk9KyzM1uEphdxKdzwwsyWbiZWakF1xLTq94aM2aRRe5\nhUsIsBPxiZet04vogg+Jh74i0fX5339fP723caX0Yp0Q9whMEnNkzoKZeEFiISU0J81RwITR9fn5\n+dbtv36ZvbVxjeji1AzPJDqfGdi+cOJlR+lFdMGhP/uXBt2Nrs1vNulVPLo4Gx6lEGA7UtOL6IJv\nF07KWD74VTq95lFydJGnEMMI6QiwrXDZ0EgvqPnnn39ih7ji37A8fxX1XnTd/nrqlOsWMoZwlcAq\nZ45Ly7i7Ey8L04th9yQSRwI2jqLrKLfepl9hev0p+fBfe2//jGHkj0TnC5SYI3EfpESXmX28Xe9w\nvunIzx5MZBiHPr7WvT9Pb7lZcph+T6eJlfN8EF0qJDpfoMQcKfuAiRciJIZxKL3s6fPTfn2k3DLz\naSC6tEh0vkCJOSL7YApPvlotvLyl12biNfQzhg2JYRxKeuvW6mxAokuRROc/9CSOo1nXjPTCwGK5\n9e9U+6JQDCAU9LgAi0eXRdKLkQdlh7n199M+vxq96nWhGEAoTmCSmOPtW5GuRJet02vopwhxEgsp\noZ/TCP+89j+D9ffz64fPX0Z0ISDR+QIl5pimyf7+DX/PxAuJJIZx6Gi+9bLJzKbVoAjHQtn7FXzy\nYCbS+aMvIQbptTtc99PL/c4D0s3RZdXSi+hCe6MH2MppdNmSXow8yJvmi11EblEvvRhAaOMRAXY0\nUJl4YRjhgs/rZZv8qjH3IrrQ0eAB9vr16/Pv95vRxTR9vq8sMvGCrmn/eNdrmszsZfNELOGSHNfv\nN7hLBhDaEjhMl2OZYy2X0vgMTul4+1qvZJWOcLLZqlu+sVmJQ9kW1PmTLsvZhou0627YKpDCaVbK\nxMvDU9e9hu4F3KhheUvUvfJTg8/AFp+f00eQXveiC1Cwd/5tcnR9beLgCNr69+5f4nDZnFsHM3tf\n+r87qGpzjsbiV3hZXuCAnzGyeU15n3K9z8A2a4bXOj798Qq8xuE2P51/ZPAZ2DzH2p5q+OF/vwA7\nbjTu18katTqegYSeBp+BAQ/h4VgL0BhNDwCQ9H+9CwAA4A4CDAAgafCTOGpYDjaU+rTEZjuVNltw\ny7a5zH/pgudNOX96scvPs118nF666809dnke5jvtW0NtBNgF61fYzSt4Zk+sR1qNzc55UGrLtZ+H\ngptdP6XFn17s8vBsLy1atj/T77rS0LtURvcaGmAJ8YLX61Vjx1dqJq1rZCwbr/HxyXmzwwxazzw8\nyR72NS3XBjMwL5Y3SmW3WXBrDdS4BEC4hILaajSzECdznc0keEgEWH/1+qx4HsybWv6/bM1jj7SH\n6P6iuWnRXmX0NeqCYYglRBeKd1ilofv6ZqVrfuxrzXj6vlzWa1E4NHI4V1L27Kbw6nachShxkidC\n9Zr5XjFPPgvRQw0NEGAAAEksIQIAJBFgAABJBBgAQBIBBgCQRIABACQRYAAASQQYAEASAQYAkESA\nAQAkEWAAAEkEGABAEgEGAJBEgAEAJBFgAABJBBgAQBIBBgCQRIABACQRYAAASQQYAEASAQYAkESA\nAQAkEWAAAEkEGABAEgEGAJBEgAEAJBFgAABJBBgAQBIBBgCQRIABACQRYAAASQQYAEASAQYAkESA\nAQAkEWAAAEn/D/hShpK0aERYAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 28, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo2()"]}, {"collapsed": false, "outputs": [], "prompt_number": 29, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}], "metadata": {}}], "nbformat": 3, "metadata": {"kernelspec": {"name": "matlab_kernel", "language": "matlab", "display_name": "Matlab"}, "language_info": {"mimetype": "text/x-matlab", "name": "matlab", "file_extension": ".m", "help_links": [{"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md", "text": "MetaKernel Magics"}]}}}