Add a fourth notebook to demonstrate the use of FAµST API to write algorithms like OMP.

gen_doc/Use_of_Faust_objects_in_algorithms.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Using the FAµST API in Algorithms\n",
+    "\n",
+    "After the little tour we've done in the previous notebooks, about the creation of Faust objects, their manipulation and their use in factorization algorithms, we shall see in this fourth notebook how the FAµST API can be deployed seamlessly in algorithms.\n",
+    "Our example, already alluded in the second notebook, we'll be the Orthogonal Matching Pursuit algorithm (OMP).\n",
+    "\n",
+    "This algorithm intervenes in the dictionary learning problem. Of course, I won't treat the theory behind but I assume the reader is already familiar to.\n",
+    "They're not so much to say so let's go straight in the code example."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 1. The OMP Algorithm Implementation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pyfaust import *\n",
+    "from numpy import zeros, copy, argmax\n",
+    "\n",
+    "def omp(y, D, niter):\n",
+    "    m,n = D.shape\n",
+    "    x = zeros((n,1))\n",
+    "    supp = []\n",
+    "    res = copy(y)\n",
+    "    i = 0\n",
+    "    K = min(*D.shape)\n",
+    "    while (len(supp) < K and i < niter):\n",
+    "        j = argmax(abs(D.T*res))\n",
+    "        supp += [j] \n",
+    "        x[supp,:] = pinv(D[:,supp])*y\n",
+    "        res = y-D[:,supp]*x[supp,:]\n",
+    "        i += 1\n",
+    "    return x"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The most important point to notice in this code is that except the import part in the header, all the code seems to be a natural numpy implementation of OMP.  \n",
+    "This is in fact the core philosophy of the FAµST API, as explained in previous notebooks and also in the API documentation, we have wanted a Faust to be seen as a numpy array (or rather as a numpy.matrix) hence this code is in fact totally compatible with the two APIs: the D function argument, which is the dictionary, can be indifferentially a pyfaust.Faust object or a numpy.matrix object.  \n",
+    "\n",
+    "In the next, we'll test this implementation in the both cases. But first, let's define a test case."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 2. The Test Case Dictionary"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For convenience, we shall set up a dictionary which garanties uniqueness of representations.\n",
+    "The dictionary will hence be the concatenation of an identity matrix and a Hadamard matrix, and because we work with Faust objects, this concatenation will be a Faust object.\n",
+    "\n",
+    "Below is the block matrix of our dictionary:  \n",
+    "$\n",
+    "D =\n",
+    "\\left[\n",
+    "\\begin{array}{c|c}\n",
+    "I_n & H_n \\\\\n",
+    "\\end{array}\n",
+    "\\right]\n",
+    "$\n",
+    "\n",
+    "$I_n$ is the identity or dirac matrix and $H_n$ the hadamard matrix, with n being a power of two.\n",
+    "\n",
+    "The condition on which the uniqueness of the representation $x$ of a vector $y$  is insured is defined by the following inequality:  \n",
+    "$ \\| x \\|_0 < (1 + 1/\\mu)/2 $ where $\\mu$ denotes the coherence of the dictionary and in case of our specially crafted dictionary $\\mu = {1 \\over \\sqrt n}$\n",
+    "\n",
+    "So let's construct the Faust of D, compute y for a unique x and test our OMP implementation to find out if we effectively retrieve the unique x as we should according to this theorem.\n",
+    "\n",
+    "Note that, for a better view and understanding you might consult this document: [Rémi Gribonval HDR](http://videos.rennes.inria.fr/hdr/gribonval/HDR2007.pdf)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[ 1.  0.  0. ...,  1.  1.  1.]\n",
+      " [ 0.  1.  0. ..., -1.  1. -1.]\n",
+      " [ 0.  0.  1. ...,  1. -1. -1.]\n",
+      " ..., \n",
+      " [ 0.  0.  0. ...,  1. -1.  1.]\n",
+      " [ 0.  0.  0. ..., -1.  1.  1.]\n",
+      " [ 0.  0.  0. ...,  1.  1. -1.]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "from pyfaust import FaustFactory as FF\n",
+    "from numpy import log2\n",
+    "n = 128\n",
+    "log2n = log2(n)\n",
+    "FD = FF.eye(n).concatenate(FF.wht(log2n), axis=1)\n",
+    "D = FD.todense()\n",
+    "print(D)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now that we have our dictionary both defined as a Faust and as a matrix, let's construct our sparsity satisfiying vector x, we'll call it $x_0$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "l0 norm of x0:  1\n"
+     ]
+    }
+   ],
+   "source": [
+    "from numpy import zeros, count_nonzero\n",
+    "from numpy.random import randn\n",
+    "from random import randint\n",
+    "from math import floor, sqrt\n",
+    "x0 = zeros((FD.shape[1], 1))\n",
+    "nnz = max(1,randint(1,floor(.5*(1+sqrt(n)))))\n",
+    "nonzero_inds = []\n",
+    "while len(nonzero_inds) < nnz:\n",
+    "    nonzero_inds += [ randint(0,n) for i in range(0,nnz) ]\n",
+    "    nonzero_inds = list(set(nonzero_inds)) # set imposes id. uniqueness\n",
+    "# we got nz indices, now build the vector x0\n",
+    "x0[nonzero_inds,0] = [ randn()*50 for i in range(0,len(nonzero_inds)) ]\n",
+    "print(\"l0 norm of x0: \", count_nonzero(x0))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It remains to compute $y$, normalizing FD and D before."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "FD = FD.normalize()\n",
+    "D = FD.todense()\n",
+    "y = FD*x0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Our test case is complete, we are fully prepared to run our OMP algorithm using a well-defined dictionary as a Faust or as numpy array we should retrieve our $x_0$ from the vector y. Let's try!"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 3. Running the Algorithm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "We succeeded to retrieve x0, OMP works!\n"
+     ]
+    }
+   ],
+   "source": [
+    "x = omp(y, FD, nnz)\n",
+    "from numpy import allclose\n",
+    "assert(allclose(x,x0))\n",
+    "print(\"We succeeded to retrieve x0, OMP works!\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We tested OMP on a Faust, go ahead and verify what I was aiming to in the first part of the notebook; is this OMP implementation really working as well on a Faust and on a numpy.matrix?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "We succeeded to retrieve x0, OMP works!\n"
+     ]
+    }
+   ],
+   "source": [
+    "x = omp(y, D, nnz)\n",
+    "assert(allclose(x,x0))\n",
+    "print(\"We succeeded to retrieve x0, OMP works!\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can conclude that the algorithm is indeed available to both numpy and Faust worlds, and we can imagine surely that other algorithms are reachable through the FAµST API. That's anyway in that purpose that the FAµST library will be extended in the future (starting from the latest version which is 2.5.1)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "***The fourth notebook is ending here***, I hope you'll be interested to try yourself writing another algorithm with the FAµST API and maybe discover any current limitation. Don't hesitate to contact us in that case, we'll appreciate any feedback!"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}