Python-for-Signal-Processing-Unpingco/Compressive_Sampling.ipynb

{
 "metadata": {
  "name": "Compressive Sampling"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Compressive sampling Overview"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In our previous discussion, we saw that imposing bandlimited-ness on our class of signals permits point-wise sampling of our signal and then later perfect reconstruction. It turns out that by imposing *sparsity* we can also obtain perfect reconstruction irrespective of whether or not we have satsified the sampling rate limits imposed by Shannon's sampling theorem. This has extremely important in practice because many signals are naturally sparse so that collecting samples at high rates only to dump most of them as the signal is compressed is expensive and wasteful. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "What Are Sparse Signals?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's carefully discuss what we mean by *sparse* in this context. A signal $f$ is sparse if it can be expressed in very few nonzero components ($\\mathbf{s}$) with respect to a given basis ($ \\mathbf{\\Psi} $ ). In other words, in matrix-vector language:\n",
      "\n",
      "$ \\mathbf{f} = \\mathbf{\\Psi} \\mathbf{s} $\n",
      "\n",
      "where  $ || \\mathbf{s} ||_0 \\leq N $ where $N$ is the length of the vector and  $|| \\cdot||_0$ counts the number of nonzero elements in $\\mathbf{s}$. Furthermore, we don't actually collect $N$ samples point-wise as we did in the Shannon sampling case. Rather, we measure $\\mathbf{f}$ indirectly as $\\mathbf{y}$ with another matrix as in:\n",
      "\n",
      "$\\mathbf{y}  = \\mathbf{\\Phi f} = \\mathbf{\\Phi} \\mathbf{\\Psi} \\mathbf{s} = \\mathbf{\\Theta s}   $\n",
      "\n",
      "where  $\\mathbf{\\Theta}$  is an $M \\times N$ matrix and $ M < N $ is the number of measurements. This setup means we have two problems to solve. First, how to design a *stable* measurement matrix $\\mathbf{\\Phi}$ and then, second, how to reconstruct $ \\mathbf{f} $ from $ \\mathbf{y} $. \n",
      "\n",
      "This may look like a standard linear algebra problem but since $ \\mathbf{\\Theta} $ has fewer rows than columns, the solution is necessarily ill-posed. This is where we inject the sparsity concept! Suppose that $f$ is $K$-sparse ( $||f||_0=K$ ), then if we somehow knew *which* $K$ columns of $ \\mathbf{\\Theta} $ matched the $K$ non-zero entries in $\\mathbf{s}$, then $\\mathbf{\\Theta}$ would be $ M \\times K $ where we could make $M > K$ and then have a stable inverse.\n",
      "\n",
      "This bit of reasoning is encapsulated in the following statement for any vector $\\mathbf{v}$ sharing the same $K$ non-zero entries as $\\mathbf{s}$, we have \n",
      "\n",
      "$$1-\\epsilon \\leq \\frac{||  \\mathbf{\\Theta v} ||_2}{|| \\mathbf{v}  ||_2} \\leq 1+\\epsilon $$ \n",
      "\n",
      "which is another way of saying that $\\mathbf{\\Theta}$ preserves the lengths of $K$-sparse vectors. Of course we don't know ahead of time which $K$ components to use, but it turns out that this condition is sufficient for a stable inverse of $\\mathbf{\\Theta}$ if it holds for any  $3K$-sparse vector $\\mathbf{v}$. This is the *Restricted Isometry Property* (RIP). Unfortunately, in order to use this sufficient condition, we would have to propose a $\\mathbf{\\Theta}$ and then check all possible combinations of nonzero entries in the $N$-length vector $\\mathbf{v}$. As you may guess, this is prohibitive.\n",
      "\n",
      "Alternatively, we can approach stability by defining *incoherence* between the measurement matrix $\\mathbf{\\Phi}$ and the sparse basis $\\mathbf{\\Psi}$ as when any of the columns of one cannot be expressed as a small subset of the columns of the other. For example, if we have delta-spikes for $\\mathbf{\\Phi}$ as the row-truncated identity matrix\n",
      "\n",
      "$$\\mathbf{\\Phi} = \\mathbf{I}_{M \\times N} $$\n",
      "\n",
      "and the discrete Fourier transform matrix for $\\mathbf{\\Psi}$ as\n",
      "\n",
      "$\\mathbf{\\Psi} = \\begin{bmatrix}\\\\\\\\\n",
      "e^{-j 2\\pi k n/N}\\\\\\\\\n",
      "\\end{bmatrix}_{N \\times N}$\n",
      "\n",
      "Then we could not write any of the columns of $\\mathbf{\\Phi}$ using just a few of the columns of $\\mathbf{\\Psi}$.\n",
      "\n",
      "It turns out that picking the measuring $M \\times N$ matrix randomly according to a Gaussian zero-mean, $1/N$ variance distribution and using the identity matrix as $\\mathbf{\\Phi}$, that the resulting $\\mathbf{\\Theta}$ matrix can be shown to satisfy RIP with a high probability. This means that we can recover $N$-length $K$-sparse signals with a high probability from just $M \\ge c K \\log (N/K)$ samples where $c$ is a small constant. Furthermore, it also turns out that we can use any orthonormal basis for $\\mathbf{\\Phi}$, not just the identity matrix, and these relations will all still hold."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Reconstructing Sparse Signals"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that we have a way, by using random matrices, to satisfy the RIP, we are ready to consider the reconstruction problem. The first impulse is to compute the least-squares solution to this problem as\n",
      "\n",
      "$$ \\mathbf{s}^* = \\mathbf{\\Theta}^T (\\mathbf{\\Theta}\\mathbf{\\Theta}^T)^{-1}\\mathbf{y} $$ \n",
      "\n",
      "But a moment's thought may convince you that since $\\mathbf{\\Theta}$ is a random matrix, most likely with lots of non-zero entries, it is highly unlikely that $\\mathbf{s}^* $ will turn out to be sparse. There is actually a deeper geometric intuition as to why this happens, but let's first consider another way of solving this so that the $\\mathbf{s}^*$ is $K$-sparse.  Suppose instead we shuffle through combinations of $K$ nonzero entries in $\\mathbf{s}$ until we satisfy the measurements $\\mathbf{y}$. Stated mathematically, this means\n",
      "\n",
      "$$ \\mathbf{s}^* = argmin || \\mathbf{s}^* ||_0  $$\n",
      "\n",
      "where\n",
      "\n",
      "$$ \\mathbf{\\Theta} \\mathbf{s}^* = \\mathbf{y} $$\n",
      "\n",
      "It can be shown that with $M=K+1$ iid Gaussian measurements, this optimization will recover a $K$-sparse signal exactly with high probability. Unfortunately, this is numerically unstable in addition to being an NP-complete problem.\n",
      "\n",
      "Thus, we need another tractable way to approach this problem. It turns out that when a signal is sparse, it usually means that the nonzero terms are highly asymmetric meaning that if there are $K$ terms, then most likely there is one term that is dominant (i.e. of much larger magnitude) and that dwarfs the other nonzero terms. Geometrically, this means that in $N$-dimensional space, the sparse signal is very close to one (or, maybe just a few)  of the axes. \n",
      "\n",
      "It turns out that one can bypass this combinatorial problem using $L_1$ minimization. To examine this, let's digress and look at the main difference between $L_2$ and $L_1$ minimization problems."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "reference:\n",
      "`http://users.ece.gatech.edu/\u0018justin/ssp2007`"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "$L_2$ vs. $L_1$ Optimization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The classic constrained least squares problem is the following:\n",
      "\n",
      "min $||\\mathbf{x}||_2^2$\n",
      "\n",
      "where $x_1 + 2 x_2 = 1$\n",
      "\n",
      "with corresponding solution illustrated below.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "from matplotlib.patches import Circle\n",
      "x1 = linspace(-1,1,10)\n",
      "fig=figure()\n",
      "ax=fig.add_subplot(111)\n",
      "ax.plot(x1,(1-x1)/2)\n",
      "ax.add_patch(Circle((0,0),1/sqrt(5),alpha=0.3))\n",
      "ax.plot(1/5,2/5,'rs')\n",
      "ax.axis('equal')\n",
      "ax.set_xlabel('$x_1$',fontsize=24)\n",
      "ax.set_ylabel('$x_2$',fontsize=24)\n",
      "ax.grid()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEVCAYAAAAl9QikAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtQVPfZB/DvLiwsyE1QLgJmjcZgFAQ0Qi40G1NzoZFo\nm7Y275vR2LFxOv0jr51Ik0zT2iTGzHQ6sTp9YzJJNJ03eW1MjdTbm2oF0ySIdxLBu8h9XcDlIrDX\n8/5hIKLclr2c8zvn+5nJ1IUDPE/PgWd/z3MuOkmSJBAREflAL3cAREQkPhYTIiLyGYsJERH5jMWE\niIh8xmJCREQ+YzEhIiKfhcodgFx0Op3cIRARCWeoq0k0vTKRJEmI/373u9/JHgPzZt7Mm3kPR9PF\nRBQ1NTVyhyAL5q0tzFtsLCZEROQzFhMBLFu2TO4QZMG8tYV5i00njdQIUymdTjdiD5CIiL4z3N9N\nrkwEUFpaKncIsmDe2sK8xcZiQkREPmObi4iIRoVtLiIiCigWEwGopafqLeatLcxbbCwmRETkM85M\niIhoVDgzISKigGIxEYBaeqreYt7awrzFxmJCREQ+48yEiIhGhTMTIiIKKBYTAailp+ot5q0tzFts\nLCZEROQzoWYmy5cvx65du5CYmIivv/560G1eeOEF7Nq1C5GRkdi8eTMyMjIG3Y4zEyIi76hmZvLM\nM89g7969Q35+9+7dOHnyJCorK7F+/XrVPHSGiEjphComBQUFGD9+/JCfLykpwdKlSwEAeXl5sNls\nsFgswQovYNTSU/UW89YW5i02oYrJSBoaGpCent7/Oi0tDfX19TJG5Lu2NsDplDsKIqLhhcodgL/d\n3M/T6XRDbrts2TKYTCYAQFxcHLKzs2E2mwF8925B7tfHjpmxfr0ZTz5ZikcfBRYsUFZ8gX7dRynx\nBOO12WxWVDzBfN1HKfFofX/3/bumpgYjEWoADwA1NTVYuHDhoAP4lStXwmw2Y8mSJQCAjIwMlJWV\nISkp6ZZtRRrAl5cDa9YAVVXAiy8CzzwDhIXJHRURaY1qBvAjKSoqwgcffAAAKC8vR1xc3KCFRDS9\nvaXYswfYuhX49FPgjjuATZsAh0PuyALr5nerWsG8tUUteQvV5vrZz36GsrIytLS0ID09HWvWrIHz\n24HCs88+i8LCQhw8eBCZmZkYN24c3n//fZkj9q/8fGDPnu9WKmvXcqVCRMogXJvLX0Rqcw2F7S8i\nCibNtLm0pm+lorX2FxEpD4uJAEbqqaq1qKill+wt5q0tasmbxURF1FpUiEj5ODNRMc5UiMifODPR\nKK5UiChYWEwE4GtPVdSiopZesreYt7aoJW8WEw0RtagQkfJxZqJhnKkQkTc4M6FB3bxSmTYNeOst\nwG6XOzIiEg2LiQAC3VPtKyoffwzs2HG9/aWEoqKWXrK3mLe2qCVvFhPql5f3XVEpKVFOUSEi5ePM\nhIZ06ND1mco333w3UwkPlzsqIpILZyY0Jnl5wO7dXKkQ0chYTAQgd09VrqIid95yYd7aopa8WUxo\n1LhSIaKhcGZCY8aZCpG2cGZCAcGVChH1YTERgNJ7qoEqKkrPO1CYt7aoJW8WE/IbrlSItIszEwoY\nzlSI1IUzE5IFVypE2sFiIgDRe6pjLSqi5z1WzFtb1JI3iwkFDVcqROrFmQnJhjMVIrFwZkKKxJUK\nkXqwmAhALT3VoQxVVP75z1K5Q5OF2vf3UJi32FhMSDFuLir/8R9cqRCJgjMTUizOVIiUhTMTEhJn\nKkTiYDERgFp6qt7qy1trRUXr+1tr1JI3iwkJQ2tFhUgknJmQsDhTIQouzkxIlbhSIVIOFhMBqKWn\n6q3R5q22osL9rS1qyZvFhFRDbUWFSCScmZBqcaZC5F+cmZAmcaVCFDwsJgJQS0/VW/7KW7Siwv2t\nLWrJm8WENEO0okIkEs5MSLM4UyHyDmcmRIPgSoXIf1hMBKCWnqq3gpW30ooK97e2qCVvFhOibymt\nqBCJhDMToiFwpkI0EGcmRGPAlQrR6LGYCEAtPVVvKSXvYBcVpeQdbMxbbCwmRKPElQrR0DgzIRoj\nzlRIazgzIQoArlSIvsNiIgC19FS9JUre/i4qouTtb8xbbCwmRH5yY1HZsYMrFdIWzkyIAoQzFVIb\nzkyIZMCZCmmJcMXk4MGDyM3NRVZWFjZs2HDL50tLSxEbG4ucnBzk5OTg1VdflSFK/1JLT9Vbasnb\n26Kilry9xbzFFip3AN5wu91Yvnw59u3bh9TUVNx99934/ve/jxkzZgzY7oEHHkBJSYlMURINrq+o\n9LW/1q5l+4vUQ6iVSUVFBaZNmwaTyQSDwYAlS5Zgx44dt2yntlmI2WyWOwRZqDXvkVYqas17JMxb\nbEIVk4aGBqSnp/e/TktLQ0NDw4BtdDodvvzyS8ycOROFhYWoqqoKdphEo8KZCqmJUG0unU434ja5\nubmoq6uDwWDAli1bUFRUhPPnzw+67bJly2AymQAAcXFxyM7O7n+X0NfHVMLrG3uqSognWK9PnDiB\n5557TjHxBOp1Xh6wenUpqquBkhIzXn65FP/5n8BjjwELFsgfX7Bea2V/3/xayb/fff+uqanBiCSB\nfPXVV9IjjzzS/3rt2rXSunXrhtze4/FI8fHxUmtr6y2fEyn1AwcOyB2CLLSa91/+ckB67DFJSk+X\npP/+b0nq7ZU7ouDQ6v4WKe/h/m4KdZ2Jy+XCnXfeif3792PSpEmYN28ePvroowEDeIvFgsTEROh0\nOpSUlOCXv/wl6uvrb/levM6ElI7XqZDSqOY6k9DQULz33ntYvHgx5syZg+XLl2PGjBnYtGkTNm3a\nBADYtm0bMjMzkZ2djW3btg06oCcSAWcqJBKhVib+JNLKpLS0tL+XqSXMeyC1r1S4v5VPNSsTIi3j\nSoWUjCsToiBzu93o7e2F3W5Hb2/vt/92wOOR4PFcPyb1eh1CQvQIDw9DeHg4jEZj//+GhIQAUP9K\nhZRnuL+bLCZEAdTT0wObzYarV9vR2NiOpqZ2dHe7ABgBhH/7v0ZIkgE6nb7/9HdJkiBJHgAO6HR2\n6HS9kKReAHZERoYiJSUWkybFIj4+DmfPxuH11404dYpFhQKLxWQQIhUTkXqq/iRi3i6XC1arFZcv\nN+P8eSu6unTQ6WKh08UhIiIWkZGxCAszDvs9vv66FJmZ5iE/73D0oru7Hd3dNgDtkCQboqIAu92E\n//mfyTh/PhwvvqgTrqiIuL/9QaS8h/u7KdRFi0RK5HA40NjYiAsXLLh4sQ1udzwMhmTExmYgNjbC\n7z8vLMyIsDAj4uKS+j9mt/egt9eCn/70JC5fduO99+7CK69E4YUXgBUrQoUqKiQmrkyIxujq1as4\nc6YG33xjgcuVhIiIZMTGTkRIiLzv0dxuF9rbrTh1qh179iTAao3F6tUu9Jx6Ds4LZ2/Z3jh9On7z\n9tsyREqi4cqEyE88Hg/q6+tx8mQN6utdCAszIT5+JgyGMLlD6xcSEor4+BQUFKQgP9+BI0eu4N13\njZh07jT+6fz8lu1/H/wQSYV4arAAbrxPjpYoKW9JklBXV4dt2w5g584mtLfPQFrafCQm3u73QvL1\n16V++14GQxjuuScNr7+egNRJHr9930BQ0v4OJrXkzZUJ0Qiam5tRXn4aTU1hiI3NQXp6vNwheU2n\n0yE6evBfdw+7veQHnJkQDaGzsxNfflmJCxdciI7OGDDwFtFnL5qx4ZuyWz7+sKEAj//xM6xcaUSY\ncrp1pEC8Ap7IC5Ik4ezZ89i69Ss0NKQhLe17wheS4aSkuPDOO52YMsWJt96S4HDIHRGJaMwrk/b2\nduzevRuHDx9GRkYGfvzjH2P8+PEAgAsXLuCjjz7C7bffjunTp2Pu3Ll+DdofRFqZiHQeuj/JkXdn\nZyf+/e8TuHTJgMTE2QgP9/+pvSMZ6TqTsdq98RcIbbz1bC7XpOl4aMWbKC+/gLKyVFgscXjpJT2e\neQZBXanwOFc+v5/NdfXqVXzve9/DqVOn+j/28ssvY/PmzXj00UcxdepUrF69Gs899xyefvppuN3u\nsUVOFEQXLlzCgQPnEBKSgfT0yXKH43eFvxr+9N8HHsjEXXfV4Pz5s/jooyysXRvZf0U92180kjGt\nTFasWIGWlhasWbMGU6dORW1tLT788ENs2rQJb7/9NhYtWgQAWLduHV588UV4PMo7i0SklQkFlsfj\nwZEjlTh8uANJSXMRHh4pd0iystu7YbEcxrhxqdi6dSqqqnQsKgQgADOTs2fPYvv27cjKysK4ceMw\nY8YMvPLKKzhy5Ag2bNiAvXv3+hQwUbDY7XZ89tmXOHLEhdTU+zRfSAAgPDwSqan3obv7Kv7rv77E\nX//qwKefXr9L8aZN4EyFBjWmYnLbbbcN+vHJkydjz5492LlzJ/bs2TOqZ7bTyNRyHrq3Ap13e3s7\nPv30c9TVJSI9fW7/3Xjl5s/rTMYqJCQU6el34/LlCbBYPsf//m87tm5FQIsKj3OxjamYpKen49y5\nc/jTn/6E1157bcDnwsLCsHHjRpw5cwafffYZCwopUmtrK7ZvPwS7fRaSk6fLHY5iJSffie7uGdi+\n/RCmT2/Dnj0IeFEhMY1pZtLT04OXXnoJb7/9NiZPnoyqqqpBt9u+fTuefvppdHV1+Ryov3Fmol1W\nqxUlJcdhNOYiJmaC3OEIob3dCofjOJ54Yg4SEhIAAOXl15+nUlUFzlQ0ImC3oO/t7YVer0fYMEdQ\na2tr/8GnJCwm2tTS0oIdO44hImIuoqPFu5JdTh0drbDbj2LRormIj//u/zsWFe3waQDf3NyM3/72\nt/j0009v+ZzRaBy2kABQZCERjVp6qt7yd95tbW0oKTkGo3GOoguJEmYmg4mJSUB4eC5KSo7g6tWr\n/R/Pz4df2l88zsU2YjFZsWIFtmzZgieffBI1NTX9H5ckCYWFhfjjH/+oyFN/iW7U3d2NXbuOICws\nBzExfIMzVjExExASko3du4+gp6dnwOf8VVRITCMWE5vNhh07dmD16tVITk7u/7hOp8P27dshSRKK\niorQ29sb0EC1TJSrY/3NX3m7XC7s21cBp3M6YmMn+uV7BlIgrn73p7i4RPT2TsW+fRWDXpA81qLC\n41xsIxaTjIwMJCcnY+3atTAaBz5uNDw8HM8//zxWrVqF1atXByxIIl989dVxNDbGY+JEk9yhqEZi\n4u2or49FefnxIbfhSkVbRiwma9aswa9+9Su8+eab6OjoGHSb+fPn45tvvvF7cHSdWnqq3vJH3qdO\nnUZlpRMpKbN8DyhIlDozudmkSVk4ccKO6upb7/d1o9EWFR7nYhuxmMTExMBgMGDVqlVISEjAvHnz\nUFxcjL179/YXl7q6Ora5SHGsVivKyuqRkjIXej1vkO1ver0eKSlzUVZWi5aWlhG350pF3UY8NXj5\n8uWoq6tDfn4+qqqq8Pnnn/cfODqdDsnJybBardi4cSN+8YtfBCVof+CpwermcrnwySelcLmyeS1J\ngNlsV2A0fo3Fix9AaOjo7x3LU4rF49N1Jnl5efjqq6/639lJkoTKykrs378f+/fvx759+/CHP/wB\nxcXF/o88gFhM1O3o0UocPgykpmbJHYom1NefQF5eCHJzM73+WhYVcfh0nUliYiIcN6xDdTodZs+e\njVWrVmHXrl04f/48zp07h507d/ovYhpALT1Vb401b6vVikOHrEhOvsu/AQWJKDOTGyUnz8ShQ5ZR\ntbtu1tf+Ki4u1WT7Sy2/3yMWk6VLl+I3v/nNkJ9PT0/H22+/jffff9+vgRGNhcvlwr/+dRKxsbMR\nEjKmx/XQGISGGhAVlYUDB06O+flFd93FmYrIRnU7lc2bN2Pnzp14/vnnkZeX1//x2tpa3H///Vi5\nciU++eQTHD16NKDB+hPbXOpUXX0WZWXXkJqaI3comlRffxQPPRSL6dOn+fy92P5SHp+fZ7Js2TL8\n5S9/QWlp6YBvFB8fD6fTiZdffhmPP/64f6IlGiOHw4Hy8kuYMOFOuUPRrAkTMlBefhFOp9Pn78Wz\nv8Qy6vMlExMTUVxcPOCW8lFRUbh06RKqqqqwZs2agARI6umpesvbvKurz8HhSBX+AVcizkz6GI3j\n0N2djOrqc15/7VD7W+1FRS2/3z6ffG80GjF9Op8HQfLq6elBRUU9Jk68Q+5QNG/ixOk4fNj/156p\nvaiIzqdb0IuMMxN1OXz4JI4fNyIlhS0uJWhqOo05cxyYMydwp2ZzphJ8fn8GPJGSOBwOnDzZhAkT\npsgdCn1rwoQpOHGi0S+zk6FwpaIsLCYCUEtP1Vujzbu2tg4uVzIMBnW8LRV5ZtLHYAiH05mIurr6\nUX/NWI9z0YuKWn6/WUxIeMeO1SA21iR3GHSTmBgTjh+vCdrPE72oiI4zExLalStX8MknZ5CWViB3\nKDSI+voy/PjHMzFhQvDvj8aZiv9xZkKqVV19GeHht8kdBg0hPNyEM2cuy/KzuVIJLhYTAailp+qt\nkfJ2u904d64F48enBCegIFHDzKRPXFwKzp61jurR3oE6zpVeVNTy+81iQsKyWq1wucYjNNQgdyg0\nBIMhDHZ7zJhuAOlvSi8qouPMhIR1+PBJnDwZg6QknhKsZM3NF5Cbey2g15yMBWcq3uPMhFRHkiSc\nOWNBXFyS3KHQCOLiknHmjEXuMG7BlYp/sZgIQC09VW8Nl7fNZsO1a+HC34drMGqamQDX79fV0RGK\n9vb2YbeT6zi/uahMmwa89RZgtwfn56vl95vFhITU1tYGSUqQOwwatQS0tbXJHcSw+orK3/4G7Nhx\nfaUSzKIiOs5MSEhffnkMZ84kYsKENLlDoVGwWmtx111tyM/PljuUUeubqZw69d1MJTxc7qjkxZkJ\nqU5DQzsiI2PlDoNGKTIyFo2NNrnD8ApXKt5hMRGAWnqq3hoqb5fLhba2XkRERAU3oCBR28wEACIi\nomG1dg/7SF+lHueBLipKzdtbLCYknOuD3JgBD2ojZdPr9QCi0dHRIXcoY8aVyvA4MyHhXL58GXv2\ntCM1VVnXLdDwGhpOoLAwHpMnT5Y7FL/Q4kyFMxNSlZ6eXuh0RrnDIK8ZYVfR23iuVAZiMRGAWnqq\n3hoq746OXhgM6i0mapyZAEBoqBGdnUM/ylfU49zXoiJq3jdjMSHhXC8mKu8nqFBYmBEdHf59LryS\naH2lIlQxOXjwIHJzc5GVlYUNGzYMus0LL7yArKws5Ofn4/Tp00GOMDDMZrPcIchiqLy7uuyqXplk\nZprlDiEgDIZwdHUN/ZdVLce5t0VFLXkLU0zcbjeWL1+Ov//97zh69CjeffddVFdXD9hm9+7dOHny\nJCorK7F+/XosW7ZMnmApoK5ds3NlIiCDQd0rk5tpbaUiTDGpqKjAtGnTYDKZYDAYsGTJEuzYsWPA\nNiUlJVi6dCkAIC8vDzabDRaL8m4w5y219FS9NVTeTqcben1IcIMJIrXOTPR6PdzuoZ9rotbjfKSi\nopa8hSkmDQ0NSE9P73+dlpaGhoaGEbepr68PWowUHB6PxGtMhKSDx6Pd0/GHKipOp9yR+Ueo3AGM\n1mj/eNx8DvRwX7ds2TKYTCYAQFxcHLKzs/v7l33vFpTw2mw2KyqeYL7uc+PnJUnCN98chF6v758v\n9L2bV8PrzEyzouLx12u324X09Ou/n97sb7W9zs8HiotLUVUF7NhhxoMPmmEwKCe+G1/3/bumpgYj\nEeaixfLycvz+97/H3r17AQCvv/469Ho9iouL+7dZuXIlzGYzlixZAgDIyMhAWVkZkpJufeYFL1oU\n13vv7UF8/AKEhAjzXogAOJ0OdHYewNKlj8gdiqJ4PIBekB6RKi5anDt3Ls6dO4eamho4HA5s3boV\nRUVFA7YpKirCBx98AOB68YmLixu0kIhGLT1Vbw2Vt16v7jcCap2ZABL0+qE7BVo9zg8eLJU7BL8Q\n5q1daGgo3nvvPSxevBgulwsrVqzAjBkzsGnTJgDAs88+i8LCQhw8eBCZmZkYN24c3n//fZmjpkAI\nCwuF2+3is98F43a7EBam3hMntE6YNpe/sc0lrpKSz3Ht2ixERY2XOxTyQmdnG2Jjq/GDH9wndyg0\nRqpocxH1iYkxwulU6cn6KuZw9CImRr0Xm2odi4kAtNpLHirvmBgjHA71Xvym1pmJ0zl8MeFxLjYW\nExJOVFQ4XC71FhO18njsGDeOdy5QK85MSDi1tbXYs6cNkyaJ8zxxAhobj+EHP0hEWlqa3KHQGHFm\nQqoSExMDQNwn9mlXx7f7jtSIxUQAaumpemvomUkM9PoueDxD3+dJZGqcmbjdboSG9iA6OnrIbXic\ni43FhISj1+sxcWIUuru5OhFFd3c7EhOjeE81FePMhIR0+PBJVFbGITHxNrlDoVGwWC4hN7cLubmZ\ncodCPuDMhFQnMTEOTqdN7jBolNxuGyZMiJU7DAogFhMBqKWn6q3h8k5ISABgDVoswaTGmYkktXy7\nz4bG41xsLCYkpKioKMTH6zk3EUBXlw2JiQaMGzdO7lAogDgzIWGdPHkKhw4ZkJIyXe5QaBiNjadx\n330SZs2aIXco5CPOTEiV0tKS4fGI/1hmtZMkC1JTk+UOgwKMxUQAaumpemukvOPj4xEZ2a26+3Sp\naWZit/cgKsqOuLi4EbflcS42FhMSlk6nw4wZyWhrq5c7FBrC1av1mDEjmdeXaABnJiS09vZ2fPjh\nEaSmzucfLIWRJAmNjfvx1FPzeBsVleDMhFQrNjYW6elhaG9X52nCIrPZLDCZIlhINILFRABq6al6\na7R5z55tQldXTUBjCSa1zEyuXatBZqZp1NvzOBcbiwkJb9KkSYiKssFu75Y7FPpWT08XoqM7kJKS\nIncoFCScmZAqnDp1Gp9/bkdq6my5QyEA9fXHYTaPw4wZvAZITTgzIdWbPn0qIiMt6OnpkjsUzevu\n7kBUlBV33HG73KFQELGYCEAtPVVveZO3wWDAPfdMRUvL6cAFFCSiz0xaW0/jvvvuQGhoqFdfx+Nc\nbCwmpBq3325CfLwNXV28m7BcOjvbMGFCJ0wmPhpAazgzIVWpra3Fzp31SE+/V+5QNKmu7gsUFd3G\n57yrFGcmpBnp6emYMsWNlpZauUPRHKu1BlOnAqmpqXKHQjJgMRGAWnqq3hpL3jqdDvffnw23+zTs\n9h7/BxUEIs5M7PZueDxncN99s8d8JwIe52JjMSHViY6Oxve+dzuuXDkpdyiaYbGcwAMPTENUVJTc\noZBMODMhVZIkCXv3/hvNzbdhwoTJcoejalZrDSZNqscjj9zH+6OpHGcmpDl97S6P5zS6uzvlDke1\nurs7oNOdxf33Z7OQaByLiQDU0lP1lq95R0dH49FHZ6K1tQJOp8M/QQWBKDMTp9OO1tbDeOSRWX5p\nb/E4FxuLCalaamoqCgomoanpKNuafuTxeNDYeARmcxomTZokdzikAJyZkOpJkoTPPz+MqqoIpKZm\nyh2OKjQ0VGLWLDvuv/9uuUOhIOLMhDRNp9Ph3ntzkZLSiitXLskdjvCuXLmA1NSryM/PkTsUUhAW\nEwGopafqLX/mHRoaiocfzkNMzEW0tFz22/cNBCXPTK5cuYS4uMtYsCDP63tvjYTHudhYTEgzIiIi\n8IMf3IPIyHO8Qn4MrNYaREdfRGHhPTAajXKHQwrDmQlpzrVr17BrVzm6uqZi4kST3OEI4cqVi4iJ\nufRtMY6UOxySyXB/N1lMSJN6enqwe/dXaG2dhJSUDLnDUSxJktDUdBoTJzbhscfuQUREhNwhkYw4\ngBecWnqq3gpk3hEREVi48H5MmdKGurrDcLtdAftZ3lLKzMTtdqG+/jCmTbPh8cfvD3gh4XEuNhYT\n0qywsDA89FA+8vLCUV//bz5D/ga9vdfQ0PBv3HNPBObPz0dYWJjcIZHCsc1FBODSpRr83/+dRVRU\nDmJjJ8odjqxstivo6TmBhx++kw+5ogE4MxkEiwndrLW1Ffv3n0Bb20QkJ9+FkBD/nvqqdC6XExZL\nFRISWvDQQzmIj4+XOyRSGM5MBKeWnqq3gp13QkICfvjDB3D33UBTUxna261B/fl95JiZ2GxX0Nxc\nhnnz9Fi8+AFZCgmPc7Fp660X0QhCQ0MxZ04WJk+24sCBSjQ0TERiYgYMBnXODJxOO65cOY2EhBY8\n9lg2JkyYIHdIJCi2uYiG4HK5cOrUGVRU1AOYgokTb1dN68vtdsFqvQCdrgZ3352GmTPv9PsV7aQ+\nnJkMgsWERqu7uxuVlWdw4oQVoaF3YOLE26DXi9kh9ng8aGmpgct1Hjk5icjMvJPXjtCocWYiOLX0\nVL2llLwjIyORn5+Dp57Kx7RpVjQ1/QtNTefgdNoD8vMCMTNxOu1obj6Lpqb9uOOOVjz11D2YNy9b\nUYVEKfs72NSSN9e1RKMUExMDs3kecnLacf78ZZw8eQAORyJiYkyIjlbmmU8dHa3o6KhBeLgVc+ak\nYtq0PMTExMgdFqkQ21xEY+R0OlFfX4/jx2tgsQAhISmIjk7CuHFxsj3CVpIkdHVdRWenBZLUjKQk\nHXJyTEhLS+NMhHzGmckgWEzIn2w2GxoamnH2rAUWix1AEiIjrxeWsLDA3mHX4ehFV9dV9PRYAFxB\ncnI47rwzGSkpSYiLiwvozyZtYTEZhEjFpLS0FGazWe4wgk7UvLu7u2GxWHDx4hU0NNjQ06MHEAsg\nFpGRcQgLi0RYmBGhoYZBv/7rr0uRmWm+5eMulxMORy/s9mvo6WmHTtcOj8eGceOASZNiMXVqEpKS\nkhQ1B/GGqPvbVyLlPdzfTa57ifwsMjISU6ZMwZQpUwBcv0Nxe3s72tpsaGysgc3Wg/b2XtjtHgDh\n0OmMAMIgSToAerS2nkVjYxR0OgmAA5LUC8CO8HA9oqKMSE2NQEpKLOLjJyM2NlPY4kHqwpUJkUzc\nbjfsdjt6e3vhcDggSRI8Hg90Ol3/f2FhYTAajQgPD0dISIjcIZPGsc01CBYTIiLvqOI6k87OTixa\ntAhZWVlYvHgxurq6Bt3OZDIhKysLOTk5mDdvXpCjDAy1nIfuLeatLcxbbMIUk1deeQX33nsvKisr\nkZ+fj1c/WXjgAAAIm0lEQVRffXXQ7XQ6HUpLS3H8+HFUVFQEOUoiIm0Sps2VkZGBsrIyJCUlobm5\nGWazGadPn75luylTpuDIkSNISEgY9vuxzUVE5B1VtLksFguSkpIAAElJSbBYLINup9PpMH/+fOTk\n5OCdd94JZohERJqlqFODFyxYgObm5ls+/tprrw143Xemy2C++OILpKSkoLq6GoWFhcjIyEBBQcGg\n2y5btgwmkwkAEBcXh+zs7P7zvfv6mEp4fWNPVQnxBOv1iRMn8NxzzykmnmC95v5WRjzBeq3k/d33\n75qaGoxEqDZXaWkpkpOT0dTUhAcffHDQNteNVq1ahdTUVPz617++5XMitblKBbqoyZ+Yt7Ywb+VT\nxanBq1evRkJCAoqLi7Fu3TrYbDasW7duwDbd3d1wu92Ijo6G1WpFQUEB/vznP+Phhx++5fuJVEyI\niJRAFcWks7MTTz/9NC5evIipU6fir3/9K6KiotDY2IgVK1Zg165duHjxIn74wx8CuP4I1p/85Cd4\n9tlnB/1+LCZERN5RRTHxN5GKiUjLYH9i3trCvJVPFWdzERGRcnFlQkREo8KVCRERBRSLiQBuPOdb\nS5i3tjBvsbGYCODEiRNyhyAL5q0tzFtsLCYCsNlscocgC+atLcxbbCwmRETkMxYTAYzmvjhqxLy1\nhXmLTdOnBhMRkXeGKhmKumtwMGm0hhIRBQTbXERE5DMWEyIi8hmLicJ8/PHHmDlzJkJCQnDs2LEh\ntzt48CByc3ORlZWFDRs2BDHCwOjs7MSiRYuQlZWFxYsXo6ura9DtTCYTsrKykJOTg3nz5gU5Sv8Z\nzf574YUXkJWVhfz8/BGf3SOKkfIuLS1FbGwscnJykJOTg1dffVWGKP1r+fLlSEpKQmZm5pDbqGJf\nS6Qo1dXV0pkzZySz2SwdPXp00G1cLpc0depU6dKlS5LD4ZBmz54tVVVVBTlS/3r++eelN954Q5Ik\nSVq3bp1UXFw86HYmk0lqbW0NZmh+N5r9t2vXLumxxx6TJEmSysvLpby8PDlC9avR5H3gwAFp4cKF\nMkUYGAcPHpSOHTsmzZo1a9DPq2Vfc2WiMBkZGZg+ffqw21RUVGDatGkwmUwwGAxYsmQJduzYEaQI\nA6OkpARLly4FACxduhSffvrpkNtKgp88MZr9d+P/H3l5ebDZbLBYLHKE6zejPW5F3783KygowPjx\n44f8vFr2NYuJgBoaGpCent7/Oi0tDQ0NDTJG5DuLxYKkpCQAQFJS0pC/TDqdDvPnz0dOTg7eeeed\nYIboN6PZf4NtU19fH7QYA2E0eet0Onz55ZeYOXMmCgsLUVVVFewwg04t+1qzpwbLacGCBWhubr7l\n42vXrsXChQtH/HpRr5EZKu/XXnttwGudTjdkjl988QVSUlJQXV2NwsJCZGRkoKCgICDxBspo99/N\n79BF3e99RhN/bm4u6urqYDAYsGXLFhQVFeH8+fNBiE5eatjXLCYy+Oc//+nT16empqKurq7/dV1d\nHdLS0nwNK+CGyzspKQnNzc1ITk5GU1MTEhMTB90uJSUFADBjxgwsXrwYFRUVwhWT0ey/m7epr69H\nampq0GIMhNHkHR0d3f/vn//85yguLkZbWxvi4+ODFmewqWVfs82lYEP1jufOnYtz586hpqYGDocD\nW7duRVFRUZCj86+ioiJs2bIFALBlyxYsWrTolm26u7vR2dkJALBardi9e/ewZ8go1Wj2X1FRET74\n4AMAQHl5OeLi4vrbgKIaTd4Wi6X/uP/HP/6BiIgIVRcSQEX7Ws7pP93q73//u5SWliYZjUYpKSlJ\nevTRRyVJkqSGhgapsLCwf7vS0lIpOztbmjVrlrR+/Xq5wvWbjo4O6YknnpAyMzOlRYsWSZ2dnZIk\nDcz7woUL0uzZs6XZs2dL8+fPl9566y05Q/bJYPvvrbfeGpBTcXGxNGvWLCkvL0/4s/X6jJT3xo0b\npZkzZ0qzZ8+Wnn76aenIkSNyhusXS5YskVJSUiSDwSClpaVJ7777rir3tWbvzUVERP7DNhcREfmM\nxYSIiHzGYkJERD5jMSEiIp+xmBARkc9YTIiIyGcsJkRBtm3bNixcuBB5eXlYvnw5urq60NjYiBUr\nVqCgoAALFy4c9LYzRErG60yIgmjz5s04ceIE1q1bh7CwMDzyyCNISEiA0+nExo0b0dzcjCeeeAJL\nly7FK6+8Ine4RKPGlQlRkFy5cgU7d+7Em2++CaPRCL1ej/z8fPztb39DQUEBUlJS8PHHH6OxsRFz\n5syRO1wir3BlQhQkb7zxBsxmM/Ly8vo/9qMf/QglJSVobW1FTEwMent70dzcDJPJJF+gRGPAYkIU\nJJIkDbi1uCRJSEpKgslkQkVFhYyREfmObS6iILn5GRWVlZVoaWnBgw8+OOzXud1u7N69G1OmTAlk\neEQ+4fNMiGSyf/9+ABi2mHz44Yc4dOgQWlpacPny5WCFRuQ1rkyIgqSrqwu1tbX9r/fv3w+dTod7\n7723/2N2ux2FhYX9r5966imsX78eCxYsCGqsRN7iyoQoCHp7e3HHHXfAYDCgtrYWdXV12LdvH5KT\nkxETE9O/3Z///GfMnz9fxkiJxoYrE6IgqK2thdVqxcqVK+HxeLBmzRr84Q9/gNVq7X9k67/+9S+U\nlZVh1apVMkdL5D2uTIiCYPr06XjppZfwySef4OjRo3jmmWfw+OOPY/z48fjpT38Ko9GIWbNmYdu2\nbdDr+R6PxMNTg4kEsHnzZixfvhwej0fuUIgGxbdARETkMxYTIiLyGYsJERH5jMWESABOpxMA4HK5\nZI6EaHAsJkQKtnPnTjz++ON4/vnnodPpkJWVhSeffBIdHR1yh0Y0AM/mIiIin3FlQkREPmMxISIi\nn7GYEBGRz1hMiIjIZywmRETkMxYTIiLyGYsJERH5jMWEiIh8xmJCREQ+YzEhIiKf/T/vDVaJTOlb\nngAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that the line is the constraint so that any solution to this problem must be on this line (i.e. satisfy the constraint). The $L_2$ solution is the one that just touches the perimeter of the circle. This is because, in $L_2$, the unit-ball has the shape of a circle and represents all solutions of a fixed $L_2$ length. Thus, the one of smallest length that intersects the line is the one that satisfies the stated minimization problem. Intuitively, this means that we *inflate* a ball at the origin and stop when it touches the contraint. The point of contact is our $L_2$ minimization solution.\n",
      "\n",
      "Now, let's do same problem in $L_1$ norm\n",
      "\n",
      "min $||\\mathbf{x}||_1=|x_1|+|x_2|$\n",
      "\n",
      "where $x_1 + 2 x_2 = 1$\n",
      "\n",
      "In this case the constant-norm unit-ball contour in the $L_1$ norm is a diamond-shape instead of a circle. Comparing the  graph below to the last shows that the solutions found are different. Geometrically, this is because the line tilts over in such a way that the inflating circular $L_2$ ball hits a point of tangency that is different from the $L_1$ ball because the $L_1$ ball creeps out mainly along the principal axes and is less influenced by the tilt of the line. This effect is much more pronounced in higher $N$-dimensional spaces where $L_1$-balls get more *spikey*. \n",
      "\n",
      "The fact that the $L_1$ problem is less sensitive to the tilt of the line is crucial since that tilt (i.e. orientation) is random due the choice of random measurement matrices. So, for this problem to be well-posed, we need to *not* be influenced by the orientation of any particular choice of random matrix and this is what casting this as a $L_1$ minimization provides. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib.patches import Rectangle\n",
      "import matplotlib.patches\n",
      "import matplotlib.transforms\n",
      "\n",
      "r=matplotlib.patches.RegularPolygon((0,0),4,1/2,pi/2,alpha=0.5)\n",
      "\n",
      "fig=figure()\n",
      "ax=fig.add_subplot(111)\n",
      "ax.plot(x1,(1-x1)/2)\n",
      "ax.plot(0,1/2,'rs')\n",
      "ax.add_patch(r)\n",
      "ax.grid()\n",
      "ax.set_xlabel('$x_1$',fontsize=24)\n",
      "ax.set_ylabel('$x_2$',fontsize=24)\n",
      "ax.axis('equal')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "(-1.0, 1.0, -0.60000000000000009, 1.2)"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEVCAYAAAAl9QikAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHqBJREFUeJzt3X1QVOfZBvBriRgSP4IwAzKCbiJxoCqKtS6ThmZ1apJS\nIZg6KW3GgnSIttPM69hRYjNJ42fNO52+frVRO4li23T8SBQi6lRNFtokhFFjTIJNMMkqYqBKpOKk\nfKjn/cPuCssu7Nlz9pzzPOf6zWTCYQ/kuXN2ubife3dxKIqigIiISIMYsxdARETiY5gQEZFmDBMi\nItKMYUJERJoxTIiISDOGCRERaTbE7AWYxeFwmL0EIiLhhHo1ia07E0VRhPjn17/+telrYN2sm3Wz\n7oHYOkxE4fV6zV6CKVi3vbBusTFMiIhIM4aJAEpKSsxegilYt72wbrE5lME2wiTlcDgG3QMkIqLb\nBvq5yc5EAB6Px+wlmIJ12wvrFhvDhIiINOM2FxERhYXbXEREFFUMEwHIsqeqFuu2F9YtNoYJERFp\nxpkJERGFhTMTIiKKKoaJAGTZU1WLddsL6xYbw4SIiDTjzISIiMLCmQkREUUVw0QAsuypqsW67YV1\ni41hQkREmgk1MyktLUV1dTWSkpLw4YcfBj1n+fLlqK6uxt13340dO3YgIyMj6HmcmRARqSPNzGTB\nggU4fPhwyNsPHjyIDz74AKdPn8aGDRuk+aMzRERWJ1SY5ObmYtSoUSFvr6qqQnFxMQDA5XKhvb0d\nra2tRi0vamTZU1WLddsL6xabUGEymObmZqSlpfmPU1NTceHCBRNXpN1XXwE9PWavgohoYEPMXoDe\nAvfzHA5HyHNLSkrgdDoBAPHx8Zg6dSrcbjeA278tmH188qQbGza4MW+eB48+Csyeba31RfvYxyrr\nMeLY7XZbaj1GHvtYZT12v96+j71eLwYj1AAeALxeL/Lz84MO4BctWgS3242ioiIAQEZGBmpqapCc\nnNzvXJEG8HV1wIoVQEMD8KtfAQsWAEOHmr0qIrIbaQbwgykoKMDOnTsBAHV1dYiPjw8aJKLp7PTg\n0CFg1y5g/37g/vuBrVuB7m6zVxZdgb+t2gXrthdZ6hZqm+tHP/oRampqcPnyZaSlpWHFihXo+e9A\nYeHChcjLy0NtbS0mT56MYcOGYfv27SavWF85OcChQ7c7lbVr2akQkTUIt82lF5G2uULh9hcRGck2\n21x24+tU7Lb9RUTWwzARwGB7qrKGiix7yWqxbnuRpW6GiURkDRUisj7OTCTGmQoR6YkzE5tip0JE\nRmGYCEDrnqqooSLLXrJarNteZKmbYWIjooYKEVkfZyY2xpkKEanBmQkFFdippKcDW7YAXV1mr4yI\nRMMwEUC091R9obJnD1BZeWv7ywqhIsteslqs215kqZthQn4u1+1QqaqyTqgQkfVxZkIhvfferZnK\nRx/dnqnceafZqyIis3BmQhFxuYCDB9mpENHgGCYCMHtP1axQMbtus7Bue5GlboYJhY2dChGFwpkJ\nRYwzFSJ74cyEooKdChH5MEwEYPU91WiFitXrjhbWbS+y1M0wId2wUyGyL85MKGo4UyGSC2cmZAp2\nKkT2wTARgOh7qpGGiuh1R4p124ssdTNMyDDsVIjkxZkJmYYzFSKxcGZClsROhUgeDBMByLKnGkqo\nUDlyxGP20kwh+/UOhXWLjWFClhEYKk8+yU6FSBScmZBlcaZCZC2cmZCQOFMhEgfDRACy7Kmq5avb\nbqFi9+ttN7LUzTAhYdgtVIhEwpkJCYszFSJjcWZCUmKnQmQdDBMByLKnqla4dcsWKrze9iJL3QwT\nkoZsoUIkEs5MSFqcqRDpizMTsiV2KkTGYZgIQJY9VbX0qlu0UOH1thdZ6maYkG2IFipEIuHMhGyL\nMxUidTgzIQqCnQqRfhgmApBlT1Uto+q2WqjwetuLLHUzTIj+y2qhQiQSzkyIQuBMhagvzkyIIsBO\nhSh8DBMByLKnqpZV6jY6VKxSt9FYt9gYJkRhYqdCFBpnJkQR4kyF7IYzE6IoYKdCdBvDRACy7Kmq\nJUrdeoeKKHXrjXWLjWFCpJPeoVJZyU6F7IUzE6Io4UyFZMOZCZEJOFMhOxEuTGprazFt2jRkZWVh\n06ZN/W73eDy45557kJ2djezsbKxevdqEVepLlj1VtWSpW22oyFK3WqxbbEPMXoAaN27cQGlpKY4e\nPYoxY8bgW9/6Fr773e8iMzOzz3kPPfQQqqqqTFolUXC+UPFtf61de3v76/+efgqdn34KAPC2t8MT\nHw8AiJswAc9s22bmsonCIlSY1NfXIz09HU6nEwBQVFSEysrKfmEi2yzE7XabvQRTyFp3sFB5NO5T\nbGus6XfuC8YvzzSyXu/ByFK3UNtczc3NSEtL8x+npqaiubm5zzkOhwPvvPMOJk6ciLy8PDQ0NBi9\nTKKw9N7+amszezVE2gjVmTgcjkHPmTZtGpqamhAbG4uKigoUFBTg7NmzQc8tKSnxdznx8fGYOnWq\n/7cE3z6mFY5776laYT1GHZ86dQqLFy+2zHqidexyASPS2uH5CnADuHVLX1Zab7SO7XK9A4+t/Pj2\nfez1ejEoRSDvvvuu8sgjj/iP165dq6xbty7k+Tdv3lQSEhKUtra2freJVPpbb71l9hJMYae6n/vO\ndxQFUBRAeeu//1YApez+h5TOTrNXZww7Xe/eRKp7oJ+bQm1zTZ8+HY2NjfB6veju7sauXbtQUFDQ\n55zW1lb/zOSNN97AXXfdhYSEBDOWqxvfbwt2Y5e6u7q6cLG51X/s7nXb5cv2eUqxXa53IFnqFmqb\na8iQIXjllVcwd+5cXL9+HWVlZcjMzMTWrVsBAAsXLsTevXvx0ksvYciQIcjKykJlZaXJqyYKraur\nC7///Z9x4Q4nSsaN7nNbZ2cb7rpPwav/24O1a2P7PPuLL34kq+Er4AXg8Xik+e1FDdnr9gXJhx8m\nY+zY7/tngl6vB06nG4qi4Ny5N5Cd3Yaf/exJvP/+UKlfUS/79Q5FpLr5CngiiwkVJL05HA6MG5eP\n999PxEsv/QXZ2d18RT1ZFjsTIoOFEyS9BXYoQ4cOBcD3/iLjsTMhsgi1QQL071C6u7sB8F2KyVoY\nJgLo/ZxvO5Gt7nCDxOv19PtcqEABboXKoUPA7t1ih4ps1ztcstTNMCEyQCQdSaCBAgUAcnLkCBUS\nE2cmRFGmR5D0FmqGEqiu7tZM5eOPOVMhfXBmQmQSvYMEGLxD8WGnQkZimAhAlj1VtUSvO9IgCTYz\nCRRuoADihIro1ztSstTNMCGKgmh0JIHUBAogTqiQmDgzIdKZEUHSW7gzlEC+mUpDw+2ZSphfSjbF\nmQmRQYwOEkB9h+Lj61R27QL277/VqWzdCoT55UR9RBwm//73v/HXv/4VS5YswbZt23DlyhX/bZ99\n9hlWr16NV199FcePH9dloXYmy56qWqLVrVeQhDMzCRRpoADWCRXRrrdeZKk7ojC5cuUKHnzwQTz5\n5JNYv349Fi1ahMzMTBw+fBgAMH78eCxbtgz/+Mc/4HK5dF0wkRWZ0ZEE0hIogHVChcQUUZgsW7YM\n6enpOHXqFDo6OvDxxx+jrKwMP/nJT7B//34AwNChQzF27FjOJXQgyjuK6k2UuvUOEqfTHfHXag0U\nwLxQEeV6602WuiMKk08//RT79u1DVlYWhg0bhszMTKxatQrHjx/Hpk2b/B0Kkeys0JEE0iNQAHYq\npE5EYTJu3Lignx87diwOHTqEAwcO4NChQ5Z4YMlAlj1Vtaxed7SCJJKZSSC9AgUwLlSsfr2jRZa6\nIwqTtLQ0NDY24ne/+x3WrFnT57ahQ4di8+bN+OSTT/C3v/2NgUJSsmJHEkjPQAHYqdDAInqdyX/+\n8x88++yz2LZtG8aOHYuGhoag5+3btw/z58/HtWvXNC9Ub3ydCUVKhCDpLdLXoQyGr1Oxn4F+bmp6\n0WJnZydiYmIGvHO2tbUhMTEx0v9E1DBMKBKiBYlPtAIFYKjYiaYXLba0tOC5557zP0urt7i4uEHv\nlFYMEtHIsqeqltXqNipI9JiZBNJ7y6s3vba/rHa9jSJL3YOGSVlZGSoqKjBv3jx4vV7/5xVFQV5e\nHn7729/i5s2b0VwjkelE7Uh6i2agAJyp2N2g21y5ubnYuHEj9uzZg+effx5xcXH+27q6urBx40bU\n1NRg7969fW6zOm5zUbhkCJLeornl1Ru3v+SjaWZSVlaGlStXIiUlJeQ5b775Jvbv34+NGzdqW6mB\nGCYUDtmCxMeoQAEYKjLRNDNZsWIFfvGLX2D9+vW4evVq0HNmzZqFjz76SNsqKSRZ9lTVMrtus4Ik\nGjOTQNHe8uot3O0vs6+3WWSpe9AwGTlyJGJjY7FkyRIkJiZixowZKC8vx+HDh/3h0tTUhM7Ozqgv\nlsgosnYkvRkZKABnKrIbdJurtLQUTU1NyMnJQUNDA/7+97/j8uXLt77Y4cDo0aNx6dIlbN68GU89\n9ZQhi9YDt7koFDsESW9Gbnn1xu0v8WiambhcLrz77ruIibnVxCiKgtOnT+PYsWM4duwYjh49ipUr\nV6K8vFz/lUcRw4SCsVuQ+JgVKABDRSSaZiZJSUl92l+Hw4EpU6ZgyZIlqK6uxtmzZ9HY2IgDBw7o\nt2LqQ5Y9VbWMrtsqQWLEzCSQ0Vtevfm2v8rLPbbc/pLl8T1omBQXF+OZZ54JeXtaWhq2bduG7du3\n67owIiNZJUjMZGagAMA3vsGZisjCejuVHTt24MCBA1i6dGmfP3Z1/vx5PPjgg1i0aBFee+01nDhx\nIqqL1RO3uciHQdKXmVtevXH7y3o0/w34kpIS/OEPf4DH4+nzjRISEtDT04Pnn38ec+bM0We1RAZi\nkPRndofiw2d/iSXst6BPSkpCeXl5nwfb8OHD8cUXX6ChoQErVqyIygJJnj1VtaJdt1WDxIyZSSAz\nAiXU9ZY9VGR5fEf090x6i4uLw4QJE/RYC5FhrBokVmKVDsVH9lARnaa3oBcZZyb2xSBRxyozlECc\nqRhP88yESBYMEvWs1qH4sFOxFoaJAGTZU1VL77pFCRIrzEwCGREokV5v0UNFlsc3w4RsQZQgsTKr\ndig+ooeK6DgzIekxSPRl1RlKIM5U9MeZCdkWg0R/Vu9QfNipGIthIgBZ9lTV0lq3qEFixZlJoGgE\nSrTu51YPFVke3wwTkpKoQSISUToUH6uHiug4MyHpMEiMJcoMJRBnKupxZkK2wSAxnmgdig87FX0x\nTAQgy56qWmrrliVIRJiZBNIjUMy6nweGSno6sGUL0NVlzH9flsc3w4SkIEuQiEzUDsXHFyq7dwOV\nlbc6FSNDRXScmZDwGCTWIuoMJZBvpvLxx7dnKnfeafaqzMWZCUmLQWI9oncoPuxU1GGYCECWPVW1\nBqtb1iARcWYSKJJAser9PNqhYtW61WKYkJBkDRKZyNKh+LBTGRhnJiQcBolYZJmhBLLjTIUzE5IG\ng0Q8snUoPuxU+mKYCECWPVW1Auu2S5DIMDMJFE6giHo/1xoqotYdiGFCQrBLkMhM1g7Fx+6dilAz\nk9raWixevBjXr19HWVkZnn766X7nLF++HNXV1bj77ruxY8cOZGRkBP1enJmIg0EiF1lnKIFknKlI\nMTO5ceMGSktL8frrr+PEiRN4+eWXcebMmT7nHDx4EB988AFOnz6NDRs2oKSkxJzFkm4YJPKRvUPx\nsVunIkyY1NfXIz09HU6nE7GxsSgqKkJlZWWfc6qqqlBcXAwAcLlcaG9vR2trqxnL1ZUse6pqHTly\nxJZBIuPMJFCwQJH1fj5YqMhStzBh0tzcjLS0NP9xamoqmpubBz3nwoULhq2R9NPV1YV9+47YLkjs\nJDBQenp6zF5SVIUKFVnKHmL2AsIV7g+TwP28gb6upKQETqcTABAfH4+pU6fC7XYDuP3bghWO3W63\npdZjxHFFRQXOnLmKb397NhwOh/+3dafz1u0yHzudbkutJ9rHqakPY//+/8Hw4TPgY/b9L5rHOTlA\nebkHDQ1AZaUbM2e6ERtrnfX1PvZ97PV6MRhhBvB1dXV44YUXcPjwYQDAb37zG8TExKC8vNx/zqJF\ni+B2u1FUVAQAyMjIQE1NDZKTk/t9Pw7grU1RFOzeXYXq6q8wbtyTuOMOOYe0dnf9ehfOnfszCguT\n8fjj9uxAb94EYgTZI5JiAD99+nQ0NjbC6/Wiu7sbu3btQkFBQZ9zCgoKsHPnTgC3wic+Pj5okIhG\nlj1VNRwOB5KSRuL730/AuXN/wY0bcg5pg7HDzAToHyQ1NTVmL8kUtbUes5egC2G2uYYMGYJXXnkF\nc+fO9T81ODMzE1u3bgUALFy4EHl5eaitrcXkyZMxbNgwbN++3eRVkxYOhwNPPFEAoArV1X9hhyIR\ndiTyEWabS2/c5hIHt7zkwiARlxTbXGRfvg7FjltesmGQyIthIgA7zkyAvnXbKVBknZkMFiS8n4uN\nYULCsFOgyIYdifw4MyHhcIYiFgaJPDgzIamwQxEHg8Q+GCYCkGVPVa2B6pY5UGSZmagNEt7PxcYw\nIWHJHCiiY0diP5yZkPA4Q7EWBom8ODMhqbFDsQ4GiX0xTAQgy56qWmrqlilQRJ2ZaA0S3s/FxjAh\nacgUKKJhR0KcmZB0OEMxFoPEPjgzIVthh2IcBgn5MEwEIMueqlpa6hY5UESZmegdJLyfi41hQtIS\nOVCsjh0JBeLMhKTHGYq+GCT2xZkJ2Ro7FP0wSCgUhokAZNlTVUvPukUKFKvOTKIdJLyfi41hQrYh\nUqBYDTsSGgxnJmQ7nKGowyAhH85MiHphhxI+BgmFi2EiAFn2VNWKZt1WDhSrzEyMDhLez8XGMCHb\nsnKgmI0dCanFmQnZHmcofTFIKBTOTIgGwA7lNgYJRYphIgBZ9lTVMrJuKwWKWTMTs4OE93OxMUyI\n/stKgWI0s4OExMeZCVEAu81QGCQULs5MiFSwU4fCICG9MEwEIMueqlpm1m1moBg1M7FakPB+LjaG\nCVEIMncoVgsSEh9nJkSDkG2GwiChSHFmQqSBTB0Kg4SihWEiAFn2VNWyUt1GBkq0ZiZWDxIrXW8j\nyVI3w4QoTCJ3KFYPEhIfZyZEKok2Q2GQkF44MyHSkUgdCoOEjMIwEYAse6pqWbnuaAaKXjMT0YLE\nytc7mmSpm2FCFCErdyiiBQmJjzMTIo2sNkNhkFC0cGZCFEVW6lAYJGQWhokAZNlTVUukuvUMlEhn\nJqIHiUjXW0+y1M0wIdKJmR2K6EFC4uPMhEhnRs9QGCRkFM5MiAxkZIfCICGrYJgIQJY9VbVErltL\noIQ7M5EtSES+3lrIUjfDhChKotmhyBYkJD7OTIiiTO8ZCoOEzMKZCZGJ9OxQGCRkVQwTAciyp6qW\nTHWrCZRQMxPZg0Sm662GLHUzTIgMoqVDkT1ISHycmRAZTO0MhUFCVsGZCZGFqOlQGCQkCmHCpKOj\nA4WFhcjKysLcuXNx7dq1oOc5nU5kZWUhOzsbM2bMMHiV0SHLnqpaMtc9UKD4ZiZ2CxKZr/dAZKlb\nmDBZtWoVHnjgAZw+fRo5OTlYvXp10PMcDgc8Hg/ef/991NfXG7xKovANFCh2CxISnzAzk4yMDNTU\n1CA5ORktLS1wu9345z//2e+8e++9F8ePH0diYuKA348zE7KKwBmKoigMErKkgX5uChMmo0aNwpUr\nVwDcevAlJCT4j3u77777MGLECMTExODnP/85ysrKgn4/hglZSe9AUZSbDBKypIF+bg4xeC0Dmj17\nNlpaWvp9fs2aNX2OHQ5HyAfZ22+/jZSUFJw5cwZ5eXnIyMhAbm5u0HNLSkrgdDoBAPHx8Zg6dSrc\nbjeA2/uYVjjuvadqhfUYdXzq1CksXrzYMuuJ5nFNTQ2SkkaisDAWDQ0fIiHhPtTU1Fhmfbze0Tu2\n8uPb97HX68VghOlMMjIy4PF4MHr0aHz55ZeYOXNm0G2u3pYsWYIxY8bgl7/8Zb/bROpMPB6P/yLb\nCeu2F9ZtfVJscy1btgyJiYkoLy/HunXr0N7ejnXr1vU55+uvv8aNGzcwYsQIXLp0Cbm5udi4cSMe\nfvjhft9PpDAhIrICKcKko6MD8+fPx+eff47x48fjT3/6E4YPH46LFy+irKwM1dXV+Pzzz/H4448D\nABITE/HEE09g4cKFQb8fw4SISB0pwkRvIoWJSG2wnli3vbBu6+Mr4ImIKKrYmRARUVjYmRARUVQx\nTATQ+znfdsK67YV1i41hIoBTp06ZvQRTsG57Yd1iY5gIoL293ewlmIJ12wvrFhvDhIiINGOYCCCc\n98WREeu2F9YtNls/NZiIiNQR4l2DjWTTDCUiigpucxERkWYMEyIi0oxhYjF79uzBxIkTcccdd+Dk\nyZMhz6utrcW0adOQlZWFTZs2GbjC6Ojo6EBhYSGysrIwd+5cXLt2Leh5TqcTWVlZyM7OxowZMwxe\npX7CuX7Lly9HVlYWcnJyBv3bPaIYrG6Px4N77rkH2dnZyM7OxurVq01Ypb5KS0uRnJyMyZMnhzxH\nimutkKWcOXNG+eSTTxS3262cOHEi6DnXr19Xxo8fr3zxxRdKd3e3MmXKFKWhocHglepr6dKlyosv\nvqgoiqKsW7dOKS8vD3qe0+lU2trajFya7sK5ftXV1cr3vvc9RVEUpa6uTnG5XGYsVVfh1P3WW28p\n+fn5Jq0wOmpra5WTJ08qkyZNCnq7LNeanYnFZGRkYMKECQOeU19fj/T0dDidTsTGxqKoqAiVlZUG\nrTA6qqqqUFxcDAAoLi7G/v37Q56rCP7kiXCuX+//Hy6XC+3t7WhtbTVjuboJ934r+vUNlJubi1Gj\nRoW8XZZrzTARUHNzM9LS0vzHqampaG5uNnFF2rW2tiI5ORkAkJycHPLB5HA4MGvWLGRnZ+OPf/yj\nkUvUTTjXL9g5Fy5cMGyN0RBO3Q6HA++88w4mTpyIvLw8NDQ0GL1Mw8lyrW371GAzzZ49Gy0tLf0+\nv3btWuTn5w/69aK+RiZU3WvWrOlz7HA4Qtb49ttvIyUlBWfOnEFeXh4yMjKQm5sblfVGS7jXL/A3\ndFGvu0846582bRqampoQGxuLiooKFBQU4OzZswaszlwyXGuGiQmOHDmi6evHjBmDpqYm/3FTUxNS\nU1O1LivqBqo7OTkZLS0tGD16NL788kskJSUFPS8lJQUAkJmZiblz56K+vl64MAnn+gWec+HCBYwZ\nM8awNUZDOHWPGDHC//FPf/pTlJeX46uvvkJCQoJh6zSaLNea21wWFmrvePr06WhsbITX60V3dzd2\n7dqFgoICg1enr4KCAlRUVAAAKioqUFhY2O+cr7/+Gh0dHQCAS5cu4eDBgwM+Q8aqwrl+BQUF2Llz\nJwCgrq4O8fHx/m1AUYVTd2trq/9+/8Ybb+Cuu+6SOkgAia61mdN/6u/1119XUlNTlbi4OCU5OVl5\n9NFHFUVRlObmZiUvL89/nsfjUaZOnapMmjRJ2bBhg1nL1c3Vq1eVxx57TJk8ebJSWFiodHR0KIrS\nt+7PPvtMmTJlijJlyhRl1qxZypYtW8xcsibBrt+WLVv61FReXq5MmjRJcblcwj9bz2ewujdv3qxM\nnDhRmTJlijJ//nzl+PHjZi5XF0VFRUpKSooSGxurpKamKi+//LKU19q2781FRET64TYXERFpxjAh\nIiLNGCZERKQZw4SIiDRjmBARkWYMEyIi0oxhQmSwvXv3Ij8/Hy6XC6Wlpbh27RouXryIsrIy5Obm\nIj8/P+jbzhBZGV9nQmSgHTt24NSpU1i3bh2GDh2KRx55BImJiejp6cHmzZvR0tKCxx57DMXFxVi1\napXZyyUKGzsTIoP861//woEDB7B+/XrExcUhJiYGOTk52L17N3Jzc5GSkoI9e/bg4sWL+OY3v2n2\ncolUYWdCZJAXX3wRbrcbLpfL/7kf/OAHqKqqQltbG0aOHInOzk60tLTA6XSat1CiCDBMiAyiKEqf\ntxZXFAXJyclwOp2or683cWVE2nGbi8gggX+j4vTp07h8+TJmzpw54NfduHEDBw8exL333hvN5RFp\nwr9nQmSSY8eOAcCAYfLqq6/ivffew+XLl3Hu3DmjlkakGjsTIoNcu3YN58+f9x8fO3YMDocDDzzw\ngP9zXV1dyMvL8x//+Mc/xoYNGzB79mxD10qkFjsTIgN0dnbi/vvvR2xsLM6fP4+mpiYcPXoUo0eP\nxsiRI/3nbdy4EbNmzTJxpUSRYWdCZIDz58/j0qVLWLRoEW7evIkVK1Zg5cqVuHTpkv9Ptr755puo\nqanBkiVLTF4tkXrsTIgMMGHCBDz77LN47bXXcOLECSxYsABz5szBqFGj8MMf/hBxcXGYNGkS9u7d\ni5gY/o5H4uFTg4kEsGPHDpSWluLmzZtmL4UoKP4KREREmjFMiIhIM4YJERFpxjAhEkBPTw8A4Pr1\n6yavhCg4hgmRhR04cABz5szB0qVL4XA4kJWVhXnz5uHq1atmL42oDz6bi4iINGNnQkREmjFMiIhI\nM4YJERFpxjAhIiLNGCZERKQZw4SIiDRjmBARkWYMEyIi0oxhQkREmjFMiIhIs/8Hhw/ZFHG+nqsA\nAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To explore this a bit, let's consider using the `cvxopt` package (Python ver 2.6 used here). This can be cast as a linear programming problem as follows:\n",
      "\n",
      "min $||\\mathbf{t}||_1 = |t_1| + |t_2|$\n",
      "\n",
      "subject to:\n",
      "\n",
      "$-t_1 < x_1 < t_1$\n",
      "\n",
      "$-t_2 < x_2 < t_2$\n",
      "\n",
      "$x_1 + 2 x_2 = 1$\n",
      "\n",
      "$t_1 > 0$\n",
      "\n",
      "$t_2 > 0$\n",
      "\n",
      "where the last two constraints are already implied by the first two and are written out just for clarity. This can be implemented and solved in `cvxopt` as the following:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from cvxopt import matrix as matrx # don't overrite numpy matrix class\n",
      "from cvxopt import solvers\n",
      "\n",
      "#t1,x1,t2,x2\n",
      "c = matrx([1,0,1,0],(4,1),'d') \n",
      "G = matrx([   [-1,  -1, 0,  0],  #column-0\n",
      "              [ 1,  -1, 0,  0],  #column-1\n",
      "              [ 0,   0, -1,-1],  #column-2\n",
      "              [ 0,   0,  1,-1],  #column-3\n",
      "           ],(4,4),'d')\n",
      "\n",
      "h = matrx([0,0,0,0],(4,1),'d') # (4,1) is 4-rows,1-column, 'd' is float type spec\n",
      "A = matrx([0,1,0,2],(1,4),'d')\n",
      "b = matrx([1],(1,1),'d')\n",
      "sol = solvers.lp(c, G, h,A,b)\n",
      "x1=sol['x'][1]\n",
      "x2=sol['x'][3]\n",
      "print 'x=%3.2f'% x1\n",
      "print 'y=%3.2f'% x2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "     pcost       dcost       gap    pres   dres   k/t\n",
        " 0:  0.0000e+00 -0.0000e+00  3e+00  3e+00  1e-16  1e+00\n",
        " 1:  2.3609e-01  2.3386e-01  5e-01  5e-01  1e-16  2e-01\n",
        " 2:  4.9833e-01  4.9734e-01  5e-02  4e-02  5e-15  1e-02\n",
        " 3:  4.9998e-01  4.9997e-01  5e-04  5e-04  2e-15  2e-04\n",
        " 4:  5.0000e-01  5.0000e-01  5e-06  5e-06  6e-16  2e-06\n",
        " 5:  5.0000e-01  5.0000e-01  5e-08  5e-08  9e-16  2e-08\n",
        "Optimal solution found.\n",
        "x=0.00\n",
        "y=0.50\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example Gaussian Random matrices"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's try out our earlier result about random Gaussian matrices and see if we can reconstruct an unknown $\\mathbf{s}$ vector using $L_1$ minimization."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.linalg\n",
      "\n",
      "def rearrange_G( x ): \n",
      "    'setup to put inequalities matrix with last 1/2 of elements as main variables'\n",
      "    n=x.shape[0]\n",
      "    return hstack([x[:,arange(0,n,2)+1], x[:,arange(0,n,2)]])\n",
      "\n",
      "K=2 # components\n",
      "Nf=128 # number of samples\n",
      "M = 12 # > K log2(Nf/K); num of measurements\n",
      "s=zeros((Nf,1)) # sparse vector we want to find\n",
      "s[0]=1 # set the K nonzero entries\n",
      "s[1]=0.5\n",
      "np.random.seed(5489) # set random seed for reproducibility\n",
      "Phi = matrix(randn(M,Nf)*sqrt(1/Nf)) # random Gaussian matrix\n",
      "y=Phi*s # measurements\n",
      "\n",
      "#-- setup L1 minimization problem -- \n",
      "\n",
      "# inequalities matrix with \n",
      "G=matrx(rearrange_G(scipy.linalg.block_diag(*[matrix([[-1,-1],[1,-1.0]]),]*Nf) ))\n",
      "# objective function row-matrix\n",
      "c=matrx(hstack([ones(Nf),zeros(Nf)]))\n",
      "# RHS for inequalities\n",
      "h = matrx([0.0,]*(Nf*2),(Nf*2,1),'d') \n",
      "# equality constraint matrix\n",
      "A = matrx(hstack([Phi*0,Phi]))\n",
      "# RHS for equality constraints \n",
      "b=matrx(y)\n",
      "\n",
      "sol = solvers.lp(c, G, h,A,b)\n",
      "\n",
      "#nonzero entries\n",
      "nze= array(sol['x']).flatten()[:Nf].round(2).nonzero()\n",
      "print array(sol['x'])[nze]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "     pcost       dcost       gap    pres   dres   k/t\n",
        " 0:  0.0000e+00 -0.0000e+00  1e+02  2e+01  1e-16  1e+00\n",
        " 1:  1.6712e-01  1.6700e-01  1e+01  1e+00  2e-16  7e-02\n",
        " 2:  1.2947e+00  1.2929e+00  4e+00  5e-01  3e-16  3e-02\n",
        " 3:  1.3785e+00  1.3745e+00  2e+00  2e-01  6e-16  8e-03"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        " 4:  1.4705e+00  1.4690e+00  5e-01  7e-02  4e-16  2e-03\n",
        " 5:  1.4976e+00  1.4972e+00  2e-01  2e-02  6e-16  7e-04\n",
        " 6:  1.4979e+00  1.4978e+00  6e-02  7e-03  3e-14  2e-04\n",
        " 7:  1.4998e+00  1.4998e+00  6e-03  8e-04  2e-14  2e-05"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        " 8:  1.5000e+00  1.5000e+00  6e-05  8e-06  3e-14  3e-07\n",
        " 9:  1.5000e+00  1.5000e+00  6e-07  8e-08  2e-14  3e-09\n",
        "Optimal solution found.\n",
        "[[ 0.99999789]\n",
        " [ 0.49999879]]\n"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That worked out! However, if you play around with this example enough with different random matrices (unset the ``seed`` statement above), you will find that it does not *always* find the correct answer. This is because the guarantees about reconstruction are all stated probabalistically (i.e. \"high-probability\"). This is another major difference between this and Shannon sampling.\n",
      "\n",
      "Let's encapulate the above $L_1$ minimization code so we can use it later."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from cStringIO import StringIO\n",
      "import sys\n",
      "\n",
      "def L1_min(Phi,y,K):\n",
      "    # inequalities matrix with \n",
      "    M,Nf = Phi.shape\n",
      "    G=matrx(rearrange_G(scipy.linalg.block_diag(*[matrix([[-1,-1],[1,-1.0]]),]*Nf) ))\n",
      "    # objective function row-matrix\n",
      "    c=matrx(hstack([ones(Nf),zeros(Nf)]))\n",
      "    # RHS for inequalities\n",
      "    h = matrx([0.0,]*(Nf*2),(Nf*2,1),'d') \n",
      "    # equality constraint matrix\n",
      "    A = matrx(hstack([Phi*0,Phi]))\n",
      "    # RHS for equality constraints \n",
      "    b=matrx(y)\n",
      "    # suppress standard output\n",
      "    old_stdout = sys.stdout\n",
      "    sys.stdout = mystdout = StringIO()\n",
      "    sol = solvers.lp(c, G, h,A,b)\n",
      "    # restore standard output\n",
      "    sys.stdout = old_stdout\n",
      "    sln = array(sol['x']).flatten()[:Nf].round(4)\n",
      "    return sln\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example: Sparse Fourier Transform"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As an additional example, let us consider  the Fourier transform and see if we can recover the sparse Fourier transform from a small set of measurements. For simplicity, we will assume that the time domain signal is real which automatically means that the Fourier transform is symmetric."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def dftmatrix(N=8): \n",
      "    'compute inverse DFT matrices'\n",
      "    n = arange(N)\n",
      "    U=matrix( exp(1j*2*pi/N*n*n[:,None] ))/sqrt(N)\n",
      "    return matrix(U)\n",
      "\n",
      "Nf=128\n",
      "K=3 # components\n",
      "M = 8 # > K log2(Nf/K); num of measurements\n",
      "s=zeros((Nf,1)) # sparse vector we want to find\n",
      "s[0]=1 # set the K nonzero entries\n",
      "s[1]=0.5\n",
      "s[-1] = 0.5 # symmetric to keep inverse Fourier transform real\n",
      "Phi = dftmatrix(Nf)[:M,:] # take M-rows\n",
      "y=Phi*s # measurements\n",
      "# have to assert the type here on my hardware\n",
      "\n",
      "sol=L1_min(Phi.real,y.real.astype(np.float64),K)\n",
      "\n",
      "print np.allclose(s.flatten(),sol)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n"
       ]
      }
     ],
     "prompt_number": 141
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(sol)\n",
      "plot(y.real)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 140,
       "text": [
        "[<matplotlib.lines.Line2D at 0x7884910>]"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD9CAYAAAC2l2x5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFjFJREFUeJzt3X9sVfX9x/HXxVuUH5vooK3e262OJtKN/lwRjCI3wkRU\nLnRoUtwUqHaVfyb7Rr/GZIl0yZDFLRvrkqUutIOJSlxYyhywDeflkiCUAYHN6qizxdsrMqH2K13F\n23v5fP8g3Fnbe9t7C/eec3g+EgL33g/n8wLl3Vc+vT9cxhgjAIDjjMt2AADA5cGABwCHYsADgEMx\n4AHAoRjwAOBQDHgAcKikA762tlZ5eXkqKSlJuObpp59WaWmp5syZo7fffvuSBwQApCfpgF+1apV2\n7dqV8PEdO3bo6NGjOnbsmDZs2KCVK1de6nwAgDQlHfBz587Vddddl/Dx7du3a8WKFZKk2bNnq7e3\nV6dOnbq0CQEAaRnTGXw4HFZBQUH8ttfrVXd395hDAQDGzj3WC3z+nQ5cLtew6xLdDwBILt13lBlT\ng/d4PAqFQvHb3d3d8ng8CdcbY2z745lnnsl6his1v52zkz/7P+yefyzGNOD9fr82b94sSdq/f7+m\nTJmivLy8MQUCAFwaSY9oli9frj179uj06dMqKChQQ0ODBgYGJEn19fW65557FAwGVVJSokmTJqml\npSUjoQEAI0s64F966aURL7B+/XqtX7/+kgWyKp/Pl+0IY2Ln/HbOLpE/2+yefyxcZqyHPKPdyOUa\n83kSAFxpxjI7easCAHAoBjwAOBQDHgAcigEPAA7FgAcAh2LAA4BDMeABwKEY8ADgUAx4AHAoBjwA\nOBQDHgAcigEPAA7FgAcAh2LAA4BDMeABwKEY8ADgUAx4AHAoBjwAOBQDHgAcKqMDno9kBYDMyeiA\nHxjI5G4AcGXL6ICPRDK5GwBc2WjwAOBQNHgAcCgaPAA4FA0eAByKBg8ADkWDBwCHosEDgEPR4AHA\noWjwAOBQNHgAcCgaPAA4FA0eABxqxAEfDAZVWVmp0tJSNTY2Dnn8k08+0YoVK1RRUaF58+aptbU1\n4bVo8ACQOe5kD8ZiMdXW1mr37t3yeDyaNWuWFixYoOLi4viaTZs2adKkSTpy5IhOnDihO++8U36/\nXy6Xa8j1aPAAkDlJG3xbW5uKiopUWFionJwc1dTUDGno1157rc6ePauBgQH19PRo4sSJww53iQYP\nAJmUdMCHw2EVFBTEb3u9XoXD4UFrli9frlgspqlTp+r222/Xli1bEl6PAQ8AmZP0iCZRE/+sX/7y\nl3K73Tp58qT+/ve/695779WJEyc0btzQrx2///1ahUIXfu3z+eTz+dIKDQBOFQgEFAgELsm1kg54\nj8ej0MWJLCkUCsnr9Q5aEwwG9cgjj2jixImaPXu2brzxRh0/flwzZswYcr277lqr733vkuQGAEf6\nfPltaGhI+1pJj2iqqqrU0dGhrq4uRSIRbd26VX6/f9Ca+fPn6w9/+IPOnz+vd999Vz09PcMOd4lv\nsgJAJiVt8G63W83NzaqurlY0GlVdXZ2Ki4vV1NQkSaqvr1dNTY3a29tVVVWladOmacOGDQmvxxk8\nAGSOyxhjMrKRy6W1a42eeSYTuwGAM7hcLqU7pnmrAgBwKN6qAAAcigYPAA5FgwcAh6LBA4BD0eAB\nwKFo8ADgUDR4AHAoGjwAOBQNHgAcigYPAA5FgwcAh6LBA4BD0eABwKFo8ADgUDR4AHAoGjwAOBQN\nHgAcigYPAA5FgwcAh6LBA4BD0eABwKFo8ADgUBkd8MZIsVgmdwSAK1dGB/z48bR4AMiUjA74nBzO\n4QEgU2jwAOBQNHgAcCgaPAA4FA0eAByKBg8ADkWDBwCHosEDgEPR4AHAoWjwAOBQGW/wDHgAyIwR\nB3wwGFRlZaVKS0vV2Ng47JqDBw/q9ttvV1lZmXw+X8JrjR/PEQ0AZIo72YOxWEy1tbXavXu3PB6P\nZs2apQULFqi4uDi+pre3VytXrtSf/vQneb1enT59OuH1aPAAkDlJG3xbW5uKiopUWFionJwc1dTU\nqLW1ddCaF198UcuWLZPX65UkTZ06NeH1aPAAkDlJB3w4HFZBQUH8ttfrVTgcHrSmo6NDPT09mjt3\nrioqKrRly5aE16PBA0DmJD2icblcI15gYGBAgUBAu3fvVn9/v775zW/qW9/6liZMmDBkbXv7WvX3\nS8ePSz6fL+l5PQBciQKBgAKBwCW5VtIB7/F4FAqF4rdDoVD8KOaigoICLVq0SPn5+ZKkqqoqBYNB\nLVy4cMj1Zs1aq1tukerqLkV0AHCez5ffhoaGtK+V9IimqqpKHR0d6urqUiQS0datW+X3+wetWbJk\nifbs2aP+/n719PToyJEjuu2224a9HmfwAJA5SRu82+1Wc3OzqqurFY1GVVdXp+LiYjU1NUmS6uvr\nNWPGDK1atUpVVVU6d+6cnnjiCU2ePHnY63EGDwCZ4zLGmIxs5HLpySeNpk6V/vd/M7EjANify+VS\numOaV7ICgENl/L1oOIMHgMygwQOAQ9HgAcChaPAA4FA0eABwKBo8ADgUDR4AHIoGDwAORYMHAIei\nwQOAQ2V8wNPgASAzMn5EQ4MHgMygwQOAQ9HgAcChaPAA4FA0eABI4vXXpbNns50iPTR4AEjiscek\n99/Pdor00OABIIlI5EI5tSMaPAAkMTBwoZzaEQ0eAJKgwY8SDR6A3dDgR4kGD8BuaPCjRIMHYDcD\nAwz4Ubn4bpLGZHJXAEiPMQz40W82TrrqKikazeSuAJCeaPTCzBqX0Ul56WQ8NufwAOwiErHvN1il\nLAx4zuEB2IWdj2ckGjwAJESDTxENHoBd0OBTRIMHYBd2fpGTRIMHgITs/CIniQYPAAnR4FN08cVO\nAGB1NPgUjR/PEQ0Ae6DBp4gGD8AuHN/gg8GgKisrVVpaqsbGxoTrDh48KLfbrW3btiW9Hg0egF3Y\nvcG7kz0Yi8VUW1ur3bt3y+PxaNasWVqwYIGKi4uHrHvqqad09913y4zwTmI0eAB24egG39bWpqKi\nIhUWFionJ0c1NTVqbW0dsq6xsVH333+/pk2bNuKGNHgAdmH3Bp90wIfDYRUUFMRve71ehcPhIWta\nW1u1evVqSZLL5Uq6IQ0egF3YvcEnPaIZaVhL0po1a7R+/Xq5XC4ZY5Ie0axdu1YdHdKLL0rXXeeT\nz+dLOTAAZEo2GnwgEFAgELgk10o64D0ej0KhUPx2KBSS1+sdtObQoUOqqamRJJ0+fVo7d+5UTk6O\n/H7/kOutXbtW77wj3X23xGwHYHXZaPA+3+Dy29DQkPa1kg74qqoqdXR0qKurSzfeeKO2bt2ql156\nadCad999N/7rVatWafHixcMO94s4gwdgF3Y/g0864N1ut5qbm1VdXa1oNKq6ujoVFxerqalJklRf\nX5/yhpzBA7ALR5/BS9K8efN05MiRQfclGuwtLS0jbkiDB2AXdm/wvJIVABKwe4PnvWgAIAEafIpo\n8ADsggafIho8ALugwaeIBg/ALmjwKaLBA7ALGnyKaPAA7IIGnyI+dBuAXdDgU8SHbgOwCxp8imjw\nAOyCBp8iGjwAu6DBp4gGD8AuaPAposEDsAsafIpo8ADsggafIho8ALugwaeIBg/ALgYGGPApocED\nsAuOaFJEgwdgFxzRpIgGD8AuaPAposEDsAsafIpo8ADsggafIho8ALugwaeIBg/ALmjwKaLBA7AL\nGnyKaPAA7IIGnyIaPAC7oMGnyO2WYjHp/PlM7wwAo3f+/IVZ5XZnO0n6Mj7gXS4+eBuA9V18HxqX\nK9tJ0pfxAS9xDg/A+ux+/i5lacDT4AFYnd3P36UsNni+0QrAymjwaaLBA7A6GnyaaPAArI4GnyYa\nPACro8GniQYPwOquiAYfDAZVWVmp0tJSNTY2Dnl8y5YtKisrU1lZmR588EH94x//GHFTGjwAq3N8\ng4/FYqqtrdW2bdt06NAhbdy4UW+99dagNV/96lcVDAZ19OhRLVy4UI8++uiIm9LgAVid4xt8W1ub\nioqKVFhYqJycHNXU1Ki1tXXQmltvvVXXXnutJOnee+9Vd3f3iJvS4AFYneMbfDgcVkFBQfy21+tV\nOBxOuP7555/XkiVLRtyUBg/A6pzQ4JO+jY4rhTdheP311/XCCy9o3759CdesXbtWknTihHTokE/z\n5/tGfX0AyKRsNfhAIKBAIHBJrpV0wHs8HoVCofjtUCgkr9c7ZN2xY8f03e9+Vzt37tSUKVMSXu/i\ngD98WJoxI83EAJAB2WrwPp9PPp8vfruhoSHtayU9oqmqqlJHR4e6uroUiUS0detW+f3+QWvee+89\nLVu2TC+88IKKiopGtSln8ACszgln8EkbvNvtVnNzs6qrqxWNRlVXV6fi4mI1NTVJkurr6/XDH/5Q\nPT09euyxxyRJOTk5amtrS7opZ/AArM4JZ/AuY4zJyEYul46fPi73OLee/J+rdecd1+iRhybpavfV\nmdgeAFKyaZP02mvS5s3ZzeFyuZTumM7oZ5Us2rJI0fNR/bvgU+0MndOaZ/8jSfri1V/U9ROu1/UT\nrlfupFzlTc7TDZNvkPeLXnm/6FXhlEIVTinUxJyJmYwL4ArmhAaf0QH/zvfekSStXi2Vll74ORKL\nqPdcrz765COd+eSM/v2ff+tU3ymd7DuptnCbftf+O534vxM60XtCBx49oLL8skxGBnCFuviJTnaW\nlU8b/OwHb4+/arxyJ+Uqd1Ju0t9z3vAhrgAyJxKhwaclnY/sG+fKyvuiAbhCOaHBZ+3tgnkWDQAr\nc0KD50O3AWAYNPg00eABWJ0TXuhEgweAYTjhaZI0eAAYBg0+TTR4AFZHg08TDR6A1dHg00SDB2B1\nNPg00eABWB0NPk00eABWR4NPEw0egNXR4NNEgwdgdTT4NNHgAVgdDT5NNHgAVkeDTxMNHoDV0eDT\nRIMHYHU0+DTR4AFYHQ0+TTR4AFZHg08TDR6A1dHg00SDB2B1NPg00eABWB0NPk00eABWR4NPEw0e\ngNXR4NNEgwdgdTT4NOXkMOABWBsNPk1XXXXh51gsG7sDQHKxmORy/XdW2VVWBrxEiwdgXU5o71IW\nB/z48XyjFYA1OeH8XaLBA8AQAwM0+DGhwQOwqkiEBj8mNHgAVkWDHyMaPACrumIafDAYVGVlpUpL\nS9XY2DjsmqefflqlpaWaM2eO3n777VFtbLcGHwgEsh1hTOyc387ZJfJnWzr5r4gGH4vFVFtbq23b\ntunQoUPauHGj3nrrrUFrduzYoaNHj+rYsWPasGGDVq5cOaqN7dbgr8T/ya3Cztkl8mdbOvmviAbf\n1tamoqIiFRYWKicnRzU1NWptbR20Zvv27VqxYoUkafbs2ert7dWpU6dG3NhuDR7AlcMpDd6d7MFw\nOKyCgoL4ba/XqwMHDoy4pru7W3l5eUk3vvpqac0a6brr0omdecePS4cOZTtF+uyc387ZJfJnWzr5\nP/rowoyyu6QD3uVyjeoixphR/b7RXs+qjh9vyHaEMbFzfjtnl8ifbenmt/nISj7gPR6PQqFQ/HYo\nFJLX6026pru7Wx6PZ8i1Pv9FAABweSU9g6+qqlJHR4e6uroUiUS0detW+f3+QWv8fr82b94sSdq/\nf7+mTJky4vEMAODyS9rg3W63mpubVV1drWg0qrq6OhUXF6upqUmSVF9fr3vuuUfBYFAlJSWaNGmS\nWlpaMhIcADACc5nt2bPHVFRUmJKSEvOLX/zicm83Zu+9957x+Xzma1/7mpk3b55paWkxxhjz8ccf\nmyVLlpiSkhKzdOlSc/bs2ewGTSIajZry8nJz3333GWPslb2vr888/PDDpry83BQXF5v9+/fbKv/z\nzz9vbr31VlNZWWkef/xxY4y1//5XrVplcnNzzcyZM+P3Jcu7YcMGU1JSYioqKszevXuzEXmQ4fI/\n8cQTZsaMGaaiosI8/vjjpre3N/6YHfJf9JOf/MS4XC5z5syZ+H2p5r+sAz4ajZrp06ebzs5OE4lE\nTFlZmWlvb7+cW47ZyZMnzZEjR4wxxnz44YcmLy/PtLe3myeffNL8+Mc/NsYYs379evPUU09lM2ZS\nP/3pT82DDz5oFi9ebIwxtsr+8MMPm40bNxpjjBkYGDC9vb22yX/mzBlTWFho+vr6TCwWM4sWLTK7\ndu2ydP5gMGgOHz48aMAkyvvmm2+asrIyE4lETGdnp5k+fbqJxWJZyX3RcPn//Oc/m1gsZmKxmHn0\n0Udtl9+YC0Vz4cKFprCwMD7g08l/WQf8vn37zMKFC+O3n332WfPss89ezi0vufvuu8/85S9/MTff\nfLP54IMPjDEXvgjcfPPNWU42vFAoZObPn2/++te/xhu8XbL39vaam266acj9dsnf399vvvKVr5hw\nOGz6+vrMvHnzzP79+y2fv7Ozc9CASZR33bp1Zv369fF1CxcuNG+88UZmww7j8/k/65VXXjHf/va3\njTH2yn///febo0ePDhrw6eS/rO9FM9xz5MPh8OXc8pJ655139Oabb2rOnDk6depU/JvHeXl5o3ox\nVzZ8//vf13PPPadx4/77n9Yu2Ts7OzVt2jStXLlSM2fOVF1dnfr7+22Tf8KECfrVr36lwsJC5efn\n67bbbtPs2bNtk/+iRHnff//9Qc+is8O/51//+tdasmSJJPvkb21tldfrVWlp6aD708l/WQe8nZ/3\n3tfXp5qaGv3sZz/T5MmTBz3mcrks+Wd79dVXlZubq4qKioRPS7VqdkmKRqM6ePCgli1bpoMHD+rT\nTz/VK6+8MmiNlfN/+OGHWr16tdrb29XV1aU33nhDr7766qA1Vs4/nJHyWvnP8qMf/Uhf+MIX9MAD\nDyRcY7X8/f39WrdunRoa/vu8/UT/lqWR8yd9Fs1YjeZ59FY0MDCgZcuW6Tvf+U78q39eXp4++OAD\n5efn6+TJk8rNzc1yyqH27dun7du3a8eOHTp37pw+/vhjPfTQQ7bILl1oJF/60pe0ePFiSdLy5cu1\nefNm5efn2yJ/W1ub5syZo6KiIknSAw88oL1799rm7/+iRHlH+5oXK/jNb36jHTt26LXXXovfZ4f8\n//rXv9TV1aWysjJJFzJ+4xvf0IEDB9LKf1kb/GieR281xhg98sgj+vrXv641a9bE7/f7/dq0aZMk\nadOmTVq6dGm2Iia0bt06hUIhdXZ26uWXX9add96p3/72t7bILkn5+fkqKirSgQMHdP78ef3xj3/U\n/PnztXjxYlvknzt3rv72t7+pp6dHn376qXbu3Km77rrLNn//FyXK6/f79fLLLysSiaizs1MdHR26\n5ZZbshl1WLt27dJzzz2n7du365prronfb4f8JSUlOnXqlDo7O9XZ2Smv16vDhw8rLy8vvfyX+PsF\nQwQCAVNeXm5mzpxpNmzYcLm3G7O9e/cal8tlysrKTHl5uSkvLzc7d+609FPdhhMIBOLPorFT9n/+\n859m9uzZZvr06Wbp0qWmr6/PVvlbWlrMHXfcYaqqqswPfvADE4vFLJ2/pqbG3HDDDWb8+PHG6/Wa\n5ubmpHl//vOfm5kzZ5ry8nITDAazmPyCi/lzcnKM1+s1GzduNEVFRebLX/5y/N/v6tWr4+utmv+z\nf/+fddNNNw16mmSq+V3G8B4CAOBEWftEJwDA5cWABwCHYsADgEMx4AHAoRjwAOBQDHgAcKj/B+Sx\nwNIXyAr5AAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 140
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Uniform Uncertainty Principle"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$\\Phi$ obeys a UUP for sets of size $K$ if\n",
      "\n",
      "<center>\n",
      "$$    0.8 \\frac{M}{N} ||f||_2^2 \\leq || \\Phi f||_2^2 \\leq 1.2 \\frac{M}{N} ||f||_2^2 $$\n",
      "</center>\n",
      "\n",
      "Measurements that satisfy this are defined as *incoherent*. Given that $f$ is $K$-sparse and we measure\n",
      "$y=\\Phi f$, then we search for the sparsest vector that explains the $y$ measurements and thus find $f$ as follows:\n",
      "    \n",
      "<center>\n",
      "    $min_f  \\\\#\\lbrace t: f(t) \\ne 0 \\rbrace $ where $\\Phi f = y$\n",
      "</center>    \n",
      "    Note that the hash mark is the size (i.e. cardinality) of the set. This means that we are looking for the fewest individual points for $f$ that satisfy the constraints.    Unfortunately, this is not practically possible, so we must use the $\\mathbb{L}_1$ norm as a proxy for sparsity.\n",
      "    \n",
      "Suppose $f$ is $K$-sparse and that $\\Phi$ obeys UUP for sets of size $4K$. Then we measure $y=\\Phi f$ and then solve\n",
      "    \n",
      "<center>\n",
      "    $min_f ||f||_1 $ where $\\Phi f = y$\n",
      "</center>    \n",
      "to recover $f$ exactly and we can use $M > K \\log N$ measurements, where the number of measurements is approximately equal to the number of active components.    Let's consider a concrete example of how this works."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example: Sampling Sinusoids"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here, we sample in the time-domain, given that we know the signal is sparse in the frequency domain. \n",
      "\n",
      "<center>\n",
      "$$ \\hat{f}(\\omega) = \\sum_{i=1}^K \\alpha_i \\delta(\\omega_i-\\omega) $$\n",
      "</center>\n",
      "\n",
      "which means that it consists of $K$-sparse nonzero elements. Therefore, the time domain signal is\n",
      "\n",
      "<center>\n",
      "$$ f(t) =  \\sum_{i=1}^K \\alpha_i e^{i \\omega_i t} $$\n",
      "</center>\n",
      "\n",
      "where the $\\alpha_i$ and $\\omega_i$ are unknown. We want solve for these unknowns by taking $M \\gt K \\log N$ samples of $f$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The problem we want to solve is \n",
      "\n",
      "$ min_g || \\hat{g} ||_{L_1}$\n",
      "\n",
      "subject to\n",
      "\n",
      "$ g(t_m)=f(t_m) $\n",
      "\n",
      "The trick here is that are minimizing in the frequency-domain while the constraints are in the time-domain. To make things easier, we will restrict our attention to real time-domain signals $f$ and we will only reconstruct the even-indexed time-samples from the signal. This means we need a way of  expressing the inverse Fourier Transform as a matrix of equality constraints. The assumption of real-valued time-domain signals implies the following symmetry in the frequency-domain:\n",
      "\n",
      "$ F(k) = F(N-k)^* $\n",
      "\n",
      "where $F$ is the Fourier transform of $f$ and the asterisk denotes complex conjugation and $k\\in \\lbrace 0,1,..N-1\\rbrace$ and $N$ is the Fourier Transform length. To make things even more tractable we will assume the time-domain signal is even, which means real-valued Fourier transform values.\n",
      "\n",
      "Suppose that $\\mathbf{U}_N$ is the $N$-point DFT-matrix. Note that we always assume $N$ is even. Since we are dealing with only real-valued signals, the transform is symmetric, so we only need half of the spectrum computed. It turns out that the even-indexed time-domain samples can be constructed as follows:\n",
      "\n",
      "$ \\mathbf{f_{even}} = \\mathbf{U}_{N/2} \\begin{bmatrix}\\\\\\\\\n",
      "F(0)+F(N/2)^* \\\\\\\\\n",
      "F(1)+F(N/2-1)^* \\\\\\\\\n",
      "F(2)+F(N/2-2)^* \\\\\\\\\n",
      "\\dots \\\\\\\\\n",
      "F(N/2-1)+F(1)^* \n",
      "\\end{bmatrix}$\n",
      "\n",
      "We can further simplify this by breaking this into real (superscript $R$)  and imaginary (superscript $I$) parts and keeping only the real part\n",
      "\n",
      "$$\\mathbf{f_{even}} = \\mathbf{U}_{N/2}^R\n",
      "\\begin{bmatrix}\\\\\\\\\n",
      "F(0)^R+F(N/2)^R \\\\\\\\\n",
      "F(1)^R+F(N/2-1)^R \\\\\\\\\n",
      "F(2)^R+F(N/2-2)^R \\\\\\\\\n",
      "\\dots \\\\\\\\\n",
      "F(N/2-1)^R+F(1)^R \n",
      "\\end{bmatrix}\n",
      "+\n",
      "\\mathbf{U}^I_N\n",
      "\\begin{bmatrix} \\\\\\\\\n",
      "-F(0)^I+F(N/2)^I  \\\\\\\\\n",
      "-F(1)^I+F(N/2-1)^I  \\\\\\\\\n",
      "-F(2)^I+F(N/2-2)^I  \\\\\\\\\n",
      "\\dots \\\\\\\\\n",
      "-F(N/2-1)^I+F(1)^I \n",
      "\\end{bmatrix}$$\n",
      "\n",
      "But we are going to force all the $F(k)^I$ to be zero in our example. Note that the second term should have a $\\mathbf{U}_{N/2}$ in it instead $\\mathbf{U}_N$ but there is something wrong with the javascript parser for that bit of TeX.\n",
      "\n",
      "Now, let's see if we can walk through to step-by-step to make sure our optimization can actually work. Note that we don't need the second term on the right with the $F^I$ terms because by our construction, $F$ is real."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def dftmatrix(N=8): \n",
      "    'compute inverse DFT matrices'\n",
      "    n = arange(N)\n",
      "    U=matrix( exp(1j*2*pi/N*n*n[:,None] ))/sqrt(N)\n",
      "    return matrix(U)\n",
      "\n",
      "def Q_rmatrix(Nf=8):\n",
      "    'implements the reordering, adding, and stacking of the matrices above'\n",
      "    Q_r=matrix(hstack([eye(Nf/2),eye(Nf/2)*0])\n",
      "               +hstack([zeros((Nf/2,1)),fliplr(eye(Nf/2)),zeros((Nf/2,Nf/2-1))]))\n",
      "    return Q_r\n",
      "\n",
      "Nf=8\n",
      "F=zeros((Nf,1)) # 8-point DFT\n",
      "F[0]= 1 # DC-term, constant signal\n",
      "n = arange(Nf/2)\n",
      "\n",
      "ft = dftmatrix(Nf).H*F # this gives the constant signal\n",
      "\n",
      "Q_r=Q_rmatrix(Nf)\n",
      "U=dftmatrix(Nf/2) #half inverse DFT matrix\n",
      "feven= U.real*Q_r*F # half the size\n",
      "print np.allclose(feven,ft[::2]) # retrieved even-numbered samples\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "False\n"
       ]
      }
     ],
     "prompt_number": 358
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# let's try this with another sparse frequency-domain signal\n",
      "F=zeros((Nf,1))  \n",
      "F[1]=1\n",
      "F[Nf-1]=1 # symmetric part\n",
      "ft = dftmatrix(Nf).H*F # this gives the constant signal\n",
      "feven= U.real*Q_r*F  # half the size\n",
      "print np.allclose(feven,ft[::2]) # retrieved even-numbered samples\n",
      "\n",
      "plot(arange(Nf),ft.real,arange(Nf)[::2],feven,'o')\n",
      "xlabel('$t$',fontsize=22)\n",
      "ylabel('$f(t)$',fontsize=22)\n",
      "title('even-numbered samples')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "False\n"
       ]
      },
      {
       "output_type": "pyout",
       "prompt_number": 359,
       "text": [
        "<matplotlib.text.Text at 0x7205970>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEgCAYAAACNV7VwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVPX+P/DXsLimgIKIgEJoSSKbGv70qkippUkutyt4\nExXjWubaLdPKxK9r3vx600r9mpCWmmG54p4OkOaSmhtYCiKbe26o7J/fHwcnUQZmYM6cWV7Px8NH\nzMyZmfeccN6ez+d1PkclhBAgIiKSiY3SBRARkWVjoyEiIlmx0RARkazYaIiISFZsNEREJCs2GiIi\nkhUbDZER2djYICMjQ5H3HjFiBKZNm6bIe+siMzMTNjY2KCsrU7oUMjA2GiIroVKpoFKplC6DrBAb\nDZEZKikpqdHzeH42KYGNhhSTl5eHwYMHw8XFBd7e3li8eLHm/gYNGuDmzZuabY8fPw4XFxeUlpYC\nAOLi4uDr6wtHR0e89NJLyMrK0mxrY2OD1atXIzAwEE2aNMHYsWOrrKOq7WNjYzFs2DDN7ceHd0JD\nQzF37lz06tULLi4uGD58OIqKivDmm2/C1dUVzz//PLKzsyu8X3JyMgIDA9GsWTNMnjy5wpd/dZ9r\n1apVCAoKQtu2bQEAW7duRWBgIBwcHNC1a1ecOnWqwj4LDg6Gq6srYmJiqmxOly5dwj//+U+0aNEC\nLi4uiIiI0Dw2YcIEtGzZEo0bN0bHjh3x888/V9g/kZGRGDNmDJo1a4bu3bvj0qVL+OKLL9C6dWu0\natUKP/30k2b70NBQzJ49Gy+88AKaNm2KAQMGVPj//Kjbt29j1KhRcHNzg7u7O6ZNm6bZ71XVSyZI\nECmgtLRUBAcHi/Hjx4vLly+L5ORk0aJFC7Fz504hhBBhYWFi+fLlmu3fffdd8dZbbwkhhNi4caPw\n8PAQu3fvFn/++acYN26c6NKli2ZblUolwsLCRFpamjh69Kho1KiR2LFjh9Zaqto+NjZWvP7665pt\nL1y4IFQqlSgtLRVCCNGjRw/RqlUroVarRXp6uvD29hbt2rUTK1euFNevXxfh4eHif/7nfyq8V2ho\nqLh586bIysoSPj4+4quvvtL5c3Xr1k389ttvoqCgQBw7dkw0btxYbNy4Udy+fVvMmTNHeHl5iaKi\nIlFYWChatmwpPv30U1FSUiLWrVsn6tSpI6ZNm1bpPnj33XfFe++9J+7fvy8KCwvF/v37NY99++23\n4s8//xQ3b94U7733nmjevLkoLCwUQggxffp0UadOHREfH6/5vL6+vmLixIni+vXrYubMmSIsLEzz\nWj169BBubm7izJkz4t69e2LAgAGa/fv4vh0wYICIiIgQFy5cEL/99pvw8/MTy5Ytq7ZeMj1sNKSI\ngwcPCkdHR3H//n3NfRMmTBAjR44UQgjx1Vdfab6gysrKhKenp0hJSRFCCPHSSy+J2bNna553/fp1\nYWdnJ7KysoQQ0hfyDz/8oHm8T58+Yt68eVprqWz7Tz75RAghfZFW1WhCQ0PF+PHjNY/HxMSI4OBg\nze3Vq1eLHj16VHivFStWaG5PnTpVvPDCCzp/rlWrVmkef/PNN0VMTIzmdklJiXB2dhZJSUkiKSlJ\n1K1bVzx48EDzuKenp9ZG884774jXX39dZGZmat1PQghx584dUb9+fXHy5EnN/vH399c8/u233wo7\nOztx8+ZNIYQQOTk5ol69eqKgoEAIITWaYcOGabbfuXOnqFOnjigrK6uwby9fvlzhswshxMKFC0XP\nnj31qpdMA4fOSBEXL15Efn4+WrRoAScnJzg5OSEuLg5Xr14FAAwaNAi//PILLl++jOTkZNjY2OBv\nf/ub5rlz587VPK9169aoU6cOcnNzNa8fGBio+dnNzQ337t0DALRr1w6NGjVCo0aNsH//fq3b5+fn\n6/xZAgICND+7urrC399fc7tZs2YV6nr8vYKCgpCXl6fz5woJCamwD1evXq3Z3tnZGQ8ePMClS5eQ\nl5eHNm3aoF69eprtg4ODtX6GDz74AB4eHujSpQv8/PywatUqzWOffvopnnvuOTg6OqJly5YoLCzE\n9evXNY8/+nldXV3h4eEBR0dHze3CwkLcuHEDgBRIePzzFxcXV3i9h5+trKwM/v7+ms83ffp0XLt2\nrdp6yfTYKV0AWSdPT080atQIly9fRp06dZ543MnJCb1798a6deuQmpqKyMhIzWMtW7bEiBEjMHny\nZJ3fT5TPg5w5c0avOhs2bIgrV65obh8/flyv51fm4dwJABw7dgzu7u4AdPtcdnZ//ZVt2bIloqKi\nsGTJkie2S0pKwrlz5/DgwQPUr19f816PNoVHNW3aFHPnzsXcuXORlJSEvn37on///jh9+jSmT5+O\nlJQUBAcH4+7du3Bzc6txqEAIUWEfHjt2DPb29nB2dtb8YwCQfj9sbGyQmpoKNzc3net1cnKqUV0k\nLx7RkCJCQkLg7e2NDz74AJmZmSgtLcXp06fx66+/arYZOnQoVq5ciR9++AFDhw7V3P/mm29iyZIl\n2LVrF4qKinD79m0kJCRofS99vxSFNKQMQDr6OHjwII4dO4bff/8dX3zxRZWvr8t7rVmzBrdu3UJ2\ndjbWr1+PIUOG1OhzxcTE4Pvvv8fGjRtx79493Lt3D4mJicjPz0eXLl3g6uqKJUuWoLi4GAkJCRUa\n5uMSEhKQk5ODsrIy2NjYoLi4GLdv30aDBg1Qp04dODg44MqVK/jggw9QWFhY7Wesyr59+5CWlob7\n9+9j2bJleO21156IXbu5uaFfv354//33kZaWhrKyMqSnpyM5ObnKesk0sdGQImxsbLB161bk5eWh\nc+fOcHFxwb/+9S/cuXNHs014eDjOnz8PNzc3tG/fXnP/gAEDMHPmTLz77rtwdnZG+/btsXPnTs3j\nj39pVXf+SFXb9+rVC1FRUXjhhRcQGRmJN954o9Ltq3qvx29HRUWhR48eCA4OxoABAxAdHV2jz9Wh\nQwesWbMGc+bMgbu7O9q0aaMZQrK3t8ePP/6Ib7/9Fu7u7ti1a5emoVXm119/RefOneHk5ITY2Fis\nWLECXl5e6NChA8aMGYOwsDB0794dfn5+8PT01OvzPr5/3nrrLYwdOxaenp4oLS3FokWLKt121apV\neOqpp9C7d284OTnhtddew+XLl6usl0yTStT0GJiISE89e/bEsGHDNM2VrINZH9FER0fD1dW1wr92\nHzd16lT4+/ujc+fOOHv2rBGrI6LK8N+21sesG83IkSOxY8cOrY9v27YNJ06cwMmTJ/HZZ59hxIgR\nxiuOiCrFZXCsj1mnzrp164bMzEytj2/evBnDhw8HIE0+37p1C1euXIGrq6uRKiSiR+3bt0/pEkgB\nZn1EU53c3NwKE5ceHh7IyclRsCIiIutj1kc0unh8PFjbYTsP54mIaqa6eTeLPqJxd3evsKBhTk6O\n5uS4SsVKf/pE99GcS2Euf6ZPn654DdZYO+tX/g/rV/aPLiy60YSHh2vOKzh48CAcHR2rnZ/xOeaD\ncRHjjFEeEZFVMOuhs8jISCQlJeH69evw9PTEjBkzUFxcDAAYPXo0+vbti+TkZLRv3x4NGzZEfHx8\nla/XJ6sPxo0dh369+hmjfCIiq8ATNsupVCqdDwNNkVqtRmhoqNJl1Ig51w6wfqWxfmXp8t3JRlPO\n3BsNEZESdPnutOg5GiIiUh4bDRERyYqNhoiIZMVGQ0REsmKjISIiWbHREBGRrNhoiIhIVmw0REQk\nKzYaIiKSFRsNERHJio2GiIhkxUZDRESyYqMhIiJZsdEQEZGs2GiIiEhWbDRERCQrNhoiIpIVGw0R\nEcmKjYaIiGTFRkNERLJioyEiIlmx0RARkazYaIiISFZsNEREJCs2GiIikhUbDRERyYqN5hFLlgAF\nBUpXQURkWdhoHpGYCPj4AAsXAvfuKV0NEZFlYKN5xNat0p/9+4GnnwbmzAFu31a6KiIi88ZG85ig\nIGD9emDfPiAtTTrCmTYNuH5d6cqIiMwTG40Wzz0HfPMNcPAgcPky8MwzwHvvST8TEZHu2Giq0bo1\nsHw5cOIEUFgoNaCxY4GsLKUrIyIyD2w0OvL0BBYtkobTGjYEAgOBN94Azp9XujIiItPGRqMnV1fg\nk0+Ac+cAd3egc2fgn/8EzpxRujIiItPERlNDTZsCM2YAGRlA+/ZAWBgweDBw7JjSlRERmRY2mlpq\n3BiYMkVqON26AeHhQN++wIEDSldGRGQaVEIIoXQRpkClUsEQu6KwEPj6a2DePMDbG/jwQ+loR6Wq\nfY1ERKZGl+9ONppyhmo0DxUXA2vXSid9OjkBH30kHemw4RCRJWGj0YOhG81DpaXADz8As2cDNjbS\nEc6gQdLPRETmjo1GD3I1moeEkJa3mTULuHsX+OADICICsLOT7S2JiGTHRqMHuRvNQ0IAe/ZIRzjZ\n2VKQICoKqFtX9rcmIjI4Nho9GKvRPColRWo4Z84AkydLJ4DWr2/UEoiIakWX707OFCioWzdgxw5g\nwwZg715pxej586WhNSIiS8FGYwI6dpSaza5dwPHjUsOZMQO4eVPpyoiIao+NxoS0by9FovfvBy5e\nlBb0nDoVuHpV6cqIiGrO7BtNcnIygoOD4e/vj8WLFz/xuFqthoODA4KCghAUFIRZs2YpUKV+nnkG\niIuTlrO5cwdo2xaYOBHIzVW6MiIi/Zl1GKC0tBTPPvss9uzZA3d3d3Tq1Alr166Fr6+vZhu1Wo3/\n/d//xebNm6t8LSXCALrKywMWLADi44F//AN4/31p1QEiIqVZfBjg8OHDaN26Nby8vGBvb4+IiAhs\n2rTpie1MtYHoqkULqdH88Qfg7CzN6QwfDpw9q3RlRETVM+tGk5ubC09PT81tDw8P5D42vqRSqXDg\nwAG0a9cOffv2RWpqqrHLNBhnZ+mEz/R0oE0boHt36QjnxAmlKyMia3Tpkm7bmfV56SodFg4LDg5G\ndnY27O3tsXLlSoSHh+O8lquVxcbGan4ODQ1FaGiogSo1LEdHae20iROBZcuAXr2AFSuA/v2VroyI\nLJ1arYZarcaFC1JaVhdmPUdz8OBBxMbGYseOHQCAuXPnwsbGBu+//36l2wsh4OzsjHPnzqFJkyYV\nHjPlOZrqHD4MvPIKsHMnEBSkdDVEZOnOngV69AC++w4IC7PwOZqOHTvi3LlzyMzMRFFREdatW4fw\n8PAK21y5ckWzE7Zs2YL69es/0WTM3fPPA0uWSNfCyclRuhoismTXrgH9+klXGu7ZU7fnmPXQmZ2d\nHeLi4jBw4ECUlJQgJiYGvr6+WLZsGQBg9OjRWL9+PZYsWQI7Ozv4+/tXGhawBIMHSxdfe+UVaWmb\nRo2UroiILE1BAfDqq0BkJDBihO7PM+uhM0My56Gzh4QARo+W4tAbN3JlaCIynLIyYOhQ6Zpaq1f/\ndakTi483U0UqFfDFF0BREfDOO0pXQ0SW5OOPpRXn4+P1v54WG42FsbcHEhKAn34CFi1SuhoisgTx\n8dLyWBs3AvXq6f98Dq5YIAcHIDER6NpVWkGAsWciqqm9e6XrZiUlAS4uNXsNztGUs4Q5mscx9kxE\ntfFojFlbwoxzNFaOsWciqqmaxJi14dCZhRs8WFqyhrFnItJVTWPM2nDorJwlDp099DD2nJsLbNrE\n2DMRafcwxgwAa9ZUnzDj0BkB+Cv2XFwMTJqkdDVEZMoexpi//lr/GLM2bDRW4mHsee9exp6JqHK1\njTFrw0EUK8LYMxFpY4gYszacoylnyXM0jzt0SAoH7NrF2DMR6RZj1oZzNFSpkBBg6VLGnonIsDFm\nbTh0ZqUYeyYiQ8eYteHQWTlrGjp7iLFnIuulb4xZGw6dUZUYeyayXnLEmLVho7FyjD0TWR+5Ysza\ncLCEGHsmsiJyxpi14RxNOWuco3kcY89Elq02MWZtOEdDemHsmchyGSPGrA2HzqgCxp6JLI+xYsza\ncOisHIfO/sLYM5HlMFSMWRsOnVGNMPZMZDmMGWPWho2GKsXYM5H5M3aMWRsOipBWD2PPXbow9kxk\nbpSIMWvDOZpynKPRjrFnIvMiR4xZG87RkEEw9kxkPpSMMWvDoTPSCWPPRKZP6RizNhw6K8ehs+oJ\nAfzrX0BeHmPPRKZG7hizNhw6I4NSqYAvvwSKihh7JjI1phBj1sbEyiFTZ28PrF/P2DORKTGVGLM2\nHPwgvTH2TGQ6TCnGrA3naMpxjkZ/D2PPO3cCwcFKV0NkfYwZY9aGczQkq5AQYMkSxp6JlGCKMWZt\nOHRGtfL3vwMZGYw9ExmTqcaYteHQWTkOndUcY89ExqNUjFkbXb47DdJorl69iqysLGRlZcHW1hbN\nmjVDs2bN4O3tDRul94KO2Ghqp7gY6NsXaNsWWLxY6WqILNdHHwH79gE//WQaCTPZGk1+fj5Wr16N\nFStW4NSpUygsLKx0OwcHB7z44ot46aWXMGTIEDz11FP6vpXRsNHU3u3bUhJt9Ghg/HilqyGyPPHx\nwKxZwMGDppMwM3ijyczMxNy5c7F161YEBASgW7duCAkJgaurK5o0aYKmTZuiqKgIN27cwJ9//onM\nzEwkJycjOTkZubm5mDJlCt5++23Y29vX+sMZGhuNYWRmSs1m2TLGnokMad8+ICJCijG3bat0NX8x\naKPZv38/xo4di2HDhmH06NFo2LChXsWcO3cO8+fPx/79+7F582a0bt1ar+fLjY3GcBh7JjIsU4gx\na2OwRrN27Vr8+uuvmDdvXq2PRjIyMjBmzBjMnz8f/v7+tXotQ2KjMaz164GJE6VDfA8PpashMl/X\nrgGdOwPTpplmwswg59GkpqaipKQECxYsMMiQ19NPP40NGzbg+++/r/Vrken6+9+leZpXXgHu3lW6\nGiLzZG4xZm0Yby7HIxrDY+yZqOZMLcasjewrA6Smptbm6WThuNozUc19/DGQlWWaqzHrq1blDx8+\nvNL7Fy5ciM8//xxlZWW1eXmyAFztmUh/D1dj3rTJNM6VqS2dG01KSgqGDx+O//u//0NaWlqV206a\nNAkdOnRAdHQ08vLyal0kmbeHqz3Pmwds2aJ0NUSmbd8+aTXmxETTOVemtnSeozlw4AAGDRqEq1ev\nAgCaNGkCW1tbfPjhh+jevTsCAwOfeM6dO3cwefJkLF261LBVy4BzNPJj7JmoaqYcY9ZGlpUB0tLS\n8PPPPyMlJQWrV6/WvIGDgwO6du2Kbt26oXv37ujUqRPs7OwwatQorFixouafwkjYaIyDsWeiypl6\njFkbWcIAvr6+iImJwapVqxAcHIyLFy/i66+/xsCBA5GamoqpU6eia9euqF+/PlxcXGT/8k5OTkZw\ncDD8/f2xWMsiW1OnToW/vz86d+6Ms2fPyloPVY2xZ6InWUqMWZtaxZs7deqEI0eOVLgvKysLarUa\nR44cQWBgIIYOHYr69evXutDKlJaW4tlnn8WePXvg7u6OTp06Ye3atfD19dVss23bNnz++efYtm0b\nDh06hAkTJuDgwYNPvBaPaIznYez5t9REOLVZhCIUoq6qLsYPHY9+vfopXR6RUSTuTsSiNYtQKApx\n9lRdPOM0Hupd/cwuYabLd2etzmyYMWPGE/e1bNkSUVFRiIqKqs1L6+Tw4cNo3bo1vLy8AAARERHY\ntGlThUazefNmTTouJCQEt27dwpUrV+Dq6ip7fVQ5lQp4ZVAiVv82AQ+80zX3p38h/cxmQ5YucXci\nJnwxAelB5b//3kDDo+nY/pNl/v7Xqnf27dvXUHXUSG5uLjw9PTW3PTw8kJubW+02ObwcpOK+/H4R\nHrySXuG+9KB0LP6O1xggy7dozaK/mky5jA6W+/tf7RHNjRs3UK9ePb0X0azOxYsX0apVq1q9hkql\n0mm7xw/rtD0vNjZW83NoaChCQ0NrWhpVo1BUfmmJgtICI1dCZHzm/PuvVquhVqv1ek61jaasrAwj\nR47EzJkz8eyzz9a0tgoSEhJw6NAhfPrpp7V6HXd3d2RnZ2tuZ2dnw+OxKNPj2+Tk5MDd3b3S13u0\n0ZC86qrqVnp/PVsLODuNqBol98339//xf4RXNoXyuGqHzlxcXLB06VKMGTMGy5cvR1FRUY0LzMnJ\nwVtvvYUjR47UuskAQMeOHXHu3DlkZmaiqKgI69atQ3h4eIVtwsPDsWrVKgDAwYMH4ejoyPkZEzB+\n6Hj4HPepcJ/tjz6IeHGcQhURGce1a8D5Q+PRLKni77/PMR+Mi7DM33+dU2d3797Fhx9+iA0bNiA6\nOhoDBw5EQEBAtcNX+fn5SElJwffff499+/bh448/RnR0tEGKB4CkpCRMnDgRJSUliImJwfjx47Fs\n2TIAwOjRowEAU6ZMQWJiIho2bIj4+PgKYYGHmDozvsTdiVj83WIUlBagnm09uNmPw/GD/ZCSAjRq\npHR1RIZXUACEhUl//l+Pir//4yLGmWUQQJYTNjMyMhAbG4uNGzdCCIFOnTrBw8MDjo6OcHBw0Fxh\n8+bNm7hw4QJOnjyJFi1a4M0338TEiRPRoEGDWn0oubDRKI+rPZMlM5fVmPVl0EZTXFxc4Xo0RUVF\nSElJwe7du5Geno4rV67g+vXrqFOnDlxdXeHq6oqAgAD06dMHfn5+tfskRsBGYxqKi4G+faVL1Wo5\n/5bILH30kbS47N69lrFQ5kMGazS5ubno2LEjnJyccOLECYNcAM3UsNGYjlu3gK5dgdGjpVUEiMxd\nfDwwa5a09JKlLJT5kMGWoNmwYQPefvtt5Ofn48aNG5r7Z86cWbsKiSrh6MjVnslyWOJqzPrSqdG0\nbNkSKpUKGRkZaN68ueb+zZs3y1YYWTcvL2DDBiA6Gjh2TOlqiGrm7FkgIkJajbltW6WrUY5OjaZ3\n797YtWsXnJyc0LNnT0yZMgUJCQkoLS2Vuz6yYiEhwJIlQHg4wMUcyNxcuwb06wd88on5LPkvF53D\nAA8ePMCCBQuQkJCA06dPa8bkGjduDD8/P7Rv377CH0dHR1kLNzTO0Ziu+fOllA5jz2QuHo0xz5ql\ndDXykiXeDAC3b9/G0aNH8frrr6N37944ffo0UlNTUVDw1/IJ7u7u8Pf3x/PPP49evXohJCQEtra2\n+n8KI2GjMV2MPZM5sdQYszayNZqHHr1MQFlZGc6fP49Tp07h1KlTOHnyJE6ePIkLFy5ACAEHBweM\nGTMG7733nkke7bDRmDbGnslcWGqMWRvZG018fDxGjhxZ5TZ3797FyZMncfz4cRw4cABqtRpxcXF4\n6aWXavq2smCjMX2MPZOps+QYszayNxp95ebmwsbGBqNGjcK2bduM9bY6YaMxD5mZQJcuwLJlQP/+\nSldD9Je9e6UrZCYlWVfCTJZLOddG586d0b179ydWWCbSFWPPZIrOnpWajLXHmLUxaqOZMmUK3Nzc\nEBkZacy3JQvD2DOZEsaYq2fUoTNTxqEz88PYMynNmmLM2pjcHI0pY6MxP4w9k5KsLcasjcnN0RAZ\nkkoFfPklUFQETJqkdDVkbT7+GMjOBr7+2nqbjK64e8is2dsD69dLiZ9Fi5SuhqxFfDywdi2wcaN1\nnCtTWxxsILPn4CCtjNulC+DtzdgzyWvvXmk15qQk6zlXprY4R1OOczTm79Ah4JVXgJ07geBgpash\nS3T2LNCjhxRjZsJMwjkasiqMPZOcGGOuOQ6dkUX5+9+BjAzpyIaxZzKUggLg1VelkzJHjFC6GvPD\nobNyHDqzHIw9kyExxlw1Dp2RVWLsmQyJMeba424ji8TYMxkCY8yGwUEFsliMPVNtMMZsOJyjKcc5\nGsvF2DPpizFm3XGOhgiMPZN+GGM2PA6dkVVg7Jl0wRizPDh0Vo5DZ5aPsWeqCmPMNcOhM6JHMPZM\nVWGMWT7cnWRVGHumyjDGLC8OHpDVYeyZHsUYs/w4R1OOczTWh7FnYoy59jhHQ1QFxp6tG2PMxsOh\nM7JqjD1bJ8aYjYtDZ+U4dGa9GHu2LowxGxaHzoh0wNizdZk2jTFmY+NuJgJjz9YiPl6a+GeM2bg4\nSEBUjrFny8YYs3I4R1OOczT0EGPPlocxZvlwjoaoBkJCgKVLGXu2FIwxK49DZ0SVGDwYSE9n7Nnc\nMcZsGjh0Vo5DZ/Q4xp7NG2PMxsGhM6JaYOzZvDHGbDq4+4mqwNizeWKM2bRwMICoGow9mxfGmE0P\n52jKcY6GqsPYs+ljjNn4OEdDZECMPZs2xphNl9kOnd29exfDhg1DRkYGfHx88M033+Cpp556Yjsv\nLy80btwYtra2sLe3x+HDhxWoliwFY8+miTFm02a2Q2eTJ0+Gs7MzJk+ejE8++QQ3b97EvHnzntjO\n29sbR48eRZMmTap8PQ6dka4YezYtjDEry6KHzjZv3ozhw4cDAIYPH46NGzdq3ZYNhAyJsWfTwhiz\n6TPb/y1XrlyBq6srAMDV1RVXrlypdDuVSoWwsDAEBQVh+fLlxiyRLNijsefp06WhGzKuoiLgP/9h\njNkcmPRBf69evXD58uUn7p89e3aF2yqVCiqVqtLX2L9/P9zc3JCWloa+ffuibdu26NatW6XbxsbG\nan4ODQ1FaGhojWsny+fgAGzfDowdCzz9NPDuu8Do0UDDhkpXZtkePABWrADmzwfatgV27GCM2ZjU\najXUarVezzHbOZq2bdtCrVajefPmuHTpEnr27ImzZ89W+Zx33nkH7u7u+Pe///3EY5yjodr47Tdg\nzhxArQYmTJCaj4OD0lVZlvx8KfW3YAHw/PPAhx9K/yVlWfQcTXh4OFauXAkAWLlyJQYMGPDENvfv\n38fdu3cBANeuXcO2bdvQvn17o9ZJ1iEwEPj+e+kkwT/+AHx8gI8+Aq5fV7oy83frFjBzpnTU+Ouv\n0nlMmzaxyZgTsz2i0RZvzsvLQ0xMDBITE5GRkYFBgwYBAJo2bYp//OMfGD16dKWvxyMaMqSMDOl8\njoQEIDoa+Pe/ATc3pasyL9euAf/9r3QU078/MHUq8OyzSldFj9Plu9NsG42hsdGQHHJypAnrb76R\nzvGYPBlo1UrpqkxbXh7w6adSimzIEGmfeXsrXRVpY9FDZ0TmwMMD+OwzIC1NOrkzOBgYNQo4d07p\nykxPZiaykqpOAAALeklEQVQwZgzg5yedq3TqFLBkCZuMJWCjITICV1dg3jypwbRsKS3QOXQocPq0\n0pUp748/gJEjgQ4dAEdHab2yhQsBd3elKyNDYaMhMqImTaTzbtLTgYAA4MUXgUGDgKNHla7M+E6d\nkoYT//Y3aaL//HkpudesmdKVkaGx0RApoHFj4P33pdBAaKi0TtfLLwP79ytdmfyOHAEGDAB695aG\nEtPTpbP7nZyUrozkwjBAOYYBSEmFhcDKldLwWqtW0jkiL7wgLXdjKVJSgNmzgdRUaYJ/1Cigfn2l\nq6LaYupMD2w0ZApKSoC1a6UhJAcHqeG88or5NhwhgD17gFmzgNxc6YJkUVFAnTpKV0aGwkajBzYa\nMiWlpcCGDdIXNCA1nEGDAFtbZevSlRDAli1S/fn5Uv1DhnCla0vERqMHNhoyRUJIl5GeNQu4fRv4\n4ANpAt1Uv7BLS6XFRufMkZrihx8CAwdyVWVLxkajBzYaMmVCSCtFz5oFXLwoDUENHw7Urat0ZZLi\nYulaMHPmAE2bSsvvvPyy+Q75ke7YaPTARkPmYv9+aVL91ClpxeiYGKBBA2VqKSwE4uOl5Xaeflpq\nMKGhbDDWhCsDEFmgrl2Bbduka7AkJUlf8J98Aty5Y7wa7t2T1iHz8ZHmYlavBn76CejZk02GnsRG\nQ2SmOnQAfvxRSnWdPCl96cfGAn/+Kd973rkDzJ0rvdfPPwObN0tzSF26yPeeZP7YaIjMnJ+fdERx\n4IB0SeM2baQ5HC0Xna2RGzekFQ18fKTzYPbulSb9g4MN9x5kudhoiCxEmzbSlSePHZMixb6+0kXY\ncnJq/pqXL0snVz7zjLSq8i+/SCtRP/ec4eomy8dGQ2RhWrUCPv8cOHNGOjHS31+6xHRGhu6vkZ0N\njB8vNZSCAukKosuXA61by1c3WS42GiIL5eYmXQvnjz+khSqff146Kz8tTftz0tOlFFtgIFCvnjRM\ntmgR4OlpvLrJ8rDREFk4Z2fpUsjp6UDbtlL8+LXXpKOUh1JTgWHDgJAQoEULqTnNnw80b65Y2WRB\neB5NOZ5HQ9bi3j1g2TJgwQJpMr9uXWnBy0mTpAuPNW6sdIVkTnjCph7YaMjaFBRIK0YXFQHR0UDD\nhkpXROaIjUYPbDRERPrjygBERKQ4NhoiIpIVGw0REcmKjYaIiGTFRkNERLJioyEiIlmx0RARkazY\naIiISFZsNEREJCs2GiIikhUbDRERyYqNhoiIZMVGQ0REsmKjISIiWbHREBGRrNhoiIhIVmw0REQk\nKzYaIiKSFRsNERHJio2GiIhkxUZDRESyYqMhIiJZsdEQEZGs2GiIiEhWbDRERCQrNhoiIpIVGw0R\nEcnKbBtNQkIC2rVrB1tbWxw7dkzrdsnJyQgODoa/vz8WL15sxAqNS61WK11CjZlz7QDrVxrrN31m\n22jat2+PDRs2oHv37lq3KS0tRXR0NH788UccPXoUK1asQFpamhGrNB5z/mU159oB860/cXci+ozs\ngxETR6DPyD5I3J2odEk1Yq77/yFzr18XdkoXUFNt27atdpvDhw+jdevW8PLyAgBERERg06ZN8PX1\nlbk6ItOWuDsRE76YgPSgdCATuOh1EelfpAMA+vXqp2xxZHHM9ohGF7m5ufD09NTc9vDwQG5uroIV\nEZmGRWsWSU3mEelB6Vj8neUOL5NyTPqIplevXrh8+fIT98+ZMwf9+/ev9vkqlUqv99N3e1MzY8YM\npUuoMXOuHTD/+qGW/rMTO6GKM7+/B+a+/829/uqYdKPZvXt3rZ7v7u6O7Oxsze3s7Gx4eHhUuq0Q\nolbvRURElbOIoTNtTaJjx444d+4cMjMzUVRUhHXr1iE8PNzI1RERWTezbTQbNmyAp6cnDh48iH79\n+uHll18GAOTl5aFfP2ky087ODnFxcRg4cCA6dOiA6OhoBgGIiIxMJax8zCg5ORkTJ05ESUkJYmJi\nMG7cOKVL0ll0dDQSExPRrFkznDp1Suly9JadnY2oqChcvXoVLi4uGDFiBEaMGKF0WTopKChAjx49\nUFhYiHr16mHIkCGYNGmS0mXprbS0FB07doSHhwe2bNmidDl68fLyQuPGjWFrawt7e3scPnxY6ZL0\ncu/ePYwZMwYnT55EYWEh4uLi0LlzZ6XL0snvv/+OiIgIze2MjAzMnDkT48ePr/wJwoqVlJQIHx8f\nceHCBVFUVCQCAgJEamqq0mXpLDk5WRw7dkz4+fkpXUqNXLp0SRw/flwIIcS1a9eEq6urWe3/e/fu\nCSGEKCgoEO3atRPnzp1TuCL9LViwQAwdOlT0799f6VL05uXlJW7cuKF0GTUWFRUlVqxYIYQQori4\nWNy6dUvhimqmtLRUNG/eXGRlZWndxmyHzgzh0fNs7O3tNefZmItu3brByclJ6TJqrHnz5ggMDAQA\nODs7o1OnTsjLy1O4Kt01aNAAAJCfn4+SkhLUrVtX4Yr0k5OTg23btuGNN94w2zCMudZ9+/ZtpKSk\nIDo6GoA0zO/g4KBwVTWzZ88e+Pj4VDiV5HFW3Wh4no3pOH/+PM6cOWM2QwcAUFZWhoCAALi6umLs\n2LFV/kUzRZMmTcJ//vMf2NiY59eASqVCWFgYgoKCsHz5cqXL0cuFCxc0w8V+fn6IiYnBgwcPlC6r\nRr777jsMHTq0ym3M8zfMQMz9vBlLkZ+fj4iICCxcuBANGzZUuhyd2djY4MSJEzh//jy+/PJLHD9+\nXOmSdLZ161Y0a9YMQUFBZntUsH//fpw4cQJr1qzBnDlzkJKSonRJOispKcGRI0cwePBgHDlyBIWF\nhUhISFC6LL0VFRVhy5YteO2116rczqobjT7n2ZA8iouLMXjwYLz++ut49dVXlS6nRry8vNC3b18k\nJSUpXYrODhw4gM2bN8Pb2xuRkZHYu3cvoqKilC5LL25ubgAAX19fDBw40KzCAB4eHmjatCn69++P\n+vXrIzIyEtu3b1e6LL1t374dHTp0gIuLS5XbWXWj4Xk2yhJCYNSoUWjXrh0mTpyodDl6uX79Om7d\nugUAuHHjBrZv34727dsrXJXu5syZg+zsbFy4cAHfffcdwsLCsGrVKqXL0tn9+/dx9+5dAMC1a9ew\nbds2s9r/zZs3R+vWrXHo0CGUlZUhMTERL774otJl6W3t2rWIjIysfkMjBRNMllqtFoGBgcLPz098\n9tlnSpejl4iICOHm5ibq1KkjPDw8RFxcnNIl6SUlJUWoVCoREBAgAgMDRWBgoNi+fbvSZenk5MmT\nIigoSPj7+4vevXuLr776SumSakytVptd6iwjI0MEBASIgIAAERYWJpYuXap0SXr7/fffRUhIiPDx\n8REDBgwQ+fn5Spekl/z8fNG0aVNx586dare1+vNoiIhIXlY9dEZERPJjoyEiIlmx0RARkazYaIiI\nSFZsNEREJCs2GiIikhUbDZGJO3DggFktr0L0ODYaIhM3aNAgfPXVV0qXQVRjbDREJiwtLQ1Xr15F\n165dlS6FqMbYaIhM2M8//wxAuvYQkblioyEyYSkpKWjSpAl8fX2VLoWoxthoiEzMtm3b0KlTJ3Tq\n1AmrV6+Gra2t5nZqaqrS5RHpjYtqEpmo7OxstGrVCv/9738xfvx4pcshqjEe0RCZqH379gEAunfv\nrnAlRLXDRkNkotRqNRwcHBAYGKh0KUS1wkZDZKLUajVjzWQR2GiITFBWVhYyMzMZayaLwEZDZILU\najWAivMzCxcuVKgaotphoyEyQUeOHNHEmgHgl19+gb29vcJVEdWMndIFENGTAgMD4ejoCDs7O9y5\ncwfLli1DXFyc0mUR1QjPoyEyQQ8ePMDbb78NlUoFe3t7fPzxx2jRooXSZRHVCBsNERHJinM0REQk\nKzYaIiKSFRsNERHJio2GiIhkxUZDRESyYqMhIiJZsdEQEZGs2GiIiEhW/x9Kt7j6wTqyuQAAAABJ\nRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 359
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can use the above cell to create more complicated real signals. You can experiment with the cell below. Just remember to impose the symmetry condition!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Nf=32 # must be even\n",
      "F=zeros((Nf,1))  \n",
      "\n",
      "# set values and corresponding symmetry conditions\n",
      "\n",
      "F[7]=1 \n",
      "F[12]=0.5\n",
      "F[9]=-0.25\n",
      "F[Nf-9]=-0.25\n",
      "F[Nf-12] = 0.5\n",
      "F[Nf-7]=1 # symmetric part\n",
      "\n",
      "Q_r=Q_rmatrix(Nf)\n",
      "U=dftmatrix(Nf/2) #half inverse DFT matrix\n",
      "ft = dftmatrix(Nf).H*F # this gives the constant signal\n",
      "feven= U.real*Q_r*F  # half the size\n",
      "print np.allclose(feven,ft[::2]) # retrieved even-numbered samples\n",
      "\n",
      "plot(arange(Nf),ft.real,arange(Nf)[::2],feven,'o')\n",
      "xlabel('$t$',fontsize=22)\n",
      "ylabel('$f(t)$',fontsize=22)\n",
      "title('even-numbered samples')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "False\n"
       ]
      },
      {
       "output_type": "pyout",
       "prompt_number": 360,
       "text": [
        "<matplotlib.text.Text at 0x73d8f10>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEgCAYAAABvi64JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXlc1WX2xz+XTVwBFUFZRKEFFQTDNDMhN0xG1PpNoaWo\nRe5LpS1OpuaSZeUIWeOMWpo62TbqhHsj4uSamrhNKoIgi4G4IQICz++Px+/1Avdyt+8K5/16+ZJ7\n73O/9/C9X57P95znnPPoGGMMBEEQBCEDDkobQBAEQTQcSHQIgiAI2SDRIQiCIGSDRIcgCIKQDRId\ngiAIQjZIdAiCIAjZINEhCBlxcHDApUuXFPnsMWPGYM6cOYp8tiVkZmbCwcEBVVVVSptCSAiJDkE0\nEHQ6HXQ6ndJmEA0cEh2C0CAVFRU2vY9qwQmlIdEhFCM3NxfPPfccPD090aFDByQlJemfb9KkCa5f\nv64fe+LECXh6eqKyshIAsGbNGgQHB8Pd3R2DBg1CVlaWfqyDgwM2bNiAsLAwtGzZElOmTKnTjrrG\nz5s3D6NGjdI/rhkCioqKwgcffIABAwbA09MT8fHxKC8vx4QJE+Dl5YXHH38c2dnZ1T4vNTUVYWFh\naNOmDd58881qQmDu91q3bh3Cw8Px6KOPAgB++uknhIWFwc3NDU8++SROnTpV7Zx169YNXl5eSEhI\nqFOo8vLy8OKLL6Jdu3bw9PREXFyc/rXp06fD398fLVq0QEREBP773/9WOz8jRozApEmT0KZNG/Tp\n0wd5eXlYsWIFgoKC0L59e/z888/68VFRUVi0aBH69euHVq1aYdiwYdW+Z0Nu3ryJl19+GW3btoWP\njw/mzJmjP+912UuoHEYQClBZWcm6devGpk2bxvLz81lqaipr164d27lzJ2OMsb59+7J//OMf+vEz\nZ85kEydOZIwxtnnzZubr68t2797NioqK2NSpU1mvXr30Y3U6Hevbty87d+4cO3bsGGvevDnbsWOH\nSVvqGj9v3jz20ksv6cdmZGQwnU7HKisrGWOMRUZGsvbt27OUlBSWnp7OOnTowDp37szWrl3LCgsL\nWWxsLHv//ferfVZUVBS7fv06y8rKYoGBgWzVqlUW/15PPfUU++2331hpaSk7fvw4a9GiBdu8eTO7\nefMmW7x4MQsICGDl5eWsrKyM+fv7s48//phVVFSwTZs2MRcXFzZnzhyj52DmzJls1qxZrKSkhJWV\nlbFffvlF/9r69etZUVERu379Ops1axbz9vZmZWVljDHG5s6dy1xcXNiXX36p/32Dg4PZjBkzWGFh\nIVuwYAHr27ev/liRkZGsbdu27MyZM+zOnTts2LBh+vNb89wOGzaMxcXFsYyMDPbbb7+xLl26sJUr\nV5q1l1A3JDqEIhw6dIi5u7uzkpIS/XPTp09nY8eOZYwxtmrVKv1kVVVVxfz8/Nj+/fsZY4wNGjSI\nLVq0SP++wsJC5uTkxLKyshhjfHL+4Ycf9K9HR0ezJUuWmLTF2PgPP/yQMcYn1bpEJyoqik2bNk3/\nekJCAuvWrZv+8YYNG1hkZGS1z1q9erX+8TvvvMP69etn8e+1bt06/esTJkxgCQkJ+scVFRWsdevW\nbN++fWzfvn2sUaNG7O7du/rX/fz8TIrO66+/zl566SWWmZlp8jwxxtitW7dY48aNWVpamv78hIaG\n6l9fv349c3JyYtevX2eMMXblyhXm6urKSktLGWNcdEaNGqUfv3PnTubi4sKqqqqqndv8/Pxqvztj\njC1btow9/fTTVtlLqA8KrxGKcPnyZRQXF6Ndu3bw8PCAh4cH1qxZgz/++AMA8Oyzz+LgwYPIz89H\namoqHBwc0Lt3b/17P/jgA/37goKC4OLigpycHP3xw8LC9D+3bdsWd+7cAQB07twZzZs3R/PmzfHL\nL7+YHF9cXGzx79K1a1f9z15eXggNDdU/btOmTTW7an5WeHg4cnNzLf69evToUe0cbtiwQT++devW\nuHv3LvLy8pCbm4uHHnoIrq6u+vHdunUz+TvMnj0bvr6+6NWrF7p06YJ169bpX/v444/RqVMnuLu7\nw9/fH2VlZSgsLNS/bvj7enl5wdfXF+7u7vrHZWVluHbtGgCezFDz979371614wm/W1VVFUJDQ/W/\n39y5c1FQUGDWXkLdOCltANEw8fPzQ/PmzZGfnw8XF5dar3t4eGDgwIHYtGkTzp49ixEjRuhf8/f3\nx5gxY/Dmm29a/Hns/rrJmTNnrLKzadOmuHr1qv7xiRMnrHq/MYS1FgA4fvw4fHx8AFj2ezk5PfiT\n9ff3x+jRo/HFF1/UGrdv3z5cuHABd+/eRePGjfWfZSgQhrRq1QoffPABPvjgA+zbtw+DBw/GkCFD\ncPr0acydOxf79+9Ht27dcPv2bbRt29bmhATGWLVzePz4cTg7O6N169b6GwOAXx8ODg44e/Ys2rZt\na7G9Hh4eNtlFyAd5OoQi9OjRAx06dMDs2bORmZmJyspKnD59Gr/++qt+zMiRI7F27Vr88MMPGDly\npP75CRMm4IsvvsCuXbtQXl6Omzdv4rvvvjP5WdZOkIyHnQFwr+TQoUM4fvw4fv/9d6xYsaLO41vy\nWRs3bsSNGzeQnZ2N77//Hi+88IJNv1dCQgK+/fZbbN68GXfu3MGdO3eQnJyM4uJi9OrVC15eXvji\niy9w7949fPfdd9XEsybfffcdrly5gqqqKjg4OODevXu4efMmmjRpAhcXF7i5ueHq1auYPXs2ysrK\nzP6OdbF3716cO3cOJSUlWLlyJf785z/XSuVu27YtYmJi8NZbb+HcuXOoqqpCeno6UlNT67SXUD8k\nOoQiODg44KeffkJubi569uwJT09PvPrqq7h165Z+TGxsLC5evIi2bdsiJCRE//ywYcOwYMECzJw5\nE61bt0ZISAh27typf73mBGauPqWu8QMGDMDo0aPRr18/jBgxAq+88orR8XV9Vs3Ho0ePRmRkJLp1\n64Zhw4Zh3LhxNv1ejz32GDZu3IjFixfDx8cHDz30kD7M5OzsjB9//BHr16+Hj48Pdu3apRc3Y/z6\n66/o2bMnPDw8MG/ePKxevRoBAQF47LHHMGnSJPTt2xd9+vRBly5d4OfnZ9XvW/P8TJw4EVOmTIGf\nnx8qKyuRmJhodOy6devQrFkzDBw4EB4eHvjzn/+M/Pz8Ou0l1I+O2eonEwRBWMnTTz+NUaNG6YWW\naHho3tNJTU1Ft27dEBoaqq/zMOTu3buIj49HeHg4IiMjsWXLFgWsJAhCgO5zGzaaTiSorKzEuHHj\nsGfPHvj4+KB79+7o378/goOD9WPWrl2Lpk2b4sSJE7h8+TL69u2L2NhYagdCEApBf3sNG02LzpEj\nRxAUFKSP5cbFxWHLli3VRMfNzQ23b9/GvXv3UFRUhCZNmtBFTxAKsXfvXqVNIBRG0+G1nJycaoua\nvr6+tWoiRowYgcrKSrRu3Rq9e/fGhg0b5DaTIAiCuI+mPR1LPJbPPvsMTk5OyMvLw6lTpxATE4PL\nly/DwaG63pL3QxAEYRvWrNNp2tPx8fGp1kwxOzsbvr6+1cakpqbixRdfRJMmTdCjRw+0a9cO58+f\nN3o8oT5Di//mzp2ruA0N0XayX/l/ZL+y/6xF06ITERGBCxcuIDMzE+Xl5di0aRNiY2OrjenXrx/+\n/e9/o6qqCpcuXUJRUZG+Qy9BEAQhL5oOrzk5OWHNmjUYPnw4KioqkJCQgODgYKxcuRIAMH78eMTF\nxeHs2bOIiIiAp6cnli9frrDVBEEQDRdNiw4AREZG1uqHNX78eP3Pbm5uDUJooqKilDbBZrRsO0D2\nKw3Zry2oI8F9dDqdTfFJgiCIhoy1c6em13QIgiAIbUGiQxAEQcgGiQ5BEAQhGyQ6BEEQhGyQ6BAE\nQRCyQaJDEARByAaJDkEQBCEbJDoEQRCEbJDoEARBELJBokMQBEHIBokOQRAEIRskOgRBEIRskOgQ\nBEEQskGiQxAEQcgGiQ5BEAQhGyQ6BEEQhGyQ6BAEQRCyQaJDEARByAaJDkEQBCEbJDoGRI+NRvLu\nZKXNIAiCqLc4KW2AmtgVsAvpK9IBADEDYhS2hiAIov5Bnk4N0sPTkfRNktJmEARB1Es0Lzqpqano\n1q0bQkNDkZRkXCyOHj2K3r17o2vXroiKijJ7zNLKUpGtJAiCIACNh9cqKysxbtw47NmzBz4+Puje\nvTv69++P4OBg/ZgbN25gzJgx2LlzJ3x9fVFYWGj2uK6OrlKaTRAE0WDRtKdz5MgRBAUFISAgAM7O\nzoiLi8OWLVuqjdm4cSOee+45+Pr6AgBat25d5zEDjwdiatxUyWwmCIJoyGja08nJyYGfn5/+sa+v\nLw4fPlxtzIULF3Dv3j089dRTKC4uxsyZM/Hiiy8aPV7glkD06NIDR385iqbOTS0KxREEQTQkUlJS\nkJKSYvP7NS06Op3O7Jh79+4hJSUFe/bsQUlJCQYMGIBnn30WjRs3rjX24vGLUphJEARRb4iKiqp2\nQz5//nyr3q/p8JqPjw+ys7P1j7Ozs/VhNAE/Pz8888wz8Pb2RseOHREREYHU1FS5TSUIgiCgcdGJ\niIjAhQsXkJmZifLycmzatAmxsbHVxgwdOhT79u1DSUkJioqKcOLECTz55JMKWUwQBNGw0XR4zcnJ\nCWvWrMHw4cNRUVGBhIQEBAcHY+XKlQCA8ePH49FHH8XYsWMRERGB0tJSzJw5E82aNVPYcoIgiIaJ\njjHGlDZCDeh0OtCpIAiCsA5r505Nh9cIgiAIbUGiQxAEQcgGiQ5BEAQhGyQ6BEEQhGyQ6BAEQRCy\nQaJDEARByAaJDkEQBCEbJDoEQRCEbJDoEARBELJBokMQBEHIBokOQRAEIRskOvWA+HggN1dpKwhC\n2+zeDSxdqrQV9R8SHY3DGPCvfwH/+5/SlhCEtjl6FNizR2kr6j8kOhrn+nXg9m3ydAjCXnJzgcxM\npa2o/5DoaBzhj4REhyDsIzcXuHwZqKpS2pL6DYmOxiHRIQhxyM0FysqAq1eVtqR+Q6KjcTIzgTZt\nSHQIwl5yc/nfEoXYpIVER+NkZgK9epHoEIQ9VFUB+fnAE0+Q6EgNiY7GuXyZi05entKWEIR2KSwE\nWrQAHnmEREdqSHQ0TmYmvzvLzeXp0wRBWE9eHtCuHRAQwG/kCOkg0dEwjHHR6dwZcHEBbtxQ2iKC\n0Ca5uUDbtlx0yNORFhIdDSOIjLs7v0ujdR2CsI3c3AeeDomOtGhedFJTU9GtWzeEhoYiKSnJ5Lij\nR4/CyckJP/74o4zWSUtmJv8j0elIdAjCHgTRad+eh9coVC0dmhadyspKjBs3Dj/++COOHTuG1atX\n49y5c0bHvfXWWxg0aBBYPbqaBNEBSHQIwh4E0WnShCcUUK2OdGhadI4cOYKgoCAEBATA2dkZcXFx\n2LJlS61xSUlJ+L//+z94enrWeTyt6VFmJr8zA0h0CMIeBNEB+N8UhdikQ9Oik5OTAz8/P/1jX19f\n5OTk1BqzZcsWTJw4EQCg0+lMHk9rC/E1PR1KmyYI2zAUHVrXkRYnpQ2wh7oERGDGjBlYsmQJdDod\nGGN1htfeeWcevL35z1FRUYiKihLJUmnIzAQiI/nP7doB+/crag5BaBYhZRog0TFHSkoKUlJSbH6/\npkXHx8cH2dnZ+sfZ2dnw9fWtNubYsWOIi4sDABQWFmL79u1wdnZGbGxsreNFR8/D8OHS2iwmly8/\n8HTatqXwGkHYQmUlX8Px8uKPAwKAtDRFTVI1NW/I58+fb9X7NS06ERERuHDhAjIzM9GuXTts2rQJ\n//znP6uNuXTpkv7nsWPHYsiQIUYFB9De3Q0lEhCE/RQUAB4evNYN4H9T//63oibVazQtOk5OTliz\nZg2GDx+OiooKJCQkIDg4GCtXrgQAjB8/3qrjaUl0btzg/aI8PPjjtm15iIAxnkJNEIRlGK7nABRe\nkxodq085xHag0+kQG8tgJPlNlfz2GzB6dPUwQMuWwPnzQOvWytmlBa5fB37/HejZU2lLpGXHDmDg\nQMBB0+lC0vPTT8DnnwPbtvHHd+4Anp78f7qBM4+wXm4pdDkaoKW7G8PQmgCF2Cxjxw4gNhYoKVHa\nEulISwOeeQbIylLaEvVT09Np2hRo1gz44w/lbKrPkOgYkJmpnVodwxodAUqbtoyCAv5v1SrTY5J3\nJyN6bDSixkQhemw0kncny2dgHVhq1wcf8P8LCmQ0TqPUFB2AanWkRNNrOlJw48aDdRI1Q56O7RQU\nAH/6E7B0KTBhwoMFZIHk3cmYvmI60sPT9c+lr+A/xwyIkdPUalhq14ULwJ49fMsLEh3z5OYC3bpV\nf05Y1+nRQwmL6jfk6RigpQVEY6JDadOWUVAADBoEdOoEfP117dcTNyZWm9gBID08HUnfmO7tJweW\n2vXRR8CkSUDHjiQ6lpCXx/92DNHSXKA1SHQM0NKFZlijI0CejmUUFvKF4tmzgSVLeJ2GIWWszOj7\nSitLZbDONJbYlZ0N/PADMG0a/x0LC+WyTrsYC69paS7QGiQ6BmjpQqPwmu0UFPAJuU8fXhD4/ffV\nX2+ka2T0fa6OrjJYZxpL7PrkE2DcOKBVK/47kqdjHhIdeSHRMUArF9qNG0BFBU+RNoRExzIE0dHp\nuLezeHH1BJJpI6ch8ERgtfcEHg/E1LipMltaHXN2FRQA69YBr7/OXyPRMU9FBfcGhW4EAlqZC7QI\nJRIYEBAA2NFSSDaE0FrNGgISHcsQRAfgacV/+QuQnMyTC4AHi/JJ3yShtLIUro6umDplqqJJBJbY\n9de/Ai+88OCunUTHPFev8ro2pxozoeG+OlSrIy4kOgZo5e7GWGgNALy9+R9RVRUVBJqispIXh7Zq\nxR8L3s6iRUBMzIMJJmZAjOIiYwxTdt28Cfztb8DRow+eI9Exj7HQGsDrdJo25eevTRv57arP0NRk\ngCA6aq/VMVajAwCNGgFubupZPD50iHdIUBNFRfwcGd7ZPvssf14LXq4pPv8cGDyYZ6wJqFV0UlLU\nU7RqSnQAqtWRChIdA9zd+f9q31fHlKcDqCfEVlkJjBrFa2HUhGFoTcDREXj7bb62o0VKSoDly4F3\n3qn+vFpFZ+rU2rYqheGWBjXRSuRDa5DoGKDTaeNCq0t01FKrs2ULF55t29TlORoTHQB46SXulR05\nIr9N9rJqFS8E7dSp+vNubsDdu0CZ8UxrRcjK4tfnjh18zURpcnNr1+gIaGEu0CIkOjXQwoVmrEZH\nQC2ezscf8yLFJk2AkyeVtuYBBQXGG6I6OwOzZmnP2ykv597k7Nm1X9Pp+NqVWsKtAL8JGTwYGDuW\ne2dKU1d4TQtzgRYh0alBQIA67sDqQu3htV9+4c0Shw/nE4zQvVcNmPJ0AODll/k61OnT8tpkD+vX\ncw8nIsL462oLsQmiM3068NVXyoeySXTkh0SnBmq/0G7e5He3QvZVTdQgOkuX8loRR0f1iY7QjcAY\njRsDM2Y8aJapdioreUcFY16OgJpEp7SUJxFERwN+fjxF/f7WV4pBoiM/JDo1UHvGiqkaHQGlO03/\n/jtw4AAwZgx/HBnJ2+wXFSlnkyF1eToA71m2cyeQnm56jFr4/vsHnRVMoaZWOPv2AaGhD4qa33gD\nSExUds3JXPaaUKtDiAeJTg3UfndTV2gNkN7TMdda/9NPgYkT+VoOALi6AlFRfCJXA+ZEp0ULbv9H\nH8lnky0wxtef/vKXuosX1eTpCKE1ga5dgS5dgI0bq4+Ta1uJ8nJes2XqemjenHu/ahHt+gIVh9ZA\nC6JjrEZHQErRMdda/+pV4NtvubdjiBBiGzFCGruswZzoAHy94eGHgffeA3x85LHLWpKTeQHwM8/U\nPU5tovPdd9WfmzmThzTHjOHiKee2Elev8vY3jo6mxwiRD3PXDGE55OnUwMODV/QrvcBpCnOejpcX\nX8Sv2TlZDMy11l+xgrdhqVnB/cwzPEVWCpusxRLRad2aT4KffCKLSVbDGO+gMHu2+RYtahGdCxf4\n9s9du1Z/vn9/vp/R9u38sZzbStSVLi2g9ptQLUKiUwO11+qYEx1nZx4zl2Kr3bpa69+5A3zxxYNm\nk4a0b8/F8NdfxbfJWiwRHYCvN3z1lTpDK/v28TWyZ581P1YtoiOE1mqKpE7HvZ2PP+aP5dxWoq71\nHAE1zwVahUTHCGq+0Oqq0RGQKsRWV2v9L78EevfmYSljxMQon8XGWN3Za4b4+PCCy337pLfLWrZt\n490e6goLCahFdJKT+TVgjOefBy5eBI4dk3dbCRIdZSDRMYKaLzRzng4gneiYaq0/+fmp+PRTXlxp\nisGD+cSjJDdv8sSGRsbntVoEBQEZGdLaZAuXLgEPPWTZWDWITnExcPAg0K+f8dednfm6ztKl8m4r\nQaKjDJRIYAS1Xmi3bvFaB2MV9YZIlTZtqrV+yY0YtG3LPQNT9OrF05Dz83k3bCUw1Y3AFB071k6K\nUAMZGdUbe9ZF69bKi85//gM8/jjPDDRFQgJfp1qyJAbLJ8uzrURuLvfO60Ktc4GW0bzopKamYsaM\nGaioqEBCQgKmTq1+R7RhwwZ8dD//tXPnzpg9eza6dOlS5zEDAoDUVKksth1zNToCUmaw1WytzxjQ\no4f5Bo7OznzReMeOBzU8cmPpeo5Ahw4PFrjVxKVL3DZLaNWKJ8VUVloWjpOCmqnSxmjeHHjlFWDZ\nMmD5cnm2lbDE0xGy12hfHfHQdHitsrIS48aNw48//ohjx45h9erVOHfuXLUxHTt2RGpqKk6ePIno\n6Gi88sorZo+r1rsbS0JrgLxdCfbv55NabKz5sUqv61i6niPQsaP6wmvXr3MBMdWRoiaOjrx7+rVr\n0tplCsbqXs8xZPp04Ouv5SsktkR0WrTgIVk1JpRoFU2LzpEjRxAUFISAgAA4OzsjLi4OW7ZsqTbm\niSeegJubGwAgJiYGV65cMXtcNYtOXTU6AnJ2ml66lGd6WXIXPWgQsHs3cO+e9HYZwxZPJzOTp9Cr\nBSG0Zs1dt5LrOqdPcy/3kUfMj23XDhg6lG9GJwd5eeZTpoEHnQkIcdC06OTk5MDPz0//2NfXFzk5\nOSbH//3vf8fQoUPNHlettTpq83TOnuVbAYwebdl4b28gMJAvKiuBtaLTpAm/FpTuZWeINaE1ASVF\nx1SqtCneeANISuJrl1JSVsbXSC1Z41PrTahWEWVN548//kBWVhaysrLg6OiINm3aoE2bNujQoQMc\nJNw3WWfF7d7evXuxfv16HDhwwOSYefPm6X9u3ToKmZlRCAuzx0Jxyczk6yfmkEt0Pv0UmDyZtwqx\nFCGLra5+YVJRUGB9EkOHDty78PWVxiZrsSaJQEDJ/mvbtgFvvWX5+C5dgPBwYMMG3vVbKvLy+LVg\nyfREolOdlJQUpNixza5NolNcXIwNGzZg9erVOHXqFMpMdOxzc3ND//79MWjQILzwwgto1qyZzYYa\nw8fHB9nZ2frH2dnZ8DUyO6SlpeHVV1/F9u3b4S5sD2oEQ9E5fpxfaGoTHUs8nTZteAy/oqL6tsxi\nkpcH/PADrzS3hsGDeabShx9KY1ddFBQAISHWvadjR+5dPPWUNDZZy6VLfGK2BqU8nRs3gBMneO89\na5g1izdeHTvWMlGwBUvWcwQCAtS37bqSREVFIcrgS50/f75V77fqK83MzMT48ePxyCOPYMuWLRg+\nfDiSk5Nx6tQp5OTkoLS0FLdu3UJGRgaOHTuGNWvWwNfXF59//jmCgoLw17/+FfdEDOhHRETgwoUL\nyMzMRHl5OTZt2oTYGivaWVlZeO6557B+/XoEBQVZfGw13t1YUhgKcKHx9OS9paQiKQkYOdK6FGQA\n6N6dp01nZUljV11YG14DHoiOWtBSeG3XLi7WQvNXS4mK4u+Rsq7LWtFR21ygZSy+D/7ll18wZcoU\njBo1Cp9++imaNm1qdJyLiwuaNWuG9u3bIzw8HMOHDwcAXLhwAR999BG6du2KrVu3WiUAJo13csKa\nNWswfPhwfcp0cHAwVt7fpGP8+PF4//33UVRUhAkTJgAAnJ2dccSCPYnVdqHdvg2UlFg+aQohNika\nVpaUAH//O3D4sPXvdXTkCQXbtwPjx4tvW13YIjodOgB790pjjy3YGl6z1iMVA0tSpY2h03Fv55NP\ngCFDxLcLINFRFGYBGzduZK+//jorLy+3ZHidpKens+joaHby5Em7jyUmNU/FDz8wNnSoQsYY4dQp\nxoKDLR8/ZAhjmzdLY8tvvzHWubPt79+4kbHYWPHssRR/f8YuXbLuPSkpjD35pDT2WEtFBWONGjF2\n965179u4kbEXXpDGJlNUVjLWpg1j6em2vb+4mDFXV8aqqsS1S+CttxhbtMiysTduMNa0qXS2aB0L\nZUSP2fDa2bNnUVFRgU8++QTOzs52i1zHjh3xr3/9C99++63dx5IStd3dWLqeIyBl2nRBQe1O0tYw\ncCDfQVLqDKWaWNuRAFBXeC0nh9fnuFrZhkyJrgTHj/PGs9Z6ZQJNm3KP584dce0SsDRdGgDc3Hgn\nbKVqneobZsNrnTp1QqdOnUT90MaNG2PhwoWiHlNstC46Umaw2TJ5G9KqFV8MX5qYjP+eS0QZK0Mj\nXSNMGzlNskr0O3d4Gry1uSzt2vFixbt3rcvSkwJbQmuAPGs6ybuTkbjxwXfZWjcNgwfb910KYily\n/hEA68JrwIP5wJ7rnuDYldt09uxZ0QVJLRjW6tSR8CYblhaGCrRrBxw6JI0ttqyN1CSwUzI+2Tod\nNwdIv1kX8KAbgbWtTBwdAX9/fv6Dg0U3yyouXVKn6BjbeK3R1nS8NwIAbP8uBbutTZywBGtFR2iH\nExEhvi0NDbsSEuPj440+v2zZMnz22WeoUlMpt5WobV8dtXk69orOxeuJ1QQHkG6zLsA+m9USYrMl\ncw3gd+eFhbwljRQY23itLDYdKWfs+y6lFEtbPB3qSiAOFovO/v37ER8fj7///e+1+pvV5LXXXsNj\njz2GcePGIVdN5dxWItzdqAFbREeKTtOAOKLj0ky+zboA+2wWCkSVxtbwWqNGPAVZqg4bpjZeK2f2\nfZdSiU53dYw4AAAgAElEQVRJCQ+Xtmxp+XvUdAOqdSwWHUdHR+zcuRMTJkxA586d0bp1a2RlZSEx\nMRG//fZbrfFPPPEEEhMT8f7774tqsJyo6UKztEZHQEpPx9rGmcaQc7MuoP54OrYuzEvpNUj1XUrV\nSUFIIrAm1KqmuUDrWCw6vXr1Qn5+Ps6cOYOVK1di8ODBKCwsxIwZM9CtWzd4eHjgT3/6Ez788EMc\nPHgQFRUVaNGihajFoHKjFpe6uJgvhFuTMebpye9sy8vFt0cMT2fayGnw3i/PZl1A/REdW9c3pGyF\nY2zjtfZH7P8upRJKa0NrAImOmFidSBAcHIzg4GAkJCTg3Llz+PHHH7F37159P55t93vXOzo6wt3d\nHUOkqu6SgYAA4L//VdoKLnzt21t3Z+bgAHh58ep/f39x7RFDdGIGxGB5CTBiahKejCpFE2fpNusC\ntB9eKynhO59amuZbEyk9HcPN/a4WlSL9nCtWJNn/XUpV1GqL6NC+OuJhd2cuPz8/jB49GqPvtxrO\nyspCSkoKjh49irCwMIwcOdJuI5VCLXc31q7nCAi1OmoUHQB4fmgMVibGYMafpas8Fygo4B2ubUHw\ndJSccDIy+DVgay8yqTPYhM39Zs8GdJ2AmAH2H1Mqm62p0RFwd+ftpYqKLN/LiDCOXaJjrNGbv79/\nNRHSMloXHSnWdSor+UZiYv3hDR7M26XIITq2CqW7O98TRoy1LFuxJ7QGyNd/bds2YMUKcY6lpvAa\n8GA+INGxD7tSpgfb0lhJQ7RsyTs1K72vjrU1OgJSZLAVFfEKbbG6VwtbHUiVzitgb0Gr0iE2WzPX\nBOToSpCTA2RnW7b9hiWoTXTUlM2qZcyKzrVr13BHgl4Ul9WwQm8GoVZHaVPV5OnYO3nX5NFHefqq\nVOndAvaGBJVOJrAncw2Qx9M5fBjo1Uu8GxKphNJeT4ewD7OiU1VVhbFjx+L3338X7UO/++47JCVJ\nUwQoNmq40NQmOmKGmHQ6vsfN6dPiHdMY9UF01B5eO33a+v2K6sLdnd+QmNiuy2bsER2lb0DrA2ZF\nx9PTE3/7298wadIk/OMf/0C5HTm4V65cwcSJE3H06FF8/PHHNh9HTtQgOtbW6AhoQXQA3oft1Clx\nj2lIeTnP/rKnnZHWw2tyiM6pU9ZvMFcXOt2DbgpiQp6Osli0ptOyZUts3rwZp06dQmBgIObOnYvf\nfvsNzIJAfHFxMbZv346xY8eid+/e6N69Oz766CO7DZcLpS+0O3f4XjpeXta/VwrRkWIxXWpPp7CQ\nL/7aswulkp4OYw3T0wHEt7u4mK/TurlZ/16l54L6gsXR1+bNmyMxMREzZszAvHnz0KdPHzDG0L17\nd/j6+sLd3R1ubm4oLy/HtWvXcP36dWRkZCAtLQ3t2rXDhAkTsGLFCjSxdhtBhVG6VseWGh0BKbY3\nkMLTCQkB7u+7Jwli2Kyk6PzxB+9w3aKF7ccQJm+p0r5LS/mE/Mgj4h5XbNGxpRuBgCA6VKtjHxaL\nzr179+Ds7IyOHTti3bp1KC8vx/79+7F7926kp6cjMzMThYWFcHFxgZeXF7y8vPDiiy9i3bp16CKm\nzy0zSt/d2LqeA/C7+9u3+YRg7R4spigosC/MY4zOnYFz53g6tqOjuMcGxBEdf38u4Pfu8fRpObE3\ntAZU359Giq0Czp3jdVAuLuIeV2zRsTW0BvDwrIMDLxmwpm8bUR2LRCcnJwcRERHw8PDAyZMn4ezs\nDBcXF/Tr1w/9+vWT2kZF0bLoODjwu7r8fNuPUZOCAvFSYgWaN+ctfi5dAh56SNxjAyI1KHUBvL15\nSrDYomsOezPXBIRWOFKIjhShNUBdogM8mA9IdGzHoij3v/71L0yePBnFxcW4ZrB93oIFCyQzTC0o\nXatja42OgNjrOlKE1wA+YUmVTCCWzUqF2OxdzxGQcl3n1KmGITpUq2M/FomOv78/dDodLl26BG9v\nb/3zW7dulcwwtaB0rY49ng6gHdGRMoNNTNFRIoNNjPAaIL3oSBFFV5voKB35qA9YJDoDBw7Erl27\n4OHhgaeffhpvv/02vvvuO1RWVkptnypQ8kJTo+hIsWWvlBlsYtncoYNyno4YoiNlVwKpwmti20yi\nozwWrem4urpix44d+OSTT/Ddd99h6dKl+nRpd3d3dOnSBSEhIdX+uathj2eRINHhMCZd/7GQEMBI\nKz9RENPT2bLF/uNYi9rDa9ev8/CzPWFgU6jR00lJEcuahonFlQuNGzfGu+++i5MnT6KoqAh79uyB\nt7c3hg0bhtLSUqxduxZTpkxBZGQkWrZsCT8/P8TExGD+/Pk4cOCApr0iJUQneXcy+sdH41rzKIyZ\nHY3k3ck2HUfMtOmbN/kulGJlwhny8MM8hHn3rvjH1nJ4rbycJ4L4+dl/LKlE5/RpnoFoTx2UKcTe\nByg31/btIQAg+2oy/vO/aESNiUL0WNv/LhsyNnVJcnNzQ9++feHj44OvvvoKAG+Xc/HiRZw6dQqn\nTp1CWloa0tLSsGPHDsyfPx9ubm6YNGkSZs2aJaoXlJqaihkzZqCiogIJCQmYOrX2xlHvvPMOkpOT\n0aRJE3z11Vd49NFHrfoMuWt1kncnY/qK6Xzf+Y7ALgDpK/ge9NbuUSKmpyPVeg7As8OCgoD//Q8I\nDxf32GJ5Z0qE17KyAB8fcdK0pdqfRqrQGiCuUDLG63Rs9XSSdyfjr8nTUTwyHfvuP2fr32VDxq7W\nfJMmTdL/7ODggIcffhgPP/wwnnvuOf3zt2/fRlpaGk6cOIEDBw6gU6dOWLNmDQYNGmTPRwMAKisr\nMW7cOOzZswc+Pj7o3r07+vfvj+DgYP2Ybdu24eTJk0hLS8Phw4cxZswYHDp0yKrPkdvTSdyYyAXH\ngPTwdCR9k2ST6IjVTFPq1v5CMoHYoiOWWLZpwz2xW7fsK9S0BrFCa4B0no5USQQArzW7cUOcGq7b\nt/n/zZvb9v7EjYnIiBDn77IhY5dDPHbsWLNjmjdvjieffBJTpkzB0qVLcezYMSQmJtrzsXqOHDmC\noKAgBAQEwNnZGXFxcdhSI+i+detWxMfHAwB69OiBGzdu4OrVq1Z9jtxpkmXMeIfD0spSq4+lFU8H\nkCZtWsz9f3Q6+UNsYmWuAdKKjlSejqMjL8o0qNSwGWE9x9ZuAmL+XTZkJIjCmqZnz57o06cPfH19\nRTleTk4O/AyC3b6+vsjJyTE75sqVK1Z9TqtWvBJdrlqdRrpGRp93dbR+McXDg9+dl5TYa5U8oiN2\nBpvY+//IHWITK3MNkEZ0GJM2vAaIZ7e9SQRi/l02ZET6U7SMt99+G5s2bcKIESNEOZ7OwluWmo1J\nTb1v3rx5+p+joqIQFRV1f/yDWh05kvKmjZyGE4vTURD1wJUPPB6IqVNqr1eZQ6d7EGKzdbtmAalF\nR4paHbFtlrtA9NIlwCBabRdSiE5ODl+Pk/K6UIvoTBs5Dekr0quFvm39u9QyKSkpSLEjhU9W0Zk8\neTImT54s2vF8fHyQnZ2tf5ydnV3Li6o55sqVK/Dx8TF6PEPRqYmwrtO1qz0WW8Yz/WLgMh7o/lsS\nmriVwtXRFVOnTLU5biyE2MQQHYPaYNFp356vl1y/zj00MZBCdM6fF+945hAzvObmxvvwlZXxLEQx\nkDK0JqAW0RH+/pK+ScL5jFKU3nbF8sW2/11qFcMbcgCYb2Wtg6yiIzYRERG4cOECMjMz0a5dO2za\ntAn//Oc/q42JjY3FZ599hri4OBw6dAju7u7wsmGfADmTCf79b6Btqxgc/jFGlG62YqVNFxRIt2AM\n8JTbzp15uOapp8Q5ptii06EDsGOHeMczh5jhNcP9aUzcd1mN1KE1QLwC0dxc+1PPYwbEIGZADG7c\n4N9LqHWJsARkXtMRGycnJ6xZswbDhw/HY489hnHjxiE4OBgrV67Eyvu98gcPHozQ0FCEhIRgxowZ\n+PLLL236LDlFZ+lSYNYs8dqni5VMIHV4DRA/xKbl8NqNG3wtUYwkCAGxQ2xSZq4JiFWrY0+6dE3c\n3YExY4Dly8U5XkNC054OAERGRuLEiRPVnhs/fny1x0uWLMGSJUvs+pyAAODAAbsOYREHD3KBePZZ\n8Y4pVtq0HKIjdgab2G17OnTga3tVVdIUQxoihNbE3LtF7LYyp04BU6aIdzxjeHoCFy/afxx7w2s1\nmTGDp/fPmWPbpnANFU17OnIil6ezdCnw+uviZVsB2vJ0xM5gE9vmJk34BCNW7VNdiBlaExDT06mo\nAH7/nYdEpUQtazo18fcHnnlG2g0IAZ4heO6ctJ8hJyQ6FhIUxO+2yoyn6ovC+fO884EF5U9WIZbo\nSF0cCjwIr1mwE7pFSGGzXCE2MQtDBcQUnYsX+Xph06biHM8UYtjMmP0tcIwxaxYPsZWXi3tcQ06c\nAIYPl+74ckOiYyHu7vyOLjVVus/49FNgwgTx/4jFEJ2SEl5oKcUGYIZ4evLebjXKrWxGCu9MrgJR\nMTPXBMQUHTmSCABxbL5xg6d2i339du3K54WNG8U9riHbtgGDB0t3fLkh0bGCmBh+AUjBH38A334r\nTXxcDNERJm859oYXc11HCtGRq0BU7eE1OdKlAXFsFju0ZsisWcDHH4vnndckOZnPPfUFEh0rGDyY\nXwBS8NlnwPPP8/5eYtOiBfdShN5TtiDHeo6AmBlsUnk6FF6TJ3MNeJDmbc+kLqXo9O/PG7Ju3y7+\nsQsL+XqOWCUEaoBExwrCw4HiYvE79d65A/ztb8Abb4h7XAGdjsey7Vn8llN0xEomkGr/HznCa5WV\nvMO0PXspGUOL4bVGjXgChz1tqKRYzxHQ6YCZM3kSkNjs2AH07ctDg/UFEh0r0Ol4torYIbYvvwR6\n9wYeekjc4xpib9q0VDuGGkOs8NrNm3x9SKzqewE5wmu5ubw+p3FjcY8rluiUlADZ2dJes4bYm+ot\nZo2OMZ5/nl8Tv/4q7nHr23oOQKJjNWKv61RU8ASCWbPEO6Yx7F3XkdPT6dSJp+JWVNh3HKls9vHh\nXY+l2HBOQIrQGiBeoeXZs3zjPTH2+bEEe+2WMrwG8PMwY4a43k5FBbBzJ7/RrU+Q6FhJ//68SLS4\nWJzj/fgj/2N44glxjmcKLYlO06bcXnsLAqWy2dGR12hcvmydLaVWdMCXInMNAFq2fLA/jT3IFVoT\nsNdDk1p0AOCVV4CffxbPCz58mF9nYrUsUgskOlbSogXw+OPAf/5j/7EYe9DyRmq0JDqAOMkEUoYE\nrUkmYIzfrc6cafnxpchcA8Tbn0auJAIBLYhO8+ZAQgKwbJk4x0tOrn+hNYBExybEymLbt49nlA0Z\nYv+xzGGv6MhRGGqIGOs6UgqlNes627fz73njRsvX1aQKrwHirOvIlS4toAXRAYBp04ANG8TZdG7b\ntvqVKi1AomMDwrqOvXn5S5fyjDWpe3gB3EW3cu+6asjt6YiRwSalUFqawcYYsGAB8P77wOjRlsf8\npQqvAeKIjpbCa5WVXOylyl4zpG1b3j3g88/tO86VK/xfjx7i2KUmSHRs4JFH+MKhPZPi6dPA8ePA\nqFHi2VUXHTrYl+ar1fCalKJjiafz889A9tVkrN4RjSN/RCHpx2hs+N68myxVeA2wX3SuXeNp/vZu\nE2AN9tgsVSagKWbO5HV39iSabN8OREfzcGh9g0THBnQ67u3YE2L7+GPefcBVpp1ufXz4nb81i9mG\nyC06Dz3EW+HcuWP7MdQQXpvxVjLuBUzH7oBdOPjIPlSM3YUpy6cjebfpi6ekhC/2S3Vnbq/oCOs5\ncnSnELDHZikF3BjBwXzdd+1a249RX9dzABIdmxk82PbU6exsYOtWYOJEcW2qC0dHfmdqS6fs8nI+\n+cuxVbeAszNPybWnu67Unk5GRt0h1tRU4NLNRPwRmV7t+Rv90/HpuiST78vM5LuoShV2tVd05A6t\nAfaLjlTrY6Z4801+Y2lL2n9ZGbB3LzBokPh2qQESHRuJigJ++41vrWwtS5cCL7/M01flIHl3MqLH\nRuNa8yi8+GZ0nXfZxigs5LbKsfZkiL0hNilFx92dbz9R14LxggWAf5DxtuTp2aZdTqnvzMXydORE\naIVjCztTknHgcjSixkQheqz1178tPPUU91S//db696am8iaiYm7epyY0v4mbUjRuDPTpA+zaBbzw\nguXvu3oV+PprXlwnB8m7kzF9xXSkh6cDAcB1ANNX8DtvS/d2lzu0JmBvBpvUdgshNmNp2YcO8XZJ\nD0U1wu9G3puT6YqbN41v/iX1nbmnJ/DLL7a//9QpYORI8eyxBFuFMnl3MpLPT0fxoHQI0dB0K69/\nW/nLX/j6TlycdTds9a3BZ03I07EDW9Z1Pv0UePFFeTJpACBxYyIXHAPSw9OR9I3p8E5NlBQde5I1\npLa7rgy2BQuAt94CZrw4DYEnAqu9Fng8EE89OhVJJr4CKTPXAPs8Hcb4dyK3p9O0Kf9sa9f4Ejcm\noniQfde/rURH8zXbLVuse199bH1jCHk6djB4MDB3ruVbFxcVAatW8U2Z5KKMGQ/vlFZanlGglOjY\nE167c4dPUlJuMGYqg+3YMeDkSeCHHwBXV37LmvRNEkorS+Hq6IqpU6Yi0C8GffoA06fzokJDLl3i\nXrRU2CM6WVl8Txq5Qz863QO7rflOxbj+bUWn497OokXAsGGWJV5cuMC7nYSFSW6eYpCnYwft2/Ot\nCI4etWx8YiK/+Pz9pbXLkEY6490uXR0tT5uTuzBUwM+Pp53aEssXuhFImWFlKoNt4ULeZULITIwZ\nEIMdq3cg5asU7Fi9AzEDYvDoo7x7sLF6DqnDa/asjyiRRCBgi1g6Mfuvf3sYOpRnjO7cadl4wcuR\nMzNQbkh07MTSLLZbt4AVK4C335beJkOmjTQe3pkaN9XiYyjl6eh03NuxJcQmh83GwmtpaXw9JyHB\n/Pv/8hfeMsUwZMQYP6YcomNLcbPcnQgMsUV0nuszDc6b7bv+7cHBAZg9m3s7llDf13MAEh27sXRd\n54svgAED5GsFLxAzIAbLJy9HdFY0IjMi4bgmGgtGL7dqEVUp0QFsD7HJJTo1PZ1Fi4DXX+f7v5gj\nJAR48kng739/8FxBAfeQjCUYiIU9+9MokbkmYIvo+HnFIMT1wfUfnRWN5VOsu/7t5fnngfx881vd\nFxcDBw/ypsL1GVrTsZNevYD0dH5ReXsbH1NSwu9o9+yR1zaBmAEx+j+ybt2Ah9pb934lRSckhHsP\n1iJHSNDfnxewVlTw9Olz53h9xerVlh/j3XeBP/2J12y5uspXUyJM4B4e1r3v9Gkuqkpgi+hkZAA9\nwmPw+efKuQ9OTjzCsXAhz3Y1xc8/87Y3Ndf46hua9XRu376NYcOGITQ0FMOHD0exkb0GsrOz8fTT\nT6Nz586IiorCV199Jbodzs7cg6lrq9pVq/jWBUrdIRpiy1bLSouOWj0dFxd+o5GdzR8vXswTA5o1\ns/wY4eH8RkAQKqkz1wRsmcDv3QPOn+f7HSmBLTbL3Y3AFKNGAf/7H3DkiOkx9T1rTUCzorNgwQL0\n6tULaWlp6NmzJxYuXFhrjLOzM5YtW4YzZ87g+++/x9tvv41z9pS4m6CudZ2yMl4M+pe/iP6xNmHL\nrpdKh9dOn7Z+/UEumwURv3iR33hMmWL9MebMAT78kHd+kGuStGUCP3+eJ3fI1cOsJrZs5KZENwJj\nuLjwLgWm1nYYI9FRPVu3bkV8fDwAID4+Hps3b641xtvbG2H3cw9bt26N7t27I9ee/v4meOYZHjq7\nd6/2a+vW8TvDiAjRP9YmLO2ObIicW1XXpGVLHm7IyrLufXKJjiDiH3wATJ5s21rM44/zfl1r18of\nXrMGJTPXANu2rJbLc7SEl1/mno4xz/3UKS5Mjzwiv11yo1nRuXr1Kry8vAAAXl5euHr1ap3jL168\niDNnzqBnz56i2+LlBQQF1a7yrqgAlizhcXu1YG14rbKSt/pRsiWHLSE2OT2dvXuBzZt5aM1W5szh\nwnX+vHo9HSUz1wDrbWZMPeE1gHuIr7/Ow7A1ERp81udUaQFVJxIMGDAA+fn5tZ5fVMNH1el00NXx\nbRUXFyMuLg7Lli1D0zoqy+bNm6f/OSoqClFRURbbKuyxY/iWTZt4d+ennrL4MJJjbXitqIjvlurs\nLJ1N5hAy2P70J8vfI6fozJnDF4rt6aXXuzev+0pJkU90rN1f6dQp+bbiMIa1olNQwL0HKTMBrWXC\nBP79nj/PG9oKbNumnhC8OVJSUpCSkmL7AZhGeeSRR1heXh5jjLHc3Fz2yCOPGB1XXl7OBg4cyJYt\nW1bn8ew9FYcPM9ap04PHlZX88c6ddh1WdEpLGXNxYezePcvGnznD2MMPS2uTOb76irGRI617T2Ag\nY7//Lo09hhw6xFiTJoz98Yf9x/r5Z8acnBgrL7f/WOZYu5axF1+07j0dOzL2v/9JY48lFBUx1qKF\n5eMPHWLssceks8dW5s9nbOzYB4+vXWOseXPGSkqUs8kerJ07NRtei42Nxdr7G1asXbsWw4YNqzWG\nMYaXX34ZnTt3xowZMyS1JyKCL3IKWwds3sxrIQYMkPRjraZRI95FwdK7XKW6ERhiS62OXOtQjz/O\nW96IcY6efpqnXcvhVVq7PlJczHffDAw0P1Yq3N15+UGZ8c42tVBLEkFNpk7l/dguX+aPd+0CIiOV\nS9CQG82Kzpw5c3Dw4EGEhobi8OHDePf+wklubi5i7pf0/vLLL1i/fj3+85//IDw8HOHh4dixY4ck\n9jg48P0vhG2sFy7kazlqjNFak0ygZOaaQKdOvCeVsUQNY5SV8clJjv1/dDq+nqe2Y5nD2kywM2eA\nRx/lNSdKodNZ18JHTUkEhnh4AK++Cnz0EX9cnzdsM4aq13Tqonnz5kYz1tq1a4fk+y0Cevfujaqq\nKtlsionh2WodOvAJcsgQ2T7aKoRkgqefNj9WDaLTuDEvxDx/nu8zYo7CQp74IPf+P1rC2vURpTPX\nBAS7fXzMj710CejeXXqbbOG117iIz54N7NjBb1IbCvRnKSIDB/JWF3Pn8otJrZOeNckEahAdwLoQ\nmxpCgmpHmLwtrX9Ssv2NIdZ4aGoNrwE8xD1qFN/mxMuLJ5E0FFQ6LWoTd3deYX7jBu+3pFa0Fl4D\n+F328eOWjVWLzWqmaVMerrJ0f5rjx9Xl6ViCWsNrArNmAQcO1P8GnzUh0RGZN9/k3aQdHZW2xDTW\nejpKFYYa8vzzwFdf8QVtc5DoWIalE/iJE7zjQmSk9DaZw9IEiHv3eF88ObcRsRZfX2D5cmDcOKUt\nkRcSHZGJiVFfxlpNtOjpdOrE16C++ML8WLXYrHYsFZ2FC/nNlBqyqyy1OTub787r4iK9TfYwcWLD\n6EJgCIlOA8TbG7h9W3tew7vvAp98wjPT6kJNNqsZSybwU6d4COjVV+WxyRyWio6aOhEQ1SHRaYDo\ndEBAgGXejpomcGH/mZUr6x6nJpvVjCUT+MKFwBtvWLY/kBxYIzpqTSJo6JDoNFAsCbExpr5MsDlz\neNfuu3dNjyHRsQxzE/jZs7wtz8SJsplkFktFR+1JBA0ZEp0GiiXJBLdu8Q4GrvJsJ28RYWG8C8Cq\nVabHqCX5Qe2YW5QXdkGto12h7JCno31IdBoolng6avUYhP1nSkuNv65Wu9VGXRP4778Du3cDkybJ\na5M5yNPRPiQ6DRRLtjhQ6+T92GPc4/nyS+Ovq9VutVFXoaWwC6ratk5u1YrXwVVW1j2OEgnUC4lO\nA8WS8JqaJ29h/5mazR8rK/mkpOT+P1rBlNdw8SLvIWjLLqhS4+jIi7CLikyPuXWLr/m1aSOfXYTl\nkOg0UDp04B2x62qDoua1kR49eO3O/UbjeoqK+P4pSjam1AqmRGfxYi44atqHxhBza1EZGfz6VmOz\nXYJEp8HSvDlfIK5rw1U1ezoA8N57fII07D6tdpvVhDHRycgAtm4Fpk1TxiZLMLeuQ6E1dUOi04Ax\nF2JT+wTeqxffCuDrrx88p3ab1YSbG0/GMAxRfvABT5H28FDOLnNYIjqUuaZeSHQaMOYy2LQwgb/3\nHk/trajgj7Vgs1oQ9qcRJvDLl4EffgAk3u/QbsyJDmWuqRsSHRFJ3p2M6LHRiBoTheix0Ujenay0\nSXViztNRW2GoMfr0Afz8gI0b+WMSHeswnMA//JC3u1F7EkZ98XS0Nl+IBS23ikTy7mRMXzEd6eHp\n+ufSV/CfYwaos3d5x468r5YptDKBv/ceMGEC35tEKzarBWECv3IF+OYbXp+jdjw9eYadKbTg6Whx\nvhAL8nREInFjYrULCADSw9OR9E2SQhaZx1ytjlYm8Kef5umxmzZpx2a1IIjORx8BL7+sjXNXl6dT\nVcWzMtXu6WhxvhAL8nREooyVGX2+tNJE2bwK0HoigYBOx3drnT6d727Zs6fSFmmH1q2BtDRgwwbe\na00L1FXUmpcHtGihrtY9xtDifCEW5OmIRCNdI6PPuzqqqHFZDfz8eMp0eXnt10pKeKFls2by22UL\n/fvzbKyfftKGUKoFT08gKQmIj+fbJmuBujwdLYTWAG3OF2JBoiMS00ZOQ+CJwGrPBR4PxNS4qQpZ\nZB4nJ8DHh2ct1UQoDNVKgZ1Ox9d27t4l0bEGT08ekpo1S2lLLKeu4lCtJBFocb4QCwqviYSw+Jf0\nTRJKK0vh6uiKqVOmqn5RUEibfuih6s9rJbRmyKBBwGuv8dodwjJ69QI+/pjvsqkVhPAaY7VvirTi\n6Wh1vhADHWN1NUJpOOh0OjTEU/Hqq0C3bjz7y5Dt24Fly4Bdu5SxiyDqokULICuL92EzJD4eiIwE\nxo1Txq6GiLVzp2bDa7dv38awYcMQGhqK4cOHo7iOvZcrKysRHh6OIUOGyGihNjCVTKBFT4doOJha\n19FKeK0ho1nRWbBgAXr16oW0tDT07NkTCxcuNDl2+fLl6NSpE3RaWaCQEVNdCbRQGEo0XEyJjlbC\na8/gOkcAAAv+SURBVA0ZzYrO1q1bER8fDwCIj4/H5s2bjY67cuUKtm3bhldeeaVBhs/MQZ4OoUWM\niU5pKX/O11cZmwjL0KzoXL16FV73czy9vLxw1US75Ndeew1Lly6Fg4Nmf1VJMVUgSqJDqBljtTqZ\nmbwMwNFREZMIC1F19tqAAQOQn59f6/lFixZVe6zT6YyGzn766Se0adMG4eHhSElJMft58+bN0/8c\nFRWFqKgoa03WHK1a8Xqc69erdxYm0SHUjKltGSi0Jj0pKSkWzaemULXo7N692+RrXl5eyM/Ph7e3\nN/Ly8tDGyDaBBw4cwNatW7Ft2zaUlpbi1q1bGD16NNatW2f0mIai01DQ6XiILSODRIfQDp6evPuA\nIZREIA81b8jnz59v1fs1G3OKjY3F2vvbRq5duxbDhg2rNWbx4sXIzs5GRkYGvvnmG/Tt29ek4DRk\njCUTkOgQaoY8He2iWdGZM2cODh48iNDQUBw+fBjvvvsuACA3NxcxMcYLrCh7zTjG1nXUvFU1QRjr\nSkA7hmoDKg69T0MtDgWAzz4DzpwBvviCPy4vB5o04f9T/gWhRo4cASZNAn799cFzYWHAqlVARIRy\ndjVEGkxxKCEeNcNrhYU8wYAEh1ArNcNrjFF4TSvQtELUqtWh9RxC7dQUnaIi/r9hMgyhTkh0CAQE\n8D5WlZX8MXUjINRO06bcu7lzhz8WvBxatlU/JDoEGjfm4bTcXP6YPB1C7eh01QtEKYlAO5DoEACq\nh9hIdAgtYBhioxod7UCiQwConkxAokNoAUPRoSQC7UCiQwAgT4fQHuTpaBMSHQJAbU+HCkMJtWNY\nIEqejnYg0SEAVO9KQJ4OoQUET6eiAsjOBtq3V9oiwhJIdAgAFF4jtIcgOleu8J9dXZW2iLAEVXeZ\nJuSjXTu+vUFJiXZFJ3l3MhI3JqKMlaGRrhGmjZyGmAHG+/AR2j9fguhQaE1bkOgQAHjLm/btgfR0\nLj6tWiltkXUk707G9BXTkR6ern8ufQX/WUsTqVzUh/Ml1OlQEoG2oPAaoadjR+D4caBFC8DZWWlr\nrCNxY2K1CRQA0sPTkfRNkkIWqZv6cL4ET4cKQ7UFiQ6hp2NH4PBhbYbWyliZ0edLK0tltkQb1Ifz\nReE1bUKiQ+jp0IG3jNei6DTSNTL6vKsjrS4boz6cL3d3vgb5++8UXtMSJDqEno4dgZMntSk600ZO\nQ+CJwGrPBR4PxNS4qQpZpG7qw/nS6fjaY1oaeTpaghIJCD0dO/KaBy0WhgqL30nfJKG0shSujq6Y\nOmWqZhbF5aa+nK/WrYEbNwBvb6UtISyFdg69T0PeOVTg5k0esnjnHWDxYqWtIQjz9O0L5OUB584p\nbUnDxdq5kzwdQo+bG9CypTbDa0TDxNOTb81BaAcSHaIaHTqQ6BDawdOTrletQaJDVGPSJOCJJ5S2\ngiAsY+hQ7dWUNXRoTec+tKZDEARhPdbOnZQyTRAEQciGZkXn9u3bGDZsGEJDQzF8+HAUFxcbHXfn\nzh3Ex8cjPDwcnTp1wqFDh2S2VB5SUlKUNsFmtGw7QPYrDdmvLTQrOgsWLECvXr2QlpaGnj17YuHC\nhUbHTZo0CZGRkThx4gTS0tIQHBwss6XyoOULV8u2A2S/0pD92kKzorN161bEx8cDAOLj47F58+Za\nY27evIn9+/dj3LhxAAAnJye4ubnJaidBEATxAM2KztWrV+Hl5QUA8PLywtWrV2uNycjIgKenJ8aM\nGYMuXbogISEBd+/eldtUgiAI4j6qzl4bMGAA8vPzaz2/aNEixMfH4/r16/rnWrZsiaKiomrjfv31\nVzz++OPYsmUL+vfvj/Hjx6N///4YPXp0rWPqdDrxfwGCIIgGQL3pSLB7926Tr3l5eSE/Px/e3t7I\ny8tDmzZtao3x9fVFq1atMGTIEADAiBEjsG7dOqOio2LtJQiCqDdoNrwWGxuLtWvXAgDWrl2LYcOG\n1Rrj7e2NoKAgHD58GFVVVUhOTkb//v3lNpUgCIK4j6rDa3Vx+/ZtjBo1CpcuXUJgYCC+/vprNGvW\nDLm5uUhISEBycjIA4Pz58xg9ejQKCwsREhKC9evXo2nTpgpbTxAE0UBhDZx9+/ax8PBwFhISwhIT\nE5U2x2rat2/PQkJCWFhYGOvevbvS5phl7NixrE2bNqxLly76527dusWGDh3KQkJC2LBhw9jt27cV\ntNA0xmyfO3cu8/HxYWFhYSwsLIxt375dQQvrJisri0VFRbFOnTqxyMhI9uWXXzLGtHP+Tdmvle/g\n7t277PHHH2ddu3ZlPXr0YJ9++iljTDvn35T91p7/Bi06FRUVLDAwkGVkZLDy8nLWtWtXdvbsWaXN\nsoqAgAB27do1pc2wmNTUVHb8+PFqE/esWbPYhx9+yBhjbMmSJeytt95Syrw6MWb7vHnz2CeffKKg\nVZaTl5fHTpw4wRhjrKCggHl5ebGzZ89q5vybsl9L38GdO3cYY4yVlpayzp07s/Pnz2vm/DNm3H5r\nz79m13TE4MiRIwgKCkJAQACcnZ0RFxeHLVu2KG2W1TANRUifeuopeHh4VHvOkporNWDMdkA759/b\n2xthYWEAgNatW6N79+7IycnRzPk3ZT+gne+gSZMmAIDi4mJUVlaiUaNGmjn/QHX7Kyoq0KgR3/bc\nmvPfoEUnJycHfn5++se+vr76i1gr6HQ69O3bF+Hh4fjHP/6htDk2YUnNlZpJSkpCp06d8PLLL+PG\njRtKm2MRFy9exJkzZ9CzZ09Nnn/B/ifut0TXyndQVVWFrl27wsvLC5MnT4a/v7+mzr+h/VOmTIG/\nvz8A685/gxad+lCb88svv+DkyZPYuHEjFi9ejP379yttkl3odDpNfS8TJ05ERkYGDh48CEdHR7zx\nxhtKm2SW4uJixMXFYdmyZWjWrFm117Rw/g3tb9q0qaa+AwcHB5w8eRIXL17E559/jhMnTlR7Xe3n\n35j91p7/Bi06Pj4+yM7O1j/Ozs6Gr6+vghZZT9u2bQEAwcHBGD58OI4cOaKwRdYj1FwBMFlzpVba\ntGkDnU4HNzc3TJ48WfXn/969e3juuefw0ksvYejQoQC0df6N2a+17wAAAgICMHjwYOzbt09T51/A\n0H5rz3+DFp2IiAhcuHABmZmZKC8vx6ZNmxAbG6u0WRZTUlKC27dvAwAKCgqwbds2hISEKGyV9VhS\nc6VW8vLyAAAVFRXYuHGjqs8/Ywwvv/wyOnfujBkzZuif18r5N2W/Vr6DwsJCfejp2rVr2L59O0JC\nQjRz/k3ZLwimxedf/PwGbZGSksLCwsJYly5d2PLly5U2xyouXbrEunbtyrp27cr69u3L/va3vylt\nklni4uJY27ZtmYuLC/P19WVr1qzRTMqoYLuzszPz9fVlq1evZqNGjWIhISHsscceY6+99hrLz89X\n2kyT7N+/n+l0Ota1a9dq6a1aOf/G7N+2bZtmvoO0tDQWHh7OQkND2cCBA9mqVasYY9pJmTZlv7Xn\nX7PFoQRBEIT2aNDhNYIgCEJeSHQIgiAI2SDRIQiCIGSDRIcgCIKQDRIdgiAIQjZIdAiCIAjZINEh\nCJVz4MABzbc3IggBEh2CUDnPPvssVq1apbQZBCEKJDoEoWLOnTuHP/74A08++aTSphCEKJDoEISK\n+e9//wuA7+VDEPUBEh2CUDH79+9Hy5YtERwcrLQpBCEKJDoEoTK2bduG7t27o3v37tiwYQMcHR31\nj8+ePau0eQRhF9TwkyBUSnZ2Ntq3b4+//vWvmDZtmtLmEIQokKdDECpl7969AIA+ffoobAlBiAeJ\nDkGolJSUFLi5uSEsLExpUwhCNEh0CEKlpKSkUKo0Ue8g0SEIFZKVlYXMzExKlSbqHSQ6BKFCUlJS\nAFRfz1m2bJlC1hCEeJDoEIQKOXr0qD5VGgAOHjwIZ2dnha0iCPtxUtoAgiBqExYWBnd3dzg5OeHW\nrVtYuXIl1qxZo7RZBGE3VKdDECrk7t27mDx5MnQ6HZydnfHee++hXbt2SptFEHZDokMQBEHIBq3p\nEARBELJBokMQBEHIBokOQRAEIRskOgRBEIRskOgQBEEQskGiQxAEQcgGiQ5BEAQhGyQ6BEEQhGz8\nP+VYaAsqXCuqAAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 360
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that we have gone through all that trouble to create the even-samples matrix, we can finally put it into the framework of the $L_1$ minimization problem:\n",
      "\n",
      "$ min_F || \\mathbf{F} ||_{L_1}$\n",
      "\n",
      "subject to\n",
      "\n",
      "$ \\mathbf{U}_{N/2}^R \\mathbf{Q}_r  \\mathbf{F}= \\mathbf{f} $"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rearrange_G( x ): \n",
      "    'setup to put inequalities matrix with first 1/2 of elements as main variables'\n",
      "    n=x.shape[0]\n",
      "    return hstack([x[:,arange(0,n,2)+1], x[:,arange(0,n,2)]])\n",
      "\n",
      "K=2 # components\n",
      "Nf=128 # number of samples\n",
      "M = 18 # > K log(N); num of measurements\n",
      "\n",
      "# setup signal DFT as F\n",
      "F=zeros((Nf,1))  \n",
      "F[1]=1\n",
      "F[2]=0.5\n",
      "F[Nf-1]=1 # symmetric parts\n",
      "F[Nf-2]=0.5\n",
      "ftime = dftmatrix(Nf).H*F # this gives the time-domain signal\n",
      "ftime = ftime.real # it's real anyway\n",
      "time_samples=[0, 2, 4, 12, 14, 16, 18, 24, 34, 36, 38, 40, 44, 46, 52, 56, 54,62]\n",
      "half_indexed_time_samples = (array(time_samples)/2).astype(int)\n",
      "Phi = dftmatrix(Nf/2).real*Q_rmatrix(Nf)\n",
      "Phi_i = Phi[half_indexed_time_samples,:]\n",
      "\n",
      "# inequalities matrix with \n",
      "G=matrx(rearrange_G(scipy.linalg.block_diag(*[matrix([[-1,-1],[1,-1.0]]),]*Nf) ))\n",
      "# objective function row-matrix\n",
      "c=matrx(hstack([zeros(Nf),ones(Nf)]))\n",
      "# RHS for inequalities\n",
      "h = matrx([0.0,]*(Nf*2),(Nf*2,1),'d') \n",
      "# equality constraint matrix\n",
      "A = matrx(hstack([Phi_i,Phi_i*0]))\n",
      "# RHS for equality constraints \n",
      "b=matrx(ftime[time_samples])\n",
      "\n",
      "sol = solvers.lp(c, G, h,A,b)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "     pcost       dcost       gap    pres   dres   k/t\n",
        " 0:  0.0000e+00 -0.0000e+00  4e+02  2e+01  3e+00  1e+00\n",
        " 1: -1.5648e+01 -1.2218e+01  2e+03  2e+01  3e+00  4e+00\n",
        " 2: -2.3184e+03 -1.7022e+03  1e+06  8e+01  1e+01  6e+02\n",
        " 3: -2.2814e+05 -1.6566e+05  1e+08  8e+01  1e+01  6e+04"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        " 4: -2.2818e+07 -1.6568e+07  1e+10  8e+01  1e+01  6e+06\n",
        " 5: -2.2818e+09 -1.6568e+09  1e+12  8e+01  1e+01  6e+08"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Certificate of dual infeasibility found.\n"
       ]
      }
     ],
     "prompt_number": 361
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import itertools as it\n",
      "\n",
      "def dftmatrix(N=8): \n",
      "    'compute inverse DFT matrices'\n",
      "    n = arange(N)\n",
      "    U=matrix( exp(1j*2*pi/N*n*n[:,None] ))/sqrt(N)\n",
      "    return matrix(U)\n",
      "\n",
      "M = 3\n",
      "\n",
      "np.random.seed(5489) # set random seed for reproducibility\n",
      "Psi= dftmatrix(128)\n",
      "Phi= randn(M,128)\n",
      "s=zeros((128,1))\n",
      "s[0]=1\n",
      "s[10]=1\n",
      "\n",
      "Theta = Phi*Psi\n",
      "y = Theta*s\n",
      "\n",
      "for i in it.combinations(range(128),2):\n",
      "   sstar=zeros((128,1))\n",
      "   sstar[array(i)]=1\n",
      "   if np.allclose(Theta*sstar,y):\n",
      "      break\n",
      "else:\n",
      "   print 'no solution'\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%qtconsole"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}