1_Sparsity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sparsity\n",
    "In this notebook we will take a look at sparse signals and note some of their properties. We will follow the following notational conventions. We will define a signal, $\\mathbf{x} \\in \\mathbb{R}^N$ to be *sparse* if $K$ elements of $\\mathbf{x}$ are non-zero, while the remaining elements are exactly zero. We will define the *sparsity* of $\\mathbf{x}$ to be the ratio of the number of non-zeros to zeros,\n",
    "$$\\rho = \\frac{K}{N}.$$\n",
    "\n",
    "Below, we will take a look at a few such signals and study their properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Including Dependencies\n",
    "%matplotlib inline         \n",
    "import scipy\n",
    "import numpy as np                  # Numpy, très important\n",
    "import matplotlib.pyplot as plt     # PyPlot for Plotting\n",
    "import scipy.fftpack as fftpk       # FFTPack for Harmonic Analy.\n",
    "import scipy.linalg as la           # For Linear Algebra\n",
    "import pywt as pw                   # PyWavelets\n",
    "import scipy.signal as sig          # For signal processing, like conv\n",
    "from DepsSparsity import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating Random K-Sparse Signals\n",
    "---\n",
    "\n",
    "For now, let us confine ourselves to synthetic examples of exactly $K$-Sparse signals. Let us take a look at how to generate such a signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RandomKSparseSignal(K,N):\n",
    "    \"\"\" Random realization of a K-sparse signal.\n",
    "    \n",
    "        Assumes that all non-zero components of the realization\n",
    "        are drawn iid from a zero-mean Gaussian with a variance \n",
    "        of ``1``. \n",
    "    \n",
    "        Input\n",
    "        ---\n",
    "        K: int\n",
    "           Number of non-zeros in the signal.\n",
    "        N: int\n",
    "           Signal dimensionality (length)\n",
    "        \n",
    "        Returns\n",
    "        ---\n",
    "        x: array_like, float\n",
    "           The K-sparse signal realization, shaped ``(n,1)``\n",
    "        s: array_like, bool\n",
    "           The support of ``x``, shaped ``(n,1)``. Has value\n",
    "           ``True`` at every non-zero location of ``x`` and\n",
    "           ``False`` everywhere else.\n",
    "           \n",
    "        \n",
    "    \"\"\"\n",
    "    x = np.random.randn(N,1)                # Generate iid Gaussian samples     \n",
    "    s = np.zeros((N,1),dtype=bool)          # Cast ``0.0`` as ``False``\n",
    "    si = np.random.permutation(N)[0:K]      # Random support locations\n",
    "    s[si] = True                            # Flag on-support values\n",
    "    x[np.invert(s)] = 0.0                   # Suppress off-support values of ``x``\n",
    "    \n",
    "    return x, s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets take a look at what what one of these $K$-sparse signals looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAFICAYAAADtSTf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm45UV95/FP3aW32900DdKoiA0yIhk0TYOaxJgcWcT1\nyaMgxiWZdsLimuiMQMyYiY6PYcsYlxEFNPPMxMcIiGOMBAUaTkaYKGu7wIgDTYuooCwX6AW6+96a\nP+oU99enz7n3nPNbqn6/er+e5z73/u49S51z61TV91f1rZ+x1goAAAAA0AxjoQsAAAAAACgOQR4A\nAAAANAhBHgAAAAA0CEEeAAAAADQIQR4AAAAANAhBHgAAAAA0yEToAozKGMO1HwAAAAAkzVprun9X\n2yBPkrjGXzwuvFC6/XbpkktCl6R6H/nIR/SRj3wkdDEASdRHxIX6iFhQFxGTIuujMXvFd5JYromC\nzM66rxRt2bIldBGAp1EfERPqI2JBXURMqqiPBHkoRMpBHgAAABATgjwUIuUgb8OGDaGLADyN+oiY\nUB8RC+oiYlJFfTR1zWszxti6lr2JPvEJl5P3938fuiQAAABAGowxPTdeYSYPhUh5Jq/dbocuAvA0\n6iNiQn1ELKiLiEkV9ZEgD4Ww1n0BAAAACIvlmijEeee55Zpf+UrokgAAAABpYLkmSpXyck0AAAAg\nJtEGecaY0zpf54YuCxaWcpDHOn/EhPqImFAfEQvqImKSbE6eMeY4Sddaay+RdGjnGBFLOcgDAAAA\nYhJlTp4x5jRJstZe0pnJu6cT8GVvQ05eRD76UZeT9/Wvhy4JAAAAkIZ+OXkTIQqzkK6Abr0ktvOI\nHDN5AAAAQByiXK7pGWPWS7rVWrspdFkwv5QvocA6f8SE+oiYUB8RC+oiYlJFfYxyJi/jOGvth/r9\nccOGDVq7dq0kadWqVVq3bp1arZakuTeP42qO7723rV//WpLiKA/HHHPMMcccc8yxP/ZiKQ/HaR97\no9x/06ZNmp6eliRt2bJF/USZkydJxpjTrbUXd34+zlq7sevv5ORF5C/+wuXkXXVV6JIAAAAAaajV\ndfKMMcdLOtcYc7cx5hFJRHORIycPAAAAiEOUQZ619lpr7Wpr7WGd79eFLhPml3KQ1z31DoREfURM\nqI+IBXURMamiPkYZ5KF+Ug7yAAAAgJhEm5O3EHLy4vIf/6PLybuOOVcAAACgErXKyUP9MJMHAAAA\nxIEgD4VIOchjnT9iQn1ETKiPiAV1ETEhJw+1kXKQBwAAAMSEnDwU4j3vcTl5/+f/hC4JAAAAkAZy\n8lAqZvIAAACAOBDkoRApB3ms80dMqI+ICfURsaAuIibk5KE2rHVfAAAAAMIiJw+FOPVUl5N3662h\nSwIAAACkgZw8lCrl5ZoAAABATAjyUIiUgzzW+SMm1EfEhPqIWFAXERNy8lAbKQd5ANB0P/+59JGP\nhC4FAGBQ5OShEG9/u8vJu+OO0CUBABTtxhulD3xAuumm0CUBAGSRk4dSsbsmADQXqzUAoF4I8lCI\nlAcArPNHTKiPKMPMjPsaFvURsaAuIibk5KE2Ug7yAKDpaOMBoF7IyUMh3vQml5N3992hSwIAKNo1\n17icvB/9KHRJAABZ5OShVJzlBYDmGnW5JgAgDII8FCLlII91/ogJ9RFlGLWNpz4iFtRFxIScPNRG\nykEewmBWAagOM3kAUC8EeShEypdQaLVaoYuQnAcekI4+OnQp4kR9RBlGPZFHfUQsqIvhTU9L73hH\n6FLEoYr6SJCHQjCThyo98YT0yCOhS4FY3H679I//GLoUzcZMHoC8HnpIuu660KVIB0EeCpFykMc6\n/+ox4Owvxfr4r/8qffOboUvRbOTkoe6oi+HRd88hJw+1kXKQh+rRUSCL+lA+3mMAedGOVIsgD4VI\nOchjnX/16Cj6S7E+Uh/KNzNDTh7qjboYHm2184MfSEcd1Sr9eQjyUIiUgzxUj/qGLAYO5Zud5T0G\nkA9ttfORj0jXXFP+8xDkoRAp767JOv/q0VH0l2J9pD6Ub9SZvBTrI+JEXQyPttrZvVv6wQ/apT9P\nI4O8bduoRFVjZgVVoqNAFu1P+ZjJA5AXfbcz6kmzYTUyyDvtNOlb3wpdirSkPMhinX/16Cj6S7E+\nUh/KN+p7nGJ9RJyoi+GlPFbMmpmRDj+8VfrzNDLI27rVfaE6fHBRJQb1yKI+lI82HkBetNVOVe9D\nI4M8KlH1Uh4AsM6/enzG+0uxPlIfyjfqe5xifUScqIvh0VY7MzPSHXe0S38egjwUIuUgD9XjM44s\n6kP5qsohAdBctNXO7t3M5I2MSlS9lIM81vlXz9e1VHd0nU+K9ZE2v3yjbrySYn1EnKiL4c3MpL0b\nuzczIx12WKv05yHIQyH40KJK/vPN5xwSbX4VmMkDkBd9t0NOXg67d7svVCfFmbzZWenTn2adfwh0\nFP2lWB8J8srn2/dh2/kU6yPiRF0Mj77bmZmR7rqrXfrzNDLIo8OvXopB3rZt0llnhS5FmugokEWb\nXz7//qbWzgMoDn23w0xeDnT41UsxR8rXM9b5V4+Oor8U6yNtfvl8Gz/s+5xifUScqIvhjboioGlm\nZqRDDmmV/jwEeShEih9cvyw4pcA2FgR5yKLNLx+fOQB50Y44zORJMsasH+V+dPjVSzXIk6TrrmsH\nLUeK6Cj6SzHvhDa/fOTkoe6oi+HRdzszM9Ldd7dLf55ogzxjzPGSLhvlvnXbeOXCC6XPfCZ0KfLx\ns1kpzWqRoxJOiicV0F+KOcFVY3AGIC/aEaeqE5MT5T/FaKy11xpjNo9y37qd1X3ggXqVt5cUB93+\nRMLLXtYKWo4U0VH0l2LeSd3a/Doa9aRWivURcaIuhkff7czMSAcf3Cr9eaIN8vKoW4e/e3e9yttL\nikEejVU4vPfIqlubX0ejbrwCAB59t0NOXg516/Drtry0lxSDPP8/Y51/9ego+kuxPtatza+jUWfy\nUqyPiBN1MTz6bmdmRrr33nbpz9PImby6BU3M5NUTOXnh0FEgiyCvfMzkAcgrxbFiL8nn5A1iw4YN\nWrt2rSRp1apVWrdunVqtlmZmpM2b22q359Zg+zM4MR7v3i3dd199ytvreNs2dzw7G0d5qjjevFmS\nWnrpS1tRlCel4zvvdMczM3GUJ6bjViu9+vjLX7b1xBOSFEd5mni8ZYsktTQ7S33kmGOORzt2gU1b\nN94orV0bvjyhjnfskJ797NHvv2nTJk1PT0uStrjGuSdjI90O0RhzsqSLJZ1qrf1aj7/bfmV/9rOl\nt71NOv/8kgtZkHe/20X0F10UuiSje+5zpfvukx5+WFq9OnRpqrFpk3TUUdLPfiYddFDo0qTlwgul\n97xH+vGPpcMPD10ahHbSSdLtt6tz4gVl+MAHpE9+UvrpT6WDDw5dGgB19LnPuTHvXXdJz39+6NKE\ns//+Lk751KeKeTxjjKy1pvv3Y8U8fPGstV+11q7uFeAtpG5Ld+q2vLSXFC+h4P9nN9zQDlqOFLHk\noz9/1i8ldWvz62jUJdIp1kfEiboYHqkWzsyMW8FXtlov1+ynbh3+7t31H6ymOOj2QV6d6lpT0FEg\nq25tfh2RgwwgL/puh5y8HOo2M9aEAUqKQZ7/n734xa2g5UgRHUV/ft1+SprQhsZu1I1XUqyPiBN1\nMTz6bmdmRjrwwFbpz9PIIK9uHT67a9aTP5FQpxMKTUFHgay6tfl1xEwegLzou52qJqOizcnLo24d\nft1mHnuZnZXGx9MaAPg69t3vtoOWI0V0FP2lmHdStza/jkadyUuxPiJO1MXw6LudmRnp5z9vl/48\nzORFoCkzeakFeeTkhUNHgazZ2bTanhCYyQOQV4qrvnqZmanmPWAmr4fbbpM+9rF8ZXjoocFv24SZ\nPGuliYm0dtf0dWzdulbQcqSIIK+/FPNOijyxd+ed0i9+UcxjNQk5eag76mJ49N1zben++7dKf65G\nBnl5g6Z77pG+9718ZTjiCGnbtsFu24Qgb3bWBXkpnZ1hJi8czgYiq8gg79Oflr429IV7mo+ZPAB5\nEeRV+x40Lsgb9Wxj1q5d7iuPxx4bPHBrynLN1II8/z+7+eZ20HKkiI6ivxTzTooM8opo/5uI6+Sh\n7qiL4dF3z732Bx5ol/5cjQvyiqhAO3e6r1FZ6wYJgwY8zOTVEzN54dBRhPfJT0rf/GboUjhFBnl5\n2/+mKuIEKoC00XdXuyqicRuvFFGB8p7J9YP/Qf+BddsoppeUN1554QtbQcuRIjqK/qrKO7njDsmY\nSp5qQczklW/UgQl5UIgFdTE8+u65177vvq3Sn4uZvB7ydvL+voOWoSkzeakFeTRW4fDehxdTMDQz\n41ZQFLHxU0yvKybM5AHIi76bnLxcirhAdVFBHss1m83/z267rR20HCmio+ivqryTmJY1FlkfCPJ6\nG3UmjzwoxIK6GF6sm6ZdcYV02WXVPJdvS3/963bpz9W4IC+mmbyUgryUL6FAoFE93vvwYgqGCPLK\nx0wegLxi7btvv919VYGcvBzqGuTFVuGHlfJM3m/8RitoOVIU69nAGFSVd8JMXlrIyUPdURfDizXI\n27mzuvGEf+0rV7ZKfy6CvB5CBHl1H6ymGOT5Olb3Wdg6irWjSElMwRBBXvmYyQOQV6x9d4ggj5y8\nEcQQ5Pmz26ks1/QbHqS28Yr/n/3wh+2g5UhRrB1FDELk5D30kHTjjZU8bU9FzuzGNEMZk5kZaWyM\nnDzsbccO6frrQ5diYdTF8GLtu6ts9/1rf+SRdunP1bggL6aNV1LZXdPn4Y0yAKizWBurFPDeh5dt\nJ9tt6dxzw5WFmbzyzcxIk5N85rC3W2+VzjordClQB7H23SGCPHLyRhDDTN6wyzXrfp282VkX4KUW\n5PnA/PnPbwUtR4pi7ShiECIn76mn3FcoBHnlm52VFi0a/j0mD6r5Qn/+B0VdDC/WvnvXrurK5J9n\naqpV+nM1biaviAqUN6JPbXdNa+eCvJR21yxi1hijibWjSEk2GAo9yCPIK5+fyUvpRB4GE/rzH4q1\nbut9DC7WvjvETB45eSOo40xe3YO8VGfyfB2788520HKkaGZGWrw4vo4iBiFy8kIP8mZmXE4wQV55\nZmdHW65JHlTzhf78D6rounj//dJb3lLoQzaeXxEQ21gxRJD32GPt0p+LIK+HFDZemZkpLodmdlYy\nxn3F9sEtk/+fEWhUzw84U6pvsYltJm+UpYS9EOT1VreZvOnp0CVIx1NPpblZ0ebNrq2oy2ciBrHm\n9lYd5FU1Xm5ckFfUxiszM6MvPazDdfKeeEL6y78s5rFSn8k79NBWyGIkKdaOIgYhcvJC70hZZJAX\n+rXEapCZvAcflD7xiT1/FyoP6nd+R9qyJchTJ2fnznrM5BVdF++5x32vw2uPRZFtdZGqbPd373Yr\nkRYvbpX+XI0L8opYxuWDtFHP5tZhueaTTxb3vKkGeczkhRNrR5ESZvLSMshM3o9/LH3pS/M/jrXS\nAw8UW7ZeHntMevzx8p8H4T//oWze7L6n+NpHFWvfXfVMXlXvQSODvLxvXlFBXsyXUHjySfe9iMYp\n1SDP/3/vuqsdtBwpirWjiEGqOXkEeeUaZPb8ySfn+havuz7efLN08snFl6/bjh3uC+UL/fkfVNFt\no5/J667z6C/WvrvqIG/xYmnr1nbpz0WQ18MoQd7DD+99/9hn8rLf80g1yNu9W5qYSOs1xyLWjiIl\nzOSlZZA82EECq8cfr2aGrVfAiXI89ZTrD1PrC5nJG16sfXeImTxy8kbg17rmzcnLfl/I9u3S4YfP\nHQ+78crMjLttlQ1kkUFeqpdQ8GdjDj64FbooyYm1o4hBqtfJKypHkyCvt1Fn8rrrYxUzbNYS5FXJ\nf/ZjD3aKbhsffNB9Jqhng4s1nz5EkDc52Sr9uRoX5IXIydu6dc8zk8PM5PnbTExUW+l9o1REZ5vy\n7pps4x9GrB1FSrLBUOjNSoqqD9YS5PUz6EzeQgPeKoK8nTvnAj2Uzwd3qW1YtH27tHo19WwYXEKB\nnLxcinjzsjvGDWL79j2fb5ggzy/5m5iodskmyzXz80He3Xe3QxcluAcfrPZ/H2tHEYOQOXkhZvL9\ncxaxvb9vxwny9jboTF53ANddH6sI8oo8iYmF1WUmr+i2cccOad9943/dMYl1FU6IIG/Hjnbpz0WQ\n18MoyzVnZ+cGG6MEeUVdyHdQBHn5cUHuOSefLH33u9U9X6wdRSr8JWayOXnWVp9b7MsyNlZMGzrs\nCb6U+BMrCwV5CwX7VQZ5dZthOf30PfP76yJ7sicV1rqx36pV9atnIcW6CoecvJoIFeRJc/+wOszk\n+U6WIG90fibvoINaoYsS3PR0tRcfjrWjiEEVOXndwVDIM/kzMy7AKyLIy7uzcpMNcgkF369k60G/\nnLwyZ33rGuR985vS/feHLsXw6jKTV2TbuHOnG7ctXx7/645JrCdoiwjy7r/ffYYX4t+DsbFWvicc\nQOOCvFAbr0h7L/UZpBLv3u0GJyzXrB8/kxdi9iI227a5r6rE2lGkoruNJMhrvkGXa0rzz9T1CgSL\nVuRJzCpV3Y4WpS5BXpF27JCWLnVjgLrVs5Bi7buLCPJuuEH67GcXvh05eTkUNZO3ePHgHb3vUPxz\nDrO7ZhNy8lLdXdOfUPjpT9uhixLc1q3uqyqxdhQxqCInr8kzeZOTBHm9DLrxirRnv9IrJy/7vQx1\nzMmztvp2tCh1CfKKbBu3b5eWLZOWLIn/dcck1r67iCDviSfc10L8e7BzZzvfEw6gkUFeEbtrLls2\n/EyeD9LqsFyz6Jm8FHfXJCdvDjN5adm1y73/2d01s9+rVHSQN0zbn5KiZ/KqCPLqNMPy1FOu/6zz\nTF5Kuaw+yGMmbzixplr4IC/PRMWgJ2nIycuhqJm8UYK87uWag/wDZ2aaEeSluFzTz+QdeGArdFGC\nmp11n4Gqg7wYO4oYVJWTNzXVzJk8grzeRp3J666PVcyy1XG5ph8cMpNXniLbRr9cc8mSetWz0GI8\nQTs768ZzxuQr19atw83kWdsa/ckG1MggL9RM3ihBXuiZvKKuk5dikDcz4xr41HPyfP2vcnDCJRTC\n6m4jmxLkdQevmDPMTN58g15m8nrz7WedZ/JiD/KKlJ3JS+l15xVj3+1XpixenK/tf+KJ4WbyyMkb\nwe7dc2cbR512DRHkcQmF+vEzefff3w5dlKD8oITlmnGoKieveyZv6dL6B3m7drkTN7OzabVlgxhk\nJq/XycOQOXl1CvJ8+1nHmbydO8N9/odRVk5enepZaDH23Tt3ujK5PLnRH2eYmbyJCUlql97PNC7I\n82/e2NjolWjYIK/fxiuD7q7Jcs16IifPCbHMKMaOIiU+GPLLW556Slq5MlyQ59ufIoK8RYvYfKWX\nQT5zgyyTZOOV3uq+XDPU5z+U7dvnlmum9LrzijHVoqgg74knXJuz0GvzJyaL6LMW0sggL+9Z3WGX\n7NR1ueby5QR5efiZvP33b4UuSlChZvJi6yhiUVVOng+Gdu50g5wVK5oxkzc5SZDXi5/JW2i55uLF\newZXva6Tl/1ehh073InaOs2w1H255ooV8S9zLjonj41XhhfTCdovfEH61a+KncnLfu/H91mTky2C\nvGH5N29iorqZvLrurrlqVXGXUPC7a6Z4CYUYGquQWK6ZHh8M+R02d+4MN8grI8jL7hwKZ9CLoe+7\n78IzefvsU/5MXlH9W1XqvFwz5EmeULiEwmhi6rv/9m+lm28udiZPGjzIqyJNq7FBXp43r+qcvFAX\nQ1+oMx5UqjN5frnmL3/ZDl2UoLZudQF+1cs1mcnrbZi8kx/8QPqbvxn+OXrN5IVarjU7O9fm521/\nmMnrb9CNV7qDq145eatXE+R1q/tMXh2Waxadk8fF0IcXU9/94IPSAw8UP5O3UF6eH/db2y593N+4\nIM+/eePjowVNvuItWdL8SygU1QmmGuT5mbyUXnMv27ZJ++/PTF4d3Xij9NWvDn+/7hmv0Dl5Re6u\n6YO82JeeVW3QSyjsu+/C18mrIshbqByx2bbNlbmuM3l1CPKK5JdrsvHKcGLpu3fulB5+2AV6Rc7k\nLV48+EweOXkjyNvhj3Imt9/GK8PurlllkLdjB0FeXv4SCqtWtUIXJaht26Q1a6oN8mLchjkWw+Sd\n3Hef+xoWOXnpGXQmr3uFSK+cvLKDvEGWjcZm69bq29Gi1GW5ZpE5eSzXHE0sffevfuW+Z2fy8rb7\nW7dKz3zmwjN5vs9asiThnDxjzEnGmOOMMWcOc7/szFhVQV6vmTy/DfdCsjl5VV9CgeWa+fiZvNSv\nk+cHJ+yuWT/33Sf98pfDD1KyM3nbtrkOqymXUCDI25u17muQSygMMpNX9ixbXZdrHnggM3l1wXLN\n0cTSdz/4oPte5HLNJ54Y7DOcfE6eMWa9JFlrN0qaNsYcNeh9i5rJGybxvl+QF/slFFat4mLoefic\nvIceaocuSlDbtkkHHMByzVgMk3fys5+57z//+XDPkT3zuXWr+xzkvZDsqIoO8riEwt5mZ13e7ULv\nca8ZtFA5eXWbyQuxIqIoO3e6IC/2Jc5F5uRll2umFNzmFUvf/cAD7v9X5HLNYWfydu9upxnkSTpF\n0qOdnzdLOn7QO4ZYrtkryBs0Vyt0kMfumqNjd00nxOAklo6i7u67TzrooOGXbGZPhmWDvKbM5LG7\n5p4G3dxmkOCqquWaRZ3ErEqIFRFFSXUmj0soDC+WvvvBB6UXvai4mTxf9/fbj5y8QayS9EjmeL9B\n75h345VRg7zs86UW5KU+k7diRSt0UYLaulV6xjPc56Cq/39MO3TFZtC8k5kZt1Tzt397+CAv2yk2\nMchjJm9Pg1xwfnbW/f+7L4+QrY/WuttUsVyzbjN5dQ3yZmbc/35qKv4gj5y88GLpux94QFq3rrgg\nb+tWd+3pFSsGn8mbmio/J2+i3IfPxSx0g2OO2ft3v/iFdMYZ7oP3+te778N46im3znrZMumzn5X+\n1/9a+D533eUCppkZ6ZJL3HU3Dj5Yuu026cIL55/dmp6WfvM33W2efFJ6xSsWriBF+MlP3Nrhm27q\n/T7O58gjpde9Tjr3XHe8davbXXHRIumv/9pdYLIK550nfeYz0v3353ucc8+V/tt/G+5xzjrLBeVT\nU+7/PMx7uHy5dMEF0rvf7f7vixa5evLOd9Yjv+/MM6VvfUv64Q/d8X33SR/9qGvcjjnGDQTL5hvU\njRuHr79FOvVUafNm6brr3PEb3+he/0I7Vhrj6u6HPrTw5/2AA6Q3v9m1M7t2SddfP//tx8eliy6S\nvvhFqdWSzjmn9+1273aPffjh0n/6T9KnPz3/42b9+tfSa1/r2sr3v999Dqam3Gu64orBH6cIW7e6\ns6eLF0sf/7hrg6XB/xdZ/nUtWyb90R+5OlakM8+Uvv1td+mKMr3hDa4ejLJzai+zs3OB/OWXu36j\n122mplw78KUvSd/5zt63sda9tytXSuef7/qRUdrwRYvc/d797t6DxXvukV79aumxxwZrH37/96UX\nvMB9boZ11lnuf/r977vj5z9fevxxN3gcxubN0qc+5Xb8G6ZNW7TI9WH/+T/vHSD+3u9JL3yhdMgh\n0ic/6d7nN7zBfW5f/nJX9rzjDf9/n5qSvvxlt2Nv7BYtcuO7d75z7/oz6GfnnnvcbZcvH20c1RTD\nfna6+27/v3jXu6od/9x/v+v3du+W3vc+Nw4fG5P+9E+lv/qr4R9v504XB6xaJX3iE9I//EP/23bH\nKa98pXTEEaO1PwuJNciblrS68/O+kh7udaMDD9ygZz1rrSRp+fJVOvzwdTr66JYOP1w6/PC2HntM\nOvroliTp1lvbkgY73n9/6f/+37b+6q8Gu/2iRdJb3tLWLbdIP/pRSx/4gPv7d7/rZnnOP3/++x98\nsHTKKW3ddJN0000t/cu/DFfeUY6///221qyRbrihpd27B7//c5/b0jveIe3c2dZBB0kf/rD7+733\ntrVkifTv/l055e0+/uAH2/rGN6Srr27p29+W7rxztMf7znda2rRJ+qd/auszn5Fe8pKF7/93fyf9\n4z+29cQT0m/9Vkuf+tTcuupBnv/YY6V/+Ie2Zmeliy5q6Q1vkL7ylbZ+9Svpssuqef9GPd68uaX/\n/b+lyy9v62Mfk17+cvf3Rx5p68ILpcMPr6Y8d9zR1v77Sxs3tmRtmPfjO9+Rbr21pR/9SDr2WFcf\nvve9lsbGpGOOaeuoo/rf//3vb+tzn5PuvLOlK6+c//le/3rpqqvc8Ute0tJ73ys99VT/27/1rW19\n6UvSV77idn495hjpqKN63/7EE117d9BBw7/+Zz7Tff5/8Qvp7rvb2mcf6RvfqO79zx53tz///b+3\n9U//JB1wQEtveYu0fPngj3foodJ117X18MPFlveaa1z7/s//LP2H/+DerzLej+98R7rssraMkd7y\nlpZ+//eLefzly6WTTmpp3br+tz/22Jae8xx3/Sf/d/83f7zPPu7/deSR0h13tPStb0kXXNDWokWD\nl+fkk9v67GeliYmWPvvZ3rdfs0b6wQ9a2rZt/se7807p4x9v6667pFe+sqU3vnHw92fTppa+/33p\nyivbeu973efpxhvdtVNf//rh3+8//MOWXvxi6V/+ZbDb+/7oD/6grRe9SPr85+f+/tOfSt/6VkuP\nPy5t2tTWVVdJ/+W/tNRuS9PT7vV+//stXX11/vrx4x+3dcAB0lVXjXb/qo797173urb+5/+Udu1q\n6eKL5/6+Y0dL117r+rMXv1j6kz+Z//Fe97qWxselT37SjQFCv76qj6emWvrEJ9x4+YgjpD/7s4Xv\nPzUl/fznbZ13njs+6STpy19u69e/li69tNry//Eft/TqV0sbN7r6+7u/29JPfzr64/n+dM2ahW//\nnOe4E6w7d7b08Y+39ZOfDNf+3HXXJm3dOi1J+sUvtqgva+3IX5L2kbTO/5znsboe9yhJp3V+PtM/\nR9dtbExe/nJr221r//2/t/bii609/XRrX/EKa085ZbD7v/Wt1n7yk9buu2+55cxretraFSusfd/7\nXHlD+dPn3NLrAAAgAElEQVQ/dc+/aJG1O3aM/jh/+ZfW/vmfu8cZ1HnnWfvBD1r7rGdZe//91l5/\n/fVDPecBB1j7N3/j/ufWWvtv/o21f/d31v7O7wz1MEFcdpm1r3mNe79mZ0OXJqwvfMHad7zD2he/\n2Nrvftfab3zD2te+1trXv97ar399/vt+6EPWvvGN1q5bt/DzvPa1bm/DZz/b2qOOsvbWW+e//fr1\n19u//Vtr16yx9uyzrT3nnMFfU1P88z9b+6pXubr6zW+GLo1z8cXW/smfWLt0qbVbt5b3PFddZe2J\nJ1r76lfH8dr7tY8f+IC1//W/Wjs+bu3OncM95oknWvuWt1j79rfnL99tt1n7m79p7YYN1n7xi8Pd\n95xz3GfsOc+x9qc/zV+WUdx5p2sfrr12z9//8IfW/sZvuPfI95VXXmltq2Xtscda+573WHvYYWHK\nHIqvi+vWWXvuudaecMKef//Wt6x95Std23HlldWXr27uuMPaF7zAfRa/9KXRHuOII6y98EJrf+/3\nii1bHVx//fVPtz/veMfw7U9WJybaK54aaWFV59IGd8ttjrKx8+szjDHfHuXxullrb/fPI2naWrup\niMctk88H8ZdwGBtz07cTA86VTk66Nd6Tk+WWM69ly1w5t21zP4eSfb/Hx0d/nOXLXQLuMMuy/LXB\n/IWTh13n7/OYFi1yx8a412EWXKAc3oEHuqVma9bUo7xl6q6DflMKv0nFfA44QPrRj9z3hRx1lPt+\n5JFzOWPzWbPGzV7s2jXY7ZtoYsK9dp/zHIOVK91SvJ07y207fT5hLP/7fu2jb0d9nzkMt9rGfc9r\nctLVk927h3+/fC59yHp2xBHS177mlmZn+RzZp56aa5emplxayJNPumX2++wTpMjB+Lq4eLFbyrt4\n8Z5/7+7bMb+lS11u7Y4d7udRjI+7OppnHFdXrVar9PZ66CCvE3hdI7ekcqM6uXPW2vMlbTLG9MkA\nGY619hJr7UZr7SVFPF7Zugd8PsgbtOLWJcibnHSvbXradRihhAzy/M57o34oFy92eRB1DfLuv999\nT53fhCK7KYUfTC2Ul3jAAdLddw8W5K1f774vXz5YnfNtya5d6Q5W/OB71654grx99nGfnZUry/2s\nxxbk9TM56QaHExPDvx/77+9yVIsI8vwJgVHqit9wzW/4ForPI8vyQd6TT861U1NT7j1PNcjz+gV5\n2b7d98/ob+nSuZMGBHmjiS7Ik3SepKOttcdYa0/QnrtgXiTpjEJKVjM+6Mju7tnEmTzJnYX+9a/D\nzuSNjc0l6eYZME1NuQT5YQJW/6H0OzINe+2d7I6E0lxwUMWGJXn54G7NmrDliIGfueveDnnQIG92\ndrD38cQT3eYmgwZt09Pt5Gfy8szOlGXlSnddwpUry32eGGaXsvq1jxMTc0HesPbf3923qJm8UT8r\nExN7ruCJid/a/6mnXPn8pjdPPuneuxSDPF8XFy92m+P0mslLud0cVhEzeRMTro7G9vmpQrvdLv2E\n5CjDylV+OWW/v49amDrzZ/TyLNfcsaMeDcvUlAvyQs/k7dzp3uc8QV6e5ZqjdgR1Xq65fLkbKDCT\n13v23gd9gwR52e/zWbZMOuGEwQcfExNzM3mpDlZiC3QkF9w9+GD5A+u6DFTz9HnPeIb7XtRM3qh1\nJcZ65mWXa+7a5fqXZcvmZvIefTS9IM+bL8jL07enZsmSuSBv2J3sPWby5lILYpnJu9cYc2qfv50l\n6bYc5amt7Bm9UZZr+oFZHRqWZcukhx4Kn5O3a1f+hmHU5Zo7d859KFPKyTPGBXgEeXNBXveFogfN\nyct+X8gwA/eDDnI5eTMzrvOsQ5tStBhz8vyAuuyZvLIHDcOaLydv1D7PB3ehZ/KyyzVjqWeev36b\nD/LGx93v/PI6Kb0gz9fFJUvccs3uwCS7SieGz07sFi1ydX/r1nwzecNMiDRJFTl5o7ytZ0u6xRhz\ntlxAt9oY8zlJL5a0XtKbCixfbfQ6qz/M9Gvdlmvec0/4mTx/Efo8pqbmrnc3qDy5JFLvIG92th5B\nnkSQ5/X6zPtLaSw0k+cHp2UEeZOT7iy1VJ/VAUXzMyzWxjN48MFd2QNrH+DOzsbz2nvJu1wz+z1v\nOfLM5GXTNGLiT0bu2OHqw9jYXJDny5pakOf1m8nzOXljY+TkDcIYF9w9+ig5eaOKbrmmtfY2Sa/s\n3PdNcsszz5AL8M621lZ8Kdw4dOfkNXV3TckFRH59fyhFzuRlvw9ictLtLur/V6Pm5PlOpE45eZL0\nohe567ukLs9yzYkJdwHvYYK8QZcR/epXLidPqk+bUjSfkxfTxisrVrjvVczk+dcew/++X/vYhJk8\nP0CLMSdvbMy9nq1b5/rKbA6V5C7cnJJsTt58u2vG8tmpA4K80bXb7dJXXozULFlrr5X0PGPMoXLB\n3aOSbrHWPlZk4eqk3yUUmra7pjQX3IWeyRvm/e1nlCBv0aI9g7xhLVrklojWcbmmJH3uc6FLEIfs\n7prZ5Zr+bwt5//sHD5aHncl7pLMd1vbtaZ6R9oNvP9CNwdiYC/SqyskbH4/ntfeSJyfPB3f77VdM\nOUbdpCebix/jINUHM76vzO4eKaU9k9cvyPOz4DF/dmKyZIm7NMyoOXkpb7wixbm75tOstZuttV/t\nXOrgMUkq6hIKddPrmllN3nhFCj+TV0SQ51/LsMs1s4PnlHLyMCe7u2b3JRQGqZcf/vDc7M5ChukI\nDjmkpe3b3c95TkbUWYw5eZKbxatid81RLwlQhn7tY57lms94hvTpTxfz+vK8X36AmncDsLL43DO/\nBNGYPQfjqQV52Zy8xx/vn5PHTN7g/AweM3nDa7Va8S3XnI8x5hC5nL3kdJ/Ra/JyzRhm8vz7m3eJ\nYxHLNYdV95w8OP1y8spYejvMhaN9/ZTq06YULdZdD/fZJ73lmv3k6fPGxqT3va+4cuS5GHrMsxCL\nF88FLX4QvXSpq4Pj4+kFeZ5/X/rl5PnLI2FhRQR5RZywr6voZvKMMbPGmJleX5LuKb6I9ZBNwM6z\nXDPWziJr2bLwS4GKahhGCVi7l2umlpMHJ09O3rCyA+KFTgb8/Ofk5GUH7jG1qStXVrdcM5bdNefL\nyYth9Yr/3I4SrI2Pu41MYqpjWT6IyZ4QXbLEDcj32Se9IC+bk5f97pGTN7ylS917NepYLPYTJWVq\nt9tPrwJ68smIcvIkbZSUzb9bJek4SdOSrs1bqDrqNeCTmjmTNzXlvkLOPBW18YpPRh92Jm/79rlr\nNQ3L73rGcs1663cJBf+3Ig3TPvjbSvVpU4rmZ/Jiy62pcrlmXXbX3HffsOUwxtWRUQZZfoAa6yyE\nD2K6Z/JmZtymT2vWhCtbSPMFeeTkDWfJktHz8aTidkqvszJPeI3SBUxba1/Z6w/GmIskXZqvSPVU\nRJAXw1nNQSxbFjYfTyp2in/58tGCPP+/GjYnz3csBHn1lucSCsMaJsh7/vNbBHkTcwPbmAKdk06S\njjmm3OeIbaC60HXyBt1htkyj5gfGPgvhB9/ZIG/JEtdO3XprHPWjStmcvOx3z392ZmbSe29GtXTp\n6Es1pfg/Q2Xy9TFPfvJCRnnIU/r9wVp7hjHmZrlr5iWlX5DXxN01p6aaFeT5mclB+XX7eXLyst9Z\nrllP3btr+v+jMcX/Lxctcu3DILsJTk66TlOqT5tStBivkydJp59e/nP4parWxv2/j+nE5qhl8ZtG\nxFTHsrIzednlmrGcAAil30ze+PjcSdeU359h5A3yUs/Jk8od/49ynbyFlmM+b8Sy1Fr2OnkTE3MV\ndpiZvKeeqkfDsmxZ2E1XpPAzedJckDZKTl72OzN59ZTdXTO7XLOsnLxBB+1btrSf/jn2gX5ZYs3J\nq0Jsu2v2ax/LPHs9rDwzedmLi8cmm5OXXa6ZZ3ldnS2Ukye5tmNigv54UEUEeTEveS6Tr49RLdc0\nxpwryfb58wmSNucqUU0VsVwz+z1msczkFZGTJ7mckGEuCpv3f0WQ1wzZz3z2Egr+b0Uaps51tzl1\naFOKFutMXhV8XYzxAt1ZfmlcDPXTX1uyibtrSr2Xa6ZsoSCPVTWDy5uTF/tnqAq+v4plueZZ8/xt\nWtLRI5al1npdQsH/fhB1CvJimMkr6hIKknTFFXMX1x2ED85GzcnrFeRxCYX6qTonL/t9Pkce2drj\nOMWtwH17nOqytMnJeE4czZeTl/0e0qhl8bMQsc6MZXPyupdrpmihnDzJtZepB8HDYLnm6Hx9LLMt\nHDVuPEHSo12/m7bWJjmLJ/W+hIL//SBi6vAWsnLl4BdxLkuRM3nD7pLZvVxzWOTkNUP37ppl5uQN\n0z5036YObUrRsjOrKX6u6rDczPeNMZzBH7Usfrlm6JOe/cy3u2bKmMkrDhuv5Ffm+H+UqnyJtXaj\ntfa2rq9kAzwp/3JNf7s6DMhe8xrp858PW4aQZ3+6P5Cj5uT5DoblmvXU7xIKZeTk+Q0BBmkf/t//\na+9xXIc2pQwTE/UIdsowORnP/32+6+Rlv4c0alliH6D223gl1pnHsg2akxdDnayLIi6hkHpOXpkn\nvEbZeOWM+f5ujEns8ppOEbtrZr/HbNEi6cADw5YhZJDXvVxz1PuTk1dvfjbZmLnZu2zQV7RBBx/d\nz12HNqUMfgOFFNVhoBpTnzfqIKtOF0PPzuTlmXlpAoK84rDxSn6xzeQt5LISHjN6KW28EoMYZvJ8\nkEZOXprGxvZcBuWXB5a1RHBycrAlwuvXt/a6X4omJtJ+7bEEHv3ax5iWa+aZyYt5c59+18lLdSZv\n0Jy8FPOYR8VyzdF15+SVMZ6d9201xlym/jtp9rKvpONylaim/NmIPJdQyH7H/GII8oq8Tl4ZS/xQ\nLl8Hs7P2ZeXkSYOfYSYnz5mYSHeDiTrkFcXU5/llvcO+Z75/j3UWIrtc0/c37K5JTl6Rht2dvBsz\neeWmFiwUgpw8wmMOExQ2RlG7a6Z4NmMUIYM8Y9zzZnPyhpnNYyavGbrrYJmXUJBcfRlkQHzHHW1J\nraePYxhEh0CQF7oUTr/2MaYgb9QlesOezK1adrnmypXu56mpdIM8XxcJ8orzR38k/eEfjn7/mGb0\nq+brY5lLhAd5W/e11j42yIMZY1ZJuiVfkerJ767Jcs1qFLm75ijyLOnwHQs5efXWK8gr6xIK0uAD\nUd/mLF0a94Way+YvI5CiOmw4E9PgbtTlrTG9hl567a753ve6JaYpI8grTt6l4b5eptpPSeXmgS70\nr7lk0ABPkqy108aYs3OWqZbIyauWz4cK1RhnP5R5c/JS3uq9znyQ1/1/DL1c8yUvaUlK+4y9VI9A\npywxDVTrcp28UcpRh+Waixfv2Vfut1/YMoXEdfLik3KQ5+tjmTnU8z7sQjtpdjPGHKK9r5+XBB/k\ndV8nr4m7a8Yg9AU085x5YXfNZvCzyT7p3OfkSWGDPF+vpqbcTF6qsu1wamIK8vqJqc8bdZAV+3LN\nJUukZcvCrnqJ0eLFrr/t9X+rw2enSWKfDa9CmTN5uaqyMWZl9kvSekkXFVO0emEmr1r+/Y5hueao\n18nz/2ty8uppbMyd1Olerhn6EgqbNrUlucFdyu1JTDtMVi2m196vfYxpcJd3Ji+G19DLqlXS6tVh\nV73ExNfFZcukAw7o3edyCYVqpTyT5+tj6Jy8PXRm6y6XC+isJP8x8T8neVH07iBv2IpLkDec0A1D\n3pm87CyDn8mLdaCA3rrrYCzLNX09Sj3ISzknrw6zETH1eU1drvn2t7tA75RT4i1jCEuXSpv7jFRZ\nrlmt0GO5GEQV5Ek6T9Khkq7ofLeS7u38fIikowsrXY1kg7zsAJ6ZvHKEbhjy5uRlN20hJ6+essGd\nP44hyHvZy1qS3HLNlNsTcvJCl8KpQ05eU5drTky4doDlmk62Li5b1vs2MX12UhD7bHiZgufk9XG8\npKOttfcaY9ZLOt1a+05JMsZcJHedvK8VWMZamJhwS7f88i2CvHKFDvLy7K7ZfV9y8uqp30ye/7lo\nw14nL/WZvJRz8uoQ4MY0uGvqck1p7jOQ6mdhWAR51Qo9lotBbDl5j1hr75Uka+1tko7J/O1sSX9R\nRMHqxl8nT3INBBuvlCt0w5D9UI6Sk0eQV3+9grwYcvJuuaUtiSAvpry0qsWUV9SvffT9ZAzlzHsJ\nhZgHqKH7ypgM0lfH9NlJQR1OlJSlipy8UYK86c5GK2s7x48YY9Z1ft5XLlcvOePj0lNP7b2Ei5m8\ncoQ+O5k3J687yGPjlfrpN5NX1tJbcvKGMzmZ5sBBqs9ANZZyNvVi6BJB3rDyrNLB8Kif8S3XvFjS\ntCRrjHmTpPMlXWeMuUbSCUp445WdO+f+UcMGef52MXR4dRC6Ych2BMPm5B10kLsgrTc25uoOS0Tq\npXu2vuycvEEHH74+kpOX7omTmF77fO1jLLOto5Zj2H4+hNAnRGMySF/Ncs1qhR7LheTrY7CNV4wx\nD1tr97h0prX2YmPMKrllmhuttY8ZYy6XdFrnJqd1P04Kuq/bNuxyzfHx/tdtwd5CNwx5PpTLlkln\nnz13zExePXXXQf//K+vSHuTkDSemQKdqdRmo1n0mz/fZMQ9QQ/eVdVOXz05TpLxc0wu5XHNfY8zn\nOtfAe5q19nxr7SnW2sc6x2dIOkzSamvtFeUUNW7dQd4oyzhi6fDqIHTHlScnrxs5efXU60RO9rIY\nRRu0ffjXf21LIsibmEj39cfUl8zXPsZSzjwziuPjcQ9Qu088p4ycvPiEHsuF5OtjmSsaBvnYnyGX\nh3epMebYfjey1m621k4XV7R68Tl5oy7XlGhchhG6YShy3T4zefXkl2VmB0/ZnXWLNmj74J9/airt\n3JKUc/JiWQa5kFjKmafvjeU19DPsqqLUkZNXLWbywm+8crTcLN2tki42xjxijDkns/EK5Cpor41X\nhmlYCfIGFzrIy3OdvG5+V0bOtNbP+PjeM3mhg7xXvKKlyUly8mIffJcppr5kvvYxlnI2OcgL3VfG\nZNCcvBjqZCpSrp9V5OQtNBzZaK29vTNLd7619jC5XLwxSbcZY+42xnywezlnivrl5DGTV47QDcOp\np0q/9VvFPBYzefXVK8grq04O0z7st5+0YkXa7Unsg+8y1aUviaWceZdrxjxAZeOV4axdKz33uaFL\nkY7QY7kYBFuuaa09ocfvNltrz7bWrpb0TrlZvi3GmG8bY95YTjHjR5BXrdAd1+tfLx18sPuZnLx0\ndQd1Zc7kDbqMqN1u6yc/kZ7xDJeXl6qUc/JiCnDnax9jKSczeWkYpK8+9VTpne8svyxwUl6uWcV1\n8vK+rQ93vq+Su3zCCcYYa61Nrjnpl5PHcs1yNKnjYrlmfXWfxS+zPp51lpudG8SKFdKrXiX99m+X\nV57YTU66z1WK6rJDYCx9Xp5AjSAPGB31M+ByzV4zc8aYfYwxZxpj7pZ0i6TTO3/aLOlsSft13ycF\nReyu+ZrXSAceWHzZmiimhiFvTh7LNeurypy8Qw6R9t9/4dv5+jgx4ZZtpir2wXeZYgmepDRy8mLo\nh/oJveolJnn7ahQvprFc1bJ9daiLoX9B0tckyRhzktxOm8dJ8sPRxyRdKukia+3t5RSxHvxM3uLF\n7niU5ZoXXlh8uZqqSQ0DQV59de+myUAqHilfJ68urz2WQDyFSyg0oa9E86S8XNMLufHKqk6u3Yyk\nyyUdLxfgXSvpTdbafa217ywrwDPGrC/jccvgd9fMs1wTg4up4yInL11VzuQNKm99bIpYAogQDjww\nnlUh89XHQw6RDjigurL0k0JOXuh2KQa0jfGJaSxXNV8fjzhCOuaYcp5jkKbJb76yWdJFki72F0Ev\nkzHmeEmfl9vYJXqLF0vbt+fbeAWDa1LDQE5efVWZk4fh1CUvrQzvfW/oEgzm0ktDl8BZvtxdcmQU\ndVmuGXMZkS5m8qTjjnNfZRjkbf2qpHOqXo5prb3WGLO5yufMY2pK2rWLIK8qMQV55OSlq9fumqGR\nd+LUZcli09WhPr797dIpp4x2X5Zr1kcd6mJqUq6fVdTHBZsma+2ITV9ali933/NcDB2Da1LDQJBX\nX72WayIOBHkYVJOXa7LxCmLWpLFcjBb62NdiqWQM/FIP39j7DRkYZJQjpo6LnLx09VquGbpOknfi\nxD74TkXT62PsyzUZRM9pel2so5SXa1ZRH+d9W621pS2XNMac1uPXm621Gwd9jA0bNmjt2rWSpFWr\nVmndunVPT3/6N6+q41tuccfj4+74ttvancFemPI0/fjGG92xf79DlyfP8diYtHVrW/ffL1Ff6nU8\nPt7S+Pjc8dhYS9bGU76Uj1eskF74wnjKw3Ezj8fHpXvvbavdjqM8vcon0b9kxVIejufq5003SYce\nGr48VR57o9x/06ZNmp6eliRt2bJF/Rhrbd8/hmaMudpa+8o+f7OxlX1yUnrpS6UbbpC+/313IeLt\n20OXqpl27ZIWLZI+9jHpwx8OXZp8Tj1Vuu466c1vls45J3RpMIwXvEA68kjpq191x897nlt6e++9\nYcsFoBq/+7vS294mvetdoUvS26OPSqtXS2efLZ17bujSAHvatEk66ijp/vulZz87dGnqyxgja+1e\n68HGQhRmEMaYkyUd0+uC7LFavnzPSyikOP1clSYtQSEnr75ivIQCgOqwXBMYHfWzXNEOR6y1X7XW\nrrbWfi10WQY1NbXnxisEeeXxA+kYGobuqfdhjY1xCYW6GhsjJw+YT9PrY+y7a/r2KHS7FIOm18U6\nSjnIq6I+Rtw01U93kJdipa1S9yxKXbHxSn31msmLbBU5gBLFvsFPyoNoxC/ljVeqwNtaIGbyqhVL\nkOeTYUfFcs36ivESCnnrI1CkptdHlmvWR9PrYh2lXD+rqI+EIQWamtrzEgoEeeVqSv4TM3n11esS\nCgDSEXtfz3JNxIyZvHLxsS/Q8uVzg7zJSbf7I8oTy0weOXnpinHjFfJOEJOm18eVK6UVK0KXor+U\nZ0q6Nb0u1lHK9ZOcvJqZmprLxzn4YOnb3w5bnqaLJcjLi+Wa9dW90UroAA9Atb74RXdSN1bGuC/a\nJsQo5SCvCnzsC5RdrmmMdNhhYcvTdLEEeeTkpat7g6UYZvLIO0FMml4fFy8O/5lfCBvBOU2vi3WU\nTXFKTRX1MfKmqV6yG6+gfLEEeXmRk1dfvXLyYh/wAUhLU/pKNI/vMxn/lIPhSIGyOXkoXywdFzl5\n6SInD5gf9TG8GNqlGFAX4xPLOC4EcvJqhpm8ajWlcWC5Zn3FeAkFAMhqSl+J5lm8WFqzJnQpmosg\nr0DZnDyUL5azk+TkpSvGSyiQd4KYUB/DI8hzqIvxWbJE2rw5dCnC4Dp5NfOsZ0kPPRS6FOloSsfF\ncs366jWT53fYBYAYxHJCFOgl5t1p646PfYHe+lbp3HNDlyIdsQR5eddV+xk8ZvLqp3vwFMNgirwT\nxIT6GF4sfWVo1EXEhJw8YB5N6bgI8uqLnDwAsWPXXyBNBHmorViCvCJy8rLfUR/k5AHzoz6Gx3Xy\nHOoiYkJOHjCPWIK8vPwZVs601g85eQBi15S+EsBwGFaitmLpuMjJS1evIC90nSTvBDGhPoYXQ65w\nDKiLiAk5ecA8Ygny8iLIq68Yl2sCQFZT+koAwyHIQ23FcnaSnLx0dW9oEMNyTfJOEBPqY3gEeQ51\nETEhJw+YR1M6LnLy6qt7eWYMQR4AZMVyQhRAtfjYo7YmJ+O4iCY5eekiJw+YH/UxvKacEM2LuoiY\nkJMHzOOSS6TDDgtdivwI8urrmc+UDjhg7piBFIDYMJMHpIkgD7V1xBGhS+CQk5euj350z+MYlmuS\nd4KYUB/DYybPoS4iJuTkAQkgJ685YgjyACCLIA9IE8NKICdy8uDFMJgi7wQxoT6Gx3JNh7qImJCT\nBySAIK85mMkDEJsYTj4BqB5BHpBT3nXVLNdsjhiCPPJOEBPqY3gEeQ51ETEhJw9IADN5zUGgDiA2\nBx8s7bdf6FIAqBpDEiAncvLgjY+HD/TIO0FMqI/hXXqp9G//behShEddREzIyQMSQJDXHDEs1wQA\nACDIA3IiJw9eDEEeeSeICfURsaAuIibk5AEJYCavOdjcAAAAxIC5AyAncvLgxXA9KvJOEBPqI2JB\nXURMyMkDEkCQ1xwxLNcEAAAgyANyIicPXgxBHnkniAn1EbGgLiIm5OQBCWAmrznIyQMAADFg7gDI\niZw8eOTkAXuiPiIW1EXEhJw8IAEs12yOGJZrAgAAEOQBOeVdV81MXnPEsFyTvBPEhPqIWFAXERNy\n8oAEEOQ1x6JF0sxM6FIAAIDUsUAMyImcPHjve5/0Z38WtgzknSAm1EfEgrqImJCTBySAnLzmWLky\ndAkAAAAkY2u6S4Axxta17EDWF74gnXaa9PWvS3/wB6FLAwAAgLowxshau9d6sGhn8owxp3V+fJ61\n9s+DFgYoEcs1AQAAUKQoF4gZY46TdK219hJJh3aOgSiRk4eYkHeCmFAfEQvqImJSRX2MMsiTdKik\n4zs/b+4cA41ETh4AAACKFH1OnjHmaklnWWs3df2enDw0wv/4H9KGDdKVV0qveU3o0gAAAKAuapeT\nJ0nGmPWSbu0O8LwNGzZo7dq1kqRVq1Zp3bp1T19c0E+Dcsxx7MdumWZbP/yh9JrXhC8PxxxzzDHH\nHHPMMcdxHm/atEnT09OSpC1btqifYDN5mY1VsjZbazdmbnOmtfaCPvdnJg9RaLfbT3/4RvH3fy/9\n8R9LV10lvepVxZULacpbH4EiUR8RC+oiYlJkfYxuJq+zqUpfxpjTfYBnjDkuG/wBTTI2tud3AAAA\nII8oc/KMMcdLukzSI5JWSzrZWntd122YyUMjfPnL0tveJl19tXTCCaFLAwAAgLqIbiZvPtbaa+WC\nO6DxuIQCAAAAisQCMSAnnxQ7KpZrokh56yNQJOojYkFdREyqqI8MK4HAmMkDAABAkaLMyRsEOXlo\nigPg0qsAAAtNSURBVMsvl045Rbr+eomNvwAAADCofjl5zOQBgTGTBwAAgCIR5AE5kZOHmJB3gphQ\nHxEL6iJiQk4ekABm8gAAAFAkcvKAwL7+dekNb5BuuEF62ctClwYAAAB1QU4eEClm8gAAAFAkgjwg\nJ3LyEBPyThAT6iNiQV1ETMjJAxLATB4AAACKRE4eENiVV0qve530ve9JL3lJ6NIAAACgLsjJAyLF\nck0AAAAUiWElkFPeddUs10SRyDtBTKiPiAV1ETEhJw9IAEEeAAAAikROHhDY1VdLJ54o3XabdNRR\noUsDAACAuiAnD4gUOXkAAAAoEsNKICdy8hAT8k4QE+ojYkFdREzIyQMSQJAHAACAIpGTBwR2/fXS\nscdKP/yhdOSRoUsDAACAuiAnD4gUOXkAAAAoEsNKICdy8hAT8k4QE+ojYkFdREzIyQMSQJAHAACA\nIpGTBwR2ww3Sy18u/fjH0uGHhy4NAAAA6oKcPCBS5OQBAACgSAwrgZzIyUNMyDtBTKiPiAV1ETEh\nJw9IAEEeAAAAikROHhDYTTdJL32ptHmzdMghoUsDAACAuiAnD4gUM3kAAAAoEkEekBM5eYgJeSeI\nCfURsaAuIibk5AEJIMgDAABAkcjJAwK7/XZp/XrpZz+TDjoodGkAAABQF+TkAZFiJg8AAABFIsgD\nciInDzEh7wQxoT4iFtRFxIScPCABBHkAAAAoEjl5QGA/+pH0whdKDzwgrVkTujQAAACoC3LygEgx\nkwcAAIAiEeQBOZGTh5iQd4KYUB8RC+oiYkJOHpCAsbE9vwMAAAB5kJMHBHbXXdILXiA9/LC0enXo\n0gAAAKAuyMkDIsVyTQAAABQp2iDPGHO8MeY4Y8y5ocsCzIecPMSEvBPEhPqIWFAXEZNkc/KMMcdJ\nOslau1HSemPMutBlAspCTh4AAACKFH1OnjHmbmvtYT1+T04eGmHzZul5z5Mef1xasSJ0aQAAAFAX\ntczJM8acKemM0OUAysRyTQAAABRpInQB5mOtvcAYc5kx5hZr7WPdf9+wYYPWrl0rSVq1apXWrVun\nVqslaW6tK8ccl32cXVc9yv1dcNfWd74jvfrV4V8Px/U+zlsfOeaY+shxE4/972IpD8dpH/vfjXL/\nTZs2aXp6WpK0ZcsW9RNsuaYx5rQev95srd1ojFkvyVprb+9svPKwtfaCrvuzXBNRaLfbT3/4RnHf\nfdJznytt3y4tXVpcuZCmvPURKBL1EbGgLiImRdbHfss1o8zJ6yzTvK0T8H1e0tXW2q913YYgD43w\ns59JBx8s7dghLVkSujQAAACoi35BXqzLNS+WdIox5lC5Gb2vLXQHoK7IyQMAAECRxkIXoBdr7WPW\n2ks6X+8KXR5gPtn11aMYG9vzO5BH3voIFIn6iFhQFxGTKuojw0ogMGbyAAAAUKQoc/IGQU4emuKB\nB6RnPlPavVsaHw9dGgAAANRFLa+TB6SAmTwAAAAUiSAPyKmonDyCPBSBvBPEhPqIWFAXERNy8oAE\nMJMHAACAIpGTBwT28MPSM54hzc6GLgkAAADqhJw8IFLGMIsHAACA4hDkATkVkZPHNfJQFPJOEBPq\nI2JBXURMyMkDamDTpk257s9MHoqUtz4CRaI+IhbURcSkivpIkAfkND09nev+BHkoUt76CBSJ+ohY\nUBcRkyrqI0EeEBhBHgAAAIpEkAfktGXLllz3Hx+XJieLKQuQtz4CRaI+IhbURcSkivpY60sohC4D\nAAAAAITU6xIKtQ3yAAAAAAB7Y7kmAAAAADQIQR4A1JQx5szMzycZY45b6HcAAKBcxpj1XccD9dFF\n9tu1C/IYtCAUY8xpna9zM79jYI0gjDHHSzqh8/N6SbLWbpQ0bYw5qtfvghUWjWaMWd9p907K/I62\nEUFk6tlpPX5HfUTpOv3zZZnjgfroovvtWgV5DFoQijHmOEnXWmsvkXRop2M4SmJgjWCyCdVvlvRo\n5+fNko6XdIqk6a7fAWX4c2vtFXJtYyWDF6CXTr3a3Klnmzt1j74albLWXivX73qnaLA+utB+u1ZB\nnnq/SUAVDtVcfdvcOX6zSv6AAr0YY47qDE68fSQ9kjneT9KqHr8DCmWMOVnSzZJkrb3AWnu7Bh/Q\nAGU4r/P90E595CQYQuvVHw/6u5HVLchj0IIgrLWXdGbxJGm9pFvk6uPDmZsxsEZVVvf43V7bJwMV\nOEbSfp3ZEb/srfTBC9BLJ6i71xjziObqG/URMai8j65bkCcxkEFAnSUet3Y6Eon6iIr1mMWT3Blp\nH/j5kw/Z3+2rPU9IAEV6yLeJmbw82kZUzhizSm7W7hxJlxhjDvF/ClcqYME+upR+eyLPnQNg0ILQ\njrPWfqjzc68PqEQdRbkONcYcKnfmeXUnl+RSuRmVjXJLia+RG9T43x3S+R1QtIcl3dv5eVrSi0Xb\niHBOk/TX1trHjTGbJZ0s6iPCW6iPLqXfrttM3qVyb4TEoAUVM8acbq29oPPzcdqzPvoPKHUUpbLW\nXtHZ5MLK5eLZzCzKcZKmrbWbev0uWKHRZF/VXJu3StJNom1EWEZybaVcgEd9RKU6ucrHGGPeKD29\njHjBPrrofttYaxe+VUQ6W+JulkuovWSh2wNFyGyH+4jc2b+TrbXX9aqP1FEAKem0eY9IOsavdKBt\nRCid3NDNklbPV/eoj2i62gV5AAAAAID+6rZcEwAAAAAwD4I8AAAAAGgQgjwAAAAAaBCCPAAAAABo\nEII8AAAAAGgQgjwAQGMZY443xjxijJnt+jqngMd+NPN4ZxZRXgAAikCQBwBoLGvttdba1ZLO7vzq\nGmvtmL+eW87H3lfSCf4w7+MBAFAUgjwAQAru7XyfLvhxHy348QAAyI0gDwAAAAAahCAPAAAAABqE\nIA8AkCRjzFmZjVM+b4w52Rhza2ejlsuMMfv0uM/x2dtIWt3nsdcbY64xxtxtjLnFGHNc5m8XZZ73\nEWPMScaYMzMbudxc4ssGACSAIA8AkCRr7fmSju4cniK3OcuZne8nS9qYvb0x5nhJV0u6W9JaSZdK\nurz7cTu3u0XSpdbawySdJukaH+hZa8/Q3EYwN1trr7DWXiDpZkkXWWtfXODLBAAkaCJ0AQAACMhv\nyGIlHWetfVzSdcaYEySdbIxZ2fmdJF0kyVpr39w5vqJzu9O7HvMiSfdYa78gd4fbjTG3dX5/WOd3\nF3Tue4Ix5jS5DVystfZdJb1OAEBCCPIAAJA2Z4I5aS74Wy3pcWPMKkmHSLqn6363ZQ+MMYd2btc9\nw/eopPXGmH2stY91fvemzvNcJGmzpPW5XwUAAGK5JgAAklteOZ9DO983d/2++9gHaid08vHuNsbc\nLbcs9FG5AFCS1An2TukcPtoVZAIAMDKCPAAAFtYdzHmruo79dfgutdYelvlaba3dz1q7qev2Z8vN\nBh5tjDmzyAIDANJFkAcAwAKstdNyAVz3pigndB37GcG9Nk/p7Oa5T/ZY0t3W2mPkgsjzjDGHdN8P\nAIBhEeQBAFLgl1s+r+v3/hIIpuv3foZu38zvTpO0ys+4GWPWS/KXRvAbqkzLzc6t71yW4dDO13ly\nG7s81rnvmyT9eWajFb9s85pel24AAGAYBHkAgMbqXNfuEUnnyu2gub5zLbpzjDEnyV0SwUo63Rjz\n7c59LpcL6Kyky40xp0qStfYKuc1Szug85jmSzug81en++nadyyG8SdIxcpdbuEXSrLX2xMzjXypp\nH2PMpZ37n9d5vkMk3WuMOba0NwUA0HjGWhu6DAAAAACAgjCTBwAAAAANQpAHAAAAAA1CkAcAAAAA\nDUKQBwAAAAANQpAHAAAAAA1CkAcAAAAADUKQBwAAAAANQpAHAAAAAA1CkAcAAAAADfL/Ab6umY4C\n7oOaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106077b10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 1024         # Desired signal dimensionality\n",
    "K = 64           # Number of non-zeros\n",
    "xsparse,xsparse_support = RandomKSparseSignal(K,N)\n",
    "ShowSignal(xsparse);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the signal is mostly a constant value of `0` with a $K$ non-zero values randomly distributed amongst it. Often in practice, it is not quite so clear that our signal of interest is $K$-sparse or even compressible. So, another way of looking at the signal is by its sorted coefficient magnitudes, or its *deacay*. Specifically, lets call the signal $\\tilde{x}$ the *sorted* veresion of $x$. Then, we have\n",
    "$$|\\tilde{x}_1| \\leq |\\tilde{x}_2| \\leq \\cdots \\leq |\\tilde{x}_N|$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowDecay(xsparse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall from the lecture that we defined a *compressible* signal as one which follows a *power-law decay*. Specifically, if we sort $x$ by the magnitude of its coefficients and call this signal $\\tilde{x}$, the original signal $x$ is said to obey a power-law decay if\n",
    "\n",
    "$$\\left|\\tilde{x}_s\\right| \\leq C s^{-\\alpha} \\quad \\text{for} \\quad s = 1,2,\\cdots$$\n",
    "\n",
    "Can we say that this $K$-sparse signal obeys a strong power-law decay? Is it essentially *compressible*?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowDecay(xsparse,scalemode='log');\n",
    "## Can we find the proper set of parameters \n",
    "## which bound the decay?\n",
    "C = 5.5         \n",
    "a = 0.5\n",
    "s = np.arange(1,N)\n",
    "plt.plot(s,C*np.power(s,-a),linewidth=2);\n",
    "## Uncomment these lines to see the decay comparison\n",
    "## in the linear domain.\n",
    "# plt.xscale('linear')\n",
    "# plt.yscale('linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets compare our $K$-Sparse signal to one which is merely compressible, but not truly sparce. In this case, lets look at an \n",
    "*iid* Laplacian signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RandomLaplaceSignal(scale,N):\n",
    "    \"\"\" Random realization of an iid Laplace signal. \n",
    "    \n",
    "        Input\n",
    "        ---\n",
    "        K: float\n",
    "           Scale parameter for the Laplacian distribution.\n",
    "        N: int\n",
    "           Signal dimensionality (length)\n",
    "        \n",
    "        Returns\n",
    "        ---\n",
    "        x: array_like, float\n",
    "           The iid Laplace signal realization, shaped ``(n,1)``                \n",
    "    \"\"\"\n",
    "    x = np.random.laplace(0.0,scale=scale,size=(N,1))    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xlaplace = RandomLaplaceSignal(np.sqrt(0.25),N)\n",
    "ShowSignal(xlaplace);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually, we can see that this Laplacian *iid* signal has values which are concentrated around 0, but with a few large outliers. Again, lets take a look at the decay profile of signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ShowDecay(xlaplace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, likewise, we can take a look at the power-law decay which this signal obeys."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowDecay(xlaplace,scalemode='log');\n",
    "## Can we find the proper set of parameters \n",
    "## which bound the decay?\n",
    "C = 4       \n",
    "a = 0.25\n",
    "s = np.arange(1,N)\n",
    "plt.plot(s,C*np.power(s,-a),linewidth=2);\n",
    "## Uncomment these lines to see the decay comparison\n",
    "## in the linear domain.\n",
    "# plt.xscale('linear')\n",
    "# plt.yscale('linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we compare the two signals, and from our previous power-law decay fitting, we can see that the $K$-sparse signal is indeed more compressible than a Laplacian signal, one which is by its essence \"compressible\". Below, we can see the direct comparison between the two decays of the two different kinds of signals, the $K$-Sparse signal, and the compressible one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "CompareDecayTwo('K-Sparse',xsparse,'Laplacian',xlaplace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stark difference is also apaprent when we look at the distribution of coefficient magnitudes in the following histogram. We belabor this point so that we can get a feel for the what it is to observe a truly \"sparse\" signal in the wild, as compared to a \"compressible\" one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "CompareHistograms('K-Sparse',xsparse,'Laplacian',xlaplace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the nice properties of compressible signals, however, is that they are never too far away from a $K$-sparse represenation. Because the power-law property of compressible signals forces the majority of information to be concentrated on a few coefficients, by retaining only the top $K$ coefficients, we are able to create a $K$-sparse approximation of a compressible signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def KSparseApproximation(x,K):\n",
    "    \"\"\" Retain the top K coefficients in  x.\n",
    "    \n",
    "    Given a signal, find the top K-magnitude coefficients.\n",
    "    Return a signal which retains these identified coefficients\n",
    "    and sets the rest to 0.\n",
    "    \n",
    "    Input\n",
    "    ---\n",
    "    x: array_like, float\n",
    "       The signal that we wish to approximate\n",
    "       \n",
    "    K: int\n",
    "       Number of coefficients to retain.\n",
    "    \n",
    "    Returns\n",
    "    ---\n",
    "    xT: array_like, float\n",
    "        The K-Sparse approximation of x\n",
    "    S:  array_like, int\n",
    "        An array of the support index locations in x.    \n",
    "    \"\"\"\n",
    "    S = np.argsort(abs(x),axis=0)[::-1][:K]\n",
    "    xT = np.zeros_like(x)\n",
    "    xT[S] = x[S]\n",
    "    \n",
    "    return xT, S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "KApprox = 50    # Number of non-zeros to retain\n",
    "xlaplace_approx = KSparseApproximation(xlaplace,KApprox)[0]\n",
    "\n",
    "ShowSignal(xlaplace)\n",
    "plt.title(\"Original Laplacian Signal\",fontsize=20)\n",
    "ShowSignal(xlaplace_approx)\n",
    "plt.title(\"$K$-Sparse Approximation\",fontsize=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When making an approximation, we should have a way of measuring the quality of that approximation. Below we define a number of metrics that we can use to compare two signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MeanSquareError(x,y):\n",
    "    \"\"\" Compute MSE between two signals.\n",
    "    \n",
    "    Returns the mean square error, which \n",
    "    is calculated by taking the squared \n",
    "    $\\ell_2$ norm of the distortion between\n",
    "    ``x`` and ``y`` and dividing by the \n",
    "    signal length.\n",
    "    \n",
    "    Inputs\n",
    "    ---\n",
    "    x: array_like, float\n",
    "       A signal.\n",
    "    y: array_like, float\n",
    "       A signal.\n",
    "       \n",
    "    Returns\n",
    "    ---\n",
    "    mseval: float\n",
    "            The Mean Square Error value.\n",
    "    \n",
    "    \"\"\"\n",
    "    mseval = np.mean(np.power(x-y,2))    \n",
    "    return mseval\n",
    "\n",
    "def Correlation(x,y):\n",
    "    \"\"\" Compute the correlation coefficient between\n",
    "        two signals.\n",
    "        \n",
    "        Inputs\n",
    "        ---\n",
    "        x: array_like, float\n",
    "           A signal.\n",
    "        y: array_like, float\n",
    "           A signal.\n",
    "\n",
    "        Returns\n",
    "        ---\n",
    "        corrcoef: float\n",
    "                  The correlation between the two vectors.\n",
    "    \"\"\"\n",
    "    n = np.size(x)\n",
    "    x = np.reshape(x,(1,n))\n",
    "    y = np.reshape(y,(1,n))\n",
    "    return np.corrcoef(x,y)[1,0]\n",
    "\n",
    "def SignalToNoise(x,y):\n",
    "    \"\"\" Compute the SNR between ``x``, the original signal,\n",
    "        and the noisy version, ``y``.\n",
    "\n",
    "        Inputs\n",
    "        ---\n",
    "        x: array_like, float\n",
    "           The original signal.\n",
    "        y: array_like, float\n",
    "           The corrupted signal\n",
    "\n",
    "        Returns\n",
    "        ---\n",
    "        snr: float\n",
    "             The signal-to-noise ratio in decibels.\n",
    "    \"\"\"\n",
    "    Noise = y - x\n",
    "    return 10*np.log10(np.var(x)/np.var(Noise))\n",
    "    \n",
    "def EnergyRatio(x,y):\n",
    "    \"\"\" Compute the energy (squared $\\ell_2$ norm) ratio\n",
    "        between two signals.\n",
    "\n",
    "        Inputs\n",
    "        ---\n",
    "        x: array_like, float\n",
    "           A signal.\n",
    "        y: array_like, float\n",
    "           A signal.\n",
    "\n",
    "        Returns\n",
    "        ---\n",
    "        er: float\n",
    "            The ratio of the signal ``x`` energy to that of \n",
    "            signal ``y``.\n",
    "    \"\"\"\n",
    "    n = np.size(x)\n",
    "    x = np.reshape(x,(n,1))\n",
    "    y = np.reshape(y,(n,1))\n",
    "    \n",
    "    xeng = np.linalg.norm(x,2)\n",
    "    xeng = xeng*xeng\n",
    "    \n",
    "    yeng = np.linalg.norm(y,2)\n",
    "    yeng = yeng*yeng\n",
    "    \n",
    "    return xeng/yeng"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We below, we can see how well our approximation from above has done in terms of the metrics we've defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(MeanSquareError(xlaplace,xlaplace_approx))\n",
    "print(Correlation(xlaplace,xlaplace_approx))    # Between [-1,1]\n",
    "print(SignalToNoise(xlaplace,xlaplace_approx))  # Measured in dB\n",
    "print(EnergyRatio(xlaplace_approx,xlaplace))    # Value > 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of particular note is the last term, the energy ratio between the approximation and the orignal signal. We term this measure the *energy-loss* of the approximation. It can be thought of as a measure of how much information about the original signal we lost through the process of our approximation.\n",
    "\n",
    "To illustrate the point that compressible signals are amenable to $K$-sparse representation, but that not all signals are, lets compare our compressible Laplacian signal to a realization of an *iid* Gaussian signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RandomGaussSignal(sig,N):\n",
    "    \"\"\" Random realization of an iid Gaussian signal. \n",
    "    \n",
    "        Input\n",
    "        ---\n",
    "        sig: float\n",
    "             Standard deviation of the generating Normal\n",
    "             distribution.\n",
    "        N: int\n",
    "           Signal dimensionality (length)\n",
    "        \n",
    "        Returns\n",
    "        ---\n",
    "        x: array_like, float\n",
    "           The iid Gaussian signal realization, shaped ``(n,1)``                \n",
    "    \"\"\"\n",
    "    x = sig*np.random.randn(N,1)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xgauss = RandomGaussSignal(1.0,N)\n",
    "ShowSignal(xgauss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is this signal compressible? It is the essense of high-entropy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowDecay(xgauss,scalemode='log');\n",
    "## Can we find the proper set of parameters \n",
    "## which bound the decay?\n",
    "C = 4       \n",
    "a = 0.1\n",
    "s = np.arange(1,N)\n",
    "plt.plot(s,C*np.power(s,-a),linewidth=2);\n",
    "## Uncomment these lines to see the decay comparison\n",
    "## in the linear domain.\n",
    "# plt.xscale('linear')\n",
    "# plt.yscale('linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the above comparison to the fitted power-law decay, we see that the Gaussian signal does not exhibit the properties of a power-law decay. To give a comparison, what happens if we take take a $K$-sparse approximation of this Gaussian signal, keeping the same number of non-zeros as in the case of the Laplacian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xgauss_approx = KSparseApproximation(xgauss,KApprox)[0]\n",
    "\n",
    "ShowSignal(xgauss)\n",
    "plt.title(\"Original Gaussian Signal\",fontsize=20)\n",
    "ShowSignal(xgauss_approx)\n",
    "plt.title(\"$K$-Sparse Approximation\",fontsize=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, using the same distortion measures as before..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(MeanSquareError(xgauss,xgauss_approx))\n",
    "print(Correlation(xgauss,xgauss_approx))    # Between [-1,1]\n",
    "print(SignalToNoise(xgauss,xgauss_approx))  # Measured in dB\n",
    "print(EnergyRatio(xgauss_approx,xgauss))    # Value > 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we can see that since the Gaussian signal has its magnitudes more evenly distributed, less energy/information is stored within the $K$ largest coefficients, and the $K$-sparse approximation of the Gaussian signal is not as close as is the case for the compressible Laplacian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gauss_ER_TestedK  = np.zeros((N-1,1))\n",
    "laplace_ER_TestedK  = np.zeros((N-1,1))\n",
    "\n",
    "gauss_SNR_TestedK  = np.zeros((N-1,1))\n",
    "laplace_SNR_TestedK  = np.zeros((N-1,1))\n",
    "\n",
    "# Loop over all the K-Sparse approximations\n",
    "TestedK = np.arange(1,N)\n",
    "for k in TestedK:\n",
    "    laplaceapprox_k = KSparseApproximation(xlaplace,k)[0]\n",
    "    gaussapprox_k = KSparseApproximation(xgauss,k)[0]\n",
    "    \n",
    "    laplace_ER_TestedK[k-1] = EnergyRatio(laplaceapprox_k,xlaplace)\n",
    "    gauss_ER_TestedK[k-1] = EnergyRatio(gaussapprox_k,xgauss)\n",
    "\n",
    "    \n",
    "CompareER('Laplace',laplace_ER_TestedK,'Gaussian',gauss_ER_TestedK,TestedK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, here, we can see very easily that compressible signals are much more amenable to $K$-sparse approximation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Sparsity in Other Bases\n",
    "It is a rare case that we find a signal of interest which is sparse in the direct (ambient) domain. Of course, such examples exist, for example, PALM and STORM acquired images for fluorescence microscopy can be seen to be sparse in direct space, as well as some astronomical images.\n",
    "\n",
    "Most often, we have a signal which is essential non-sparse, but does contain some pattern or regulariaty. Finding a *sparse basis* often amounts to exploiting the frequency content of the signal (i.e. harmonic analysis) to find a compact representation of the information contained in the signal. \n",
    "\n",
    "Lets take a look at an example from image processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy.ndimage as img      # An image processing library\n",
    "\n",
    "## Load the test image\n",
    "TestImage = img.imread(\"elephants.jpg\")\n",
    "plt.imshow(TestImage);\n",
    "plt.title(\"A Test Image\",fontsize=18);\n",
    "\n",
    "## Converting the image to grayscale, floating point, and normalized\n",
    "X = np.array(TestImage,dtype=float)\n",
    "X = X[:,:,0]/3 + X[:,:,1]/3 + X[:,:,2]/3\n",
    "X = X / la.norm(X)\n",
    "\n",
    "plt.matshow(X);\n",
    "plt.gray();\n",
    "plt.title(\"Grayscale Version\",fontsize=18);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Here we have a nice test image. Is this image sparse? Is it compressible? It is not immediately clear from looking at it in the image domain. However, intuitively, we can see that there should be a concise description of this image. Very broadly, the image consists of just 3 different textures, the sky, the ground, and the elephants. Within these textures the textures are more or less consistent. So, the information cotent of the image, beyond describing these particular textures, is the description of the boundaries, or edges, between the textures.\n",
    "\n",
    "Note that, above, we have converted the image to grayscale. This makes the example simpler. Otherwise, we would need to consider three different color channels in all of our operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discrete Cosine Transform\n",
    "---\n",
    "\n",
    "One of the commonly used transforms in image compression (ala the JPEG compression standard) is the *Discrete Cosine Transform*. One can think of the DCT as a real version of the Fourier Transform. The DCT comes in a few different flavors, but \"the\" DCT is known as the Type 2 DCT, or DCT-II. The DCT-II of a signal is given as\n",
    "\n",
    "$$X_j = \\sum_{i=0}^{N-1} x_i \\cos \\left[\\frac{\\pi}{N}(i+0.5) j\\right],$$\n",
    "\n",
    "where the coefficients $X_j$ refer to our DCT-II coefficients. Like the Fourier transform, the DCT-II builds up an exact representation of $x$ via waveforms of varying frequency, here, cosines. Also like the Fourier transform, we can see that the DCT-II is a *linear* transform, that is, all coefficients are defined according to a sum over the time-domain coefficients and vice-versa. Finally, the DCT-II is also an *orthonormal* transform, so its transpose serves as its own inverse,\n",
    "\n",
    "$$\\Psi_{\\text DCT}^T \\Psi_{\\text DCT} = I$$.\n",
    "\n",
    "Lets take a look at the DCT-II coefficients of this test image. This requires us to apply the DCT-II in both the horizontal and vertical directions (the same procedure which is used to compute the 2D FFT)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Xdct = fftpk.dct(fftpk.dct(X.T).T)  # Apply to both directions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Out of curiosity, what do these coefficients look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.matshow(Xdct);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The scale of the values of the DCT-II is very widely varying, so to *see* the coefficients visually, we usually plot the magnitudes on a log scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.matshow(np.log10(np.abs(Xdct)));\n",
    "plt.jet()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cursorially, we see that the largest magnitudes are bunching towards the top-left corner of the 2D set of DCT-II coefficients. This region corresponds to the lowest-frequency terms in our DCT-II transform space. The implication is that the low-pass portion of the DCT-II space contains the majority of the information about the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.matshow(np.log10(np.abs(Xdct[0:50,0:50])));\n",
    "plt.jet()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even zooming further onto the low-pass coefficeints, we see that the information is further clustered into the lower $10\\times 7$ set of frequencies. What happens if we retain these but none of the others?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Xdct_lopass = np.zeros(np.shape(Xdct))\n",
    "Xdct_lopass[0:6,0:9] = Xdct[0:6,0:9]\n",
    "Xrec =  fftpk.idct(fftpk.idct(Xdct_lopass.T,norm='ortho').T,norm='ortho')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.matshow(Xrec)\n",
    "plt.gray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we can see that these coefficeints contain information about the darker regions (the elephants), the ground, and the sky. However, this is not enough information to be useful. What about if we keep the $50\\times 50$ region?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Xdct_lopass = np.zeros(np.shape(Xdct),dtype=float)\n",
    "Xdct_lopass[0:49,0:49] = Xdct[0:49,0:49]\n",
    "# Xdct_lopass = np.array(Xdct)\n",
    "plt.matshow(np.log10(np.abs(Xdct_lopass)))\n",
    "plt.jet()\n",
    "\n",
    "Xrec =  fftpk.idct(fftpk.idct(Xdct_lopass.T,norm='ortho').T,norm='ortho')\n",
    "plt.matshow(Xrec)\n",
    "plt.gray();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we have a recognizable image. How many coefficeints did we retain, here? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print \"%f%% Coefficients Retained\" % (100*(50.0*50.0)/np.size(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we can see that even with a very small percentage of the coefficients, we still retain much of the information. The image is still recognizable as three elephants in a line on a field with a horizon in the background. Adding more frequency components will allow us to get closer and closer to the original image. \n",
    "\n",
    "Just how compressible is this image? We can look at the weight decay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowDecay(Xdct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does it fit a power-law?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowDecay(Xdct,scalemode='log');\n",
    "## Can we find the proper set of parameters \n",
    "## which bound the decay?\n",
    "C = 1100       \n",
    "a = 0.9\n",
    "s = np.arange(1,np.size(X))\n",
    "plt.plot(s,C*np.power(s,-a),linewidth=2);\n",
    "## Uncomment these lines to see the decay comparison\n",
    "## in the linear domain.\n",
    "# plt.xscale('linear')\n",
    "# plt.yscale('linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, in fact, we see that it does indeed obey a power-law decay quite well, hence the use of DCT-II for $K$-sparse representation of images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wavelet Transform\n",
    "---\n",
    "Yet another tool which is used for image compression are *wavelets* (ala JPEG2000). The study of wavelets is a huge topic and is somewhat beyond the scope of this notebook. The best way to think about wavelets is that this representation exists half way between the time- and frequency-basis. Recall that when we use a harmonic analysis such as the Fourier transform or the DCT, we seek to represent our signal as a superposition of many periodic (*infinitely long*) sinusoids. These sinusoids are completely non-local. As such, a single sinusoid will tell us some information about every time sample, and vice-versa. \n",
    "\n",
    "In the case of wavelets, (né «onde-lettes»), a signal is represented by *localized* sinusoids. These localized shapes can exist at any location **and** at any scale. Because of this, wavlet representations of a signal are naturally *multi-scale*, capturing the structure of the signal across many different levels of detail. \n",
    "\n",
    "One can represent a signal *exactly* with a wavelet representation, just as is the case with the FFT and DCT-II, and this is accomplished by clever self-referential properties of the wavelet construction. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Levels = 5      # Number of levels of wavelet decomposition\n",
    "W = pw.wavedec2(X,'db1',level=Levels,mode='sym');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets take a look at the wavelet coefficients. Generally, this is done by arranging the coefficients in the same spatial orientation as the original image, concentrating the lower-frequency information towards the top left corner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ShowWaveletCoeffs(W,Magnitude=False,cmapname='jet');\n",
    "plt.title(\"Wavelet Coefficients\",fontsize=18);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we suppress the baseband and measure the *magnitude* of the wavelet coefficients so we can get a better view of detail (non-baseband) coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowWaveletCoeffs(W,SuppressBaseband=True,Magnitude=True,cmapname='PuRd')\n",
    "plt.title(\"Wavelet Coefficients (sans Baseband)\",fontsize=18);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the extreme left-hand corner we have the coefficients which are known as the *baseband*. These coefficients correspond to the idea of the DC coeffcient of a frequency transform such as the FFT or DCT. An image is *decomposed* with wavlets in a layer-by-layer fashion. At each layer, three channels of coefficients, **horizontal**, **vertical**, and **vertical**, are calculated, along with a baseband. Each successive level of wavelet decomposition performs the same calculation on the baseband of the previous level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also see that the wavelet coefficients are themselves quite compressible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowWaveletDecay(W)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets see how well we can approximate this image using a $K$ sparse representation of the wavelet coefficients. As comparison, lets take the same $K$ value as our earlier DCT experiment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "K = 50*50     # Match the same number of coeffs as DCT\n",
    "WK = WaveletTopKApproximation(W,K)\n",
    "print \"%f%% Coefficients Retained\" % (100.0*(float(K))/np.size(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we view the thresholded wavelet coefficients. Unfortunately, the scale of the plot generated by `PyPlot` doesn't give us a very good view."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowWaveletCoeffs(WK,SuppressBaseband=False,Magnitude=False,cmapname='jet');\n",
    "plt.title(\"Thresholded Wavelet Coefficients\",fontsize=18);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, lets see what the recovered spatial-domain image looks like when we drop so many wavelet coefficients from the representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "XRec = pw.waverec2(WK,'db1')\n",
    "\n",
    "plt.matshow(XRec);\n",
    "plt.gray();\n",
    "plt.title(\"Recovered From Few Wavelet Coeffs\",fontsize=18);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}