pyLDPC-Presentation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\"> Projet Stats'App: Implémentation de LDPC </h1>\n",
    "\n",
    "\n",
    "## update:  02/24/16 - v.0.7.0\n",
    "\n",
    "L'acronyme LDPC désigne les \"low-density parity-check codes\". Ces codes sont en fait des cas particuliers de codes correcteurs, qui sont nécessaires pour assurer une bonne transmission des données lors de transfert pouvant induire certaines perturbations. \n",
    "\n",
    "La partie codage consiste en la transformation du message composé de k bits en un mot codé sur n bits. On a donc n-k bits supplémentaires (dits de redondance), présents afin d'améliorer le décodage du code. Dans la suite, on notera G la matrice de codage. L'utilisation d'une telle matrice pour le codage vient du fait que le code est linéaire.  \n",
    "\n",
    "> \n",
    "Le code est un sous-espace vectoriel de dimension k dit code. Une suite de bits x appartient au code si et seulement x s'écrit v.G, où v est le message de taille (1,k). Ou encore, si x n'appartient n'a pas de composante dans l'orthogonal du code, et donc en notant H une matrice formée par les lignes (e1,e2, ..en-k) (avec (e1, ..) base de l'orthogonal) , on a H.x' = 0. \n",
    "> \n",
    "La particularité des codes LDPC réside dans la construction de cette matrice H où le nombre de 1 doit croitre linéairement avec la taille du code. Le nombre de 1 dans la matrice H est donc dans ce cas relativement faible et la matrice est dite creuse. Trouver une telle matrice à partir d'un code n'est en général pas évident. On commmencera donc par construire une telle matrice de parité puis on déterminera la matrice génératrice G qui lui correspond. \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On importe tous les modules dont on aura besoin: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from scipy import sparse \n",
    "from scipy.sparse import csr_matrix\n",
    "import cv2\n",
    "from time import time\n",
    "from IPython.display import display, HTML\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## 1- Construction des matrices de décodage / codage\n",
    " ### 1- 1 Construction de la matrice de parité H\n",
    " \n",
    " #### Matrice de parité régulière \n",
    " \n",
    "On présente ci-dessous la méthode de construction de Gallager d'une matrice de parité régulière avec la taille du message k, la taille du mot codé n, le nombre de 1 par ligne (resp. colonne)  d_c  (resp. d_v).\n",
    "\n",
    " Cette matrice doit donc vérifier: (n-k).d_c = n.d_v \n",
    "\n",
    "  > <h3 align=\"center\">  Méthode de Gallager </h3>\n",
    "  \n",
    "  > Dans l'image ci-dessous, K représente le nombre de 1 par ligne d_c, et w_c (c pour column) représente le nombre de 1 par colonne d_v.\n",
    "  \n",
    "  > <img src=\"Images/CallagerH.png\">\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RegularH(n,d_v,d_c):\n",
    "\n",
    "    \"\"\" ------------------------------------------------------------------------------\n",
    "\n",
    "    Builds a regular Parity-Check Matrix H (n,d_v,d_c) following Gallager's algorithm : \n",
    "\n",
    "    ----------------------------------------------------------------------------------\n",
    "\n",
    "    Paramaeters:\n",
    "\n",
    "     n: Number of columns (Same as number of coding bits)\n",
    "     d_v: number of ones per column (number of parity-check equations including a certain variable) \n",
    "     d_c: number of ones per row (number of variables participating in a certain parity-check equation);  \n",
    "\n",
    "    ----------------------------------------------------------------------------------\n",
    "\n",
    "    Errors: \n",
    "\n",
    "     The number of ones in the matrix is the same no matter how we calculate it (rows or columns), therefore, if m is \n",
    "     the number of rows in the matrix: \n",
    "\n",
    "     m*d_c = n*d_v with m < n (because H is a decoding matrix) => Parameters must verify:\n",
    "\n",
    "\n",
    "     0 - all integer parameters\n",
    "     1 - d_v < d_v\n",
    "     2 - d_c divides n \n",
    "\n",
    "    ---------------------------------------------------------------------------------------\n",
    "\n",
    "     Returns: 2D-array (shape = (m,n))\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    if  n%d_c:\n",
    "        raise ValueError('d_c must divide n. Help(RegularH) for more info.')\n",
    "\n",
    "    if d_c <= d_v: \n",
    "        raise ValueError('d_c must be greater than d_v. Help(RegularH) for more info.')\n",
    "\n",
    "    m = (n*d_v)// d_c\n",
    "    a=m//d_v\n",
    "    Set=np.zeros((a,n),dtype=int)  \n",
    "    \n",
    "\n",
    "    # Filling the first set with consecutive ones in each row of the set \n",
    "\n",
    "    for i in range(a):     \n",
    "        for j in range(i*d_c,(i+1)*d_c): \n",
    "            Set[i,j]=1\n",
    "\n",
    "    #Create list of Sets and append the first reference set\n",
    "    Sets=[]\n",
    "    Sets.append(Set.tolist())\n",
    "\n",
    "    #Create remaining sets by permutations of the first set's columns: \n",
    "    i=1\n",
    "    for i in range(1,d_v):\n",
    "        newSet = np.transpose(np.random.permutation(np.transpose(Set))).tolist()\n",
    "        Sets.append(newSet)\n",
    "\n",
    "    #Returns concatenated list of sest:\n",
    "    H = np.concatenate(Sets)\n",
    "    return H\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "La matrice de parité régulière (12,3,4):\n",
      "\n",
      " [[1 1 1 1 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 1 1 1 1 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 1 1 1 1]\n",
      " [0 0 0 1 0 0 0 1 1 0 0 1]\n",
      " [0 1 0 0 0 1 1 0 0 1 0 0]\n",
      " [1 0 1 0 1 0 0 0 0 0 1 0]\n",
      " [0 1 0 0 1 0 0 0 1 0 1 0]\n",
      " [0 0 1 0 0 1 0 1 0 1 0 0]\n",
      " [1 0 0 1 0 0 1 0 0 0 0 1]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Test de la fonction \n",
    "RegularH_test = 1\n",
    "\n",
    "if RegularH_test:\n",
    "    Htest = RegularH(12,3,4)\n",
    "    print(\"La matrice de parité régulière (12,3,4):\\n\\n\",Htest )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### 1-2 Construction de la Matrice de codage G à partir de H\n",
    " \n",
    "\n",
    "### I - Le problème:\n",
    "Soit une matrice de parité H donnée. Trouver la matrice de codage correspondante c'est résoudre l'équation en G avec G' transposée de G:\n",
    "                                                \n",
    "                                                H.G' = 0 (1)\n",
    "                                                \n",
    "                          Sous contrainte rang(G) = n - rang(H) (1')\n",
    "                         \n",
    "La contrainte (1') vient de:\n",
    "\n",
    "<img src=\"Equations/Construction1.png\">\n",
    "\n",
    "> <font color=#A44057> rang(G) = dim(Ker(H)) </font> \n",
    "\n",
    "Par ailleurs, le théorème du rang appliqué à H donne: \n",
    "> <font color=#A44057> n = rang(H) + dim(ker(H)) </font> \n",
    " \n",
    " D'où (1')\n",
    " \n",
    " \n",
    " ### II -  Solutions \n",
    " \n",
    " \n",
    " \n",
    " On propose deux méthodes pour résoudre l'équation (1):\n",
    " \n",
    " - 1- En effectuant des opérations sur les lignes et les colonnes de H en utilisant les matrices équivalentes:\n",
    " on btient dans ce cas une matrice génératrice quelconque. \n",
    " \n",
    " - 2- En <b> changeant H </b> (en permutant ses colonnes) pour trouver une matrice de codage associée <b> systématique </b> , c'est à dire sous la forme: \n",
    " \n",
    "<img src=\"Equations/Construction2.png\">\n",
    " \n",
    " \n",
    " X représente la partie redondante de l'information. \n",
    "\n",
    "> <h3 align=\"center\" font=\"b\"> Solution 1: Matrice de Codage quelconque </h3>\n",
    "\n",
    "### La théorie \n",
    "\n",
    " En posant rang(G) = k, n-k est le rang H. Il existe P et Q inversibles de tailles respectives m (on rappelle que H a m lignes) et n telles que  P.H.Q = J où J est la matrice de taille m,n de rang n-k de la forme: \n",
    "                                        \n",
    " $$ \\begin{pmatrix} I_{n-k} & 0_{n-k,k} \\\\ 0_{m-(n-k),n-k} & 0_{m-(n-k),k} \\end{pmatrix} $$   \n",
    "                                       \n",
    "<img src=\"Equations/Construction3.png\">\n",
    "\n",
    " On résout l'équation (2): \n",
    " \n",
    "<img src=\"Equations/Construction4.png\">\n",
    "\n",
    " où Y1 est de taille (n-k,k) , et Y2 de taille (k,k). \n",
    "         \n",
    " \n",
    "Le produit bloc par bloc de (2) donne : Y1 = 0 , Y2 quelconque. Or (1') impose que G soit de rang k. comme Y1 = 0 et Y2 est de taille k, Y2 doit être inversible, il suffit de prendre:  Y2 = I_k.  \n",
    "\n",
    " On trouve G = (Q.Y)' Où  \n",
    " \n",
    " <img src=\"Equations/Construction5.png\">\n",
    "\n",
    "\n",
    "\n",
    "###  L'Algorithme: \n",
    "\n",
    "On a besoin de plusieurs fonctions: \n",
    "\n",
    "- BinaryProduct: Fonction qui effectue le produit de deux matrices dans M_{p,q}(F_2) .\n",
    "\n",
    "- GaussJordan: Fonction qui effectue les opérations sur les lignes d'une matrice A (forme échelonée réduite).\n",
    "   \n",
    "                                   |1 0 0 * * * * |\n",
    "                                   |0 1 0 * * * * |\n",
    "             Exemple:              |0 0 1 1 * * * |\n",
    "                                   |0 0 0 0 1 * * |\n",
    "                              \n",
    " Et renvoie aussi la matrice P telle que P.A = reduced(A) en effectuant les même opérations sur l'identité.\n",
    "                                   \n",
    " Cette fonction effectue un double pivot de Gauss: On parcourt les colonnes et dès qu'un pivot est trouvé, tous les 1 au-dessus et au-dessous sont annulés par addition de lignes. voir l'algorithme de Gauss-Jordan. https://fr.wikipedia.org/wiki/%C3%89limination_de_Gauss-Jordan. \n",
    " \n",
    "Cette fonction sera appliquée à H pour les lignes et à sa transposée pour les colonnes. \n",
    " \n",
    "- CodingMatrix: Fonction qui calcule la matrice Q en effectuant un Gauss-Jordan double puis renvoie G = (Q.Y)'.\n",
    "- BinaryRank: Fonction qui calcule le rang d'une matrice. (Basée sur Gauss-Jordan).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### La fonction BinaryProduct: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def BinaryProduct(X,Y):\n",
    "        \n",
    "    \"\"\" Binary Matrices or Matrix-vector product in Z/2Z. Works with scipy.sparse.csr_matrix matrices X,Y too.\"\"\"\n",
    " \n",
    "    A = X.dot(Y)\n",
    "    \n",
    "    if type(A)!=scipy.sparse.csr_matrix:\n",
    "        return A%2\n",
    "    \n",
    "    return A.toarray()%2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### La fonction GaussJordan:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def GaussJordan(MATRIX,change=0):\n",
    "\n",
    "    \"\"\" \n",
    "    Description:\n",
    "\n",
    "    Performs the row reduced echelon form of MATRIX and returns it.\n",
    "\n",
    "    If change = 1, all changes in the MATRIX's rows are applied to identity matrix P: \n",
    "\n",
    "    Let A be our parameter MATRIX. refA the reduced echelon form of A. P is the square invertible matrix:\n",
    "\n",
    "    P.A = Aref.\n",
    "\n",
    "    -------------------------------------------------\n",
    "    Parameters: \n",
    "\n",
    "    MATRIX: 2D-Array. \n",
    "    change : boolean (default = 0)\n",
    "\n",
    "    ------------------------------------------------\n",
    "\n",
    "    change = 0  (default)\n",
    "     >>> Returns 2D-Array Row Reduced Echelon form of Matrix\n",
    "\n",
    "    change = 1 \n",
    "    >>> Returns Tuple of 2D-arrays (refMATRIX, P) where P is described above.\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    A = np.copy(MATRIX)\n",
    "    m,n = A.shape\n",
    "    \n",
    "\n",
    "    if change:\n",
    "        P = np.identity(m).astype(int)\n",
    "\n",
    "    pivot_old = -1 \n",
    "    for j in range(n):            \n",
    "\n",
    "        filtre_down = A[pivot_old+1:m,j]\n",
    "        pivot = np.argmax(filtre_down)+pivot_old+1\n",
    "        \n",
    "\n",
    "        if A[pivot,j]:\n",
    "            pivot_old+=1 \n",
    "            if pivot_old != pivot:\n",
    "                aux = np.copy(A[pivot,:])\n",
    "                A[pivot,:] = A[pivot_old,:]\n",
    "                A[pivot_old,:] = aux\n",
    "                if change:\n",
    "                    aux = np.copy(P[pivot,:])\n",
    "                    P[pivot,:] = P[pivot_old,:]\n",
    "                    P[pivot_old,:] = aux\n",
    "\n",
    "            \n",
    "\n",
    "            for i in range(m):\n",
    "                if i!=pivot_old and A[i,j]:\n",
    "                    if change:\n",
    "                        P[i,:] = abs(P[i,:]-P[pivot_old,:])\n",
    "                    A[i,:] = abs(A[i,:]-A[pivot_old,:])\n",
    "\n",
    "\n",
    "        if pivot_old == m-1:\n",
    "            break\n",
    "\n",
    " \n",
    "    if change:    \n",
    "        return A,P \n",
    "    return A \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "La matrice Hrowreduced P.H: \n",
      "\n",
      "[[1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0]\n",
      " [0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1]\n",
      " [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]\n",
      " [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1]\n",
      " [0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 1 0]\n",
      " [0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0]\n",
      " [0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1]\n",
      " [0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 1 1]\n",
      " [0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1 0 1 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 1]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 1]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]] \n",
      " \n",
      "P.H = Hrowreduced ? True\n"
     ]
    }
   ],
   "source": [
    "####################### Test de la fonction:\n",
    "\n",
    "GaussJordan_test = 1\n",
    "\n",
    "if GaussJordan_test:\n",
    "    \n",
    "    Htest=RegularH(20,4,5)\n",
    "    \n",
    "    #### On test GaussJordan sur Htest avec l'argument 1 pour obtenir la matrice de passage P:\n",
    "    Hrowreduced, P = GaussJordan(Htest,1)\n",
    "    print(\"\\nLa matrice Hrowreduced P.H: \\n\\n{} \\n \".format(Hrowreduced))\n",
    "\n",
    "    #On vérife Hreducedtest = P.H \n",
    "    print(\"P.H = Hrowreduced ?\",(BinaryProduct(P,Htest)==Hrowreduced).all())\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Le double Gauss-Jordan\n",
    "Prenons une matrice de taille m,n avec n < m. Effectuer des opérations sur les lignes de H revient à multiplier à gauche H par une matrice P inversible (bien déterminée). Effectuer des opérations sur les colonnes de H c'est effectuer des opérations sur les lignes de sa transposée, et donc multiplier à gauche tH par une matrice tQ inversible. \n",
    "\n",
    "On commence par effectuer les opérations sur les colonnes, c'est-à-dire, on \"GaussJordan\" la transposée de H. Puis on GausseJordan la transposée de la matrice obtenue. Exemple:\n",
    "                     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "La matrice H:\n",
      " [[1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]\n",
      " [1 1 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0]\n",
      " [0 0 0 1 1 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0]\n",
      " [0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 1 1]\n",
      " [0 0 1 0 0 0 0 0 1 1 0 0 1 0 0 0 1 0 0 0]\n",
      " [0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 0 1 1 0]\n",
      " [1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 1 0 0 0]\n",
      " [0 0 1 1 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0]\n",
      " [0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 1]]\n",
      "\n",
      "La matrice transposée de H:\n",
      " [[1 0 0 0 1 0 0 0 0 1 0 0]\n",
      " [1 0 0 0 1 0 0 0 0 0 0 1]\n",
      " [1 0 0 0 0 0 0 1 0 0 1 0]\n",
      " [1 0 0 0 0 1 0 0 0 0 1 0]\n",
      " [1 0 0 0 0 1 0 0 1 0 0 0]\n",
      " [0 1 0 0 0 1 0 0 0 1 0 0]\n",
      " [0 1 0 0 0 0 1 0 0 0 0 1]\n",
      " [0 1 0 0 0 0 1 0 1 0 0 0]\n",
      " [0 1 0 0 0 0 0 1 0 1 0 0]\n",
      " [0 1 0 0 0 0 0 1 1 0 0 0]\n",
      " [0 0 1 0 1 0 0 0 0 0 0 1]\n",
      " [0 0 1 0 0 1 0 0 0 0 1 0]\n",
      " [0 0 1 0 0 0 0 1 0 0 1 0]\n",
      " [0 0 1 0 1 0 0 0 0 0 0 1]\n",
      " [0 0 1 0 1 0 0 0 0 0 1 0]\n",
      " [0 0 0 1 0 1 0 0 0 1 0 0]\n",
      " [0 0 0 1 0 0 0 1 0 1 0 0]\n",
      " [0 0 0 1 0 0 1 0 1 0 0 0]\n",
      " [0 0 0 1 0 0 1 0 1 0 0 0]\n",
      " [0 0 0 1 0 0 1 0 0 0 0 1]]\n",
      "\n",
      "RowReduced appliqué à transposée de H (comme çi on manipulait les colonnes de H):\n",
      " [[1 0 0 0 0 0 0 1 0 0 0 1]\n",
      " [0 1 0 0 0 0 0 1 0 0 0 1]\n",
      " [0 0 1 0 0 0 0 1 0 0 0 1]\n",
      " [0 0 0 1 0 0 0 1 0 0 0 1]\n",
      " [0 0 0 0 1 0 0 1 0 0 0 0]\n",
      " [0 0 0 0 0 1 0 1 0 0 0 0]\n",
      " [0 0 0 0 0 0 1 1 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 1 0 0 1]\n",
      " [0 0 0 0 0 0 0 0 0 1 0 1]\n",
      " [0 0 0 0 0 0 0 0 0 0 1 1]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0]]\n",
      "\n",
      "Puis on transpose cette matrice:\n",
      " [[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]\n",
      " [1 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0]]\n",
      "\n",
      " et on lui applique GaussJordan :\n",
      " [[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]]\n",
      "(12, 20)\n"
     ]
    }
   ],
   "source": [
    "Exemple = 1\n",
    "if Exemple:    \n",
    "    Htest = RegularH(20,3,5)\n",
    "    tHtest = np.transpose(Htest)\n",
    "    print(\"\\nLa matrice H:\\n\", Htest)\n",
    "    print(\"\\nLa matrice transposée de H:\\n\", tHtest)\n",
    "\n",
    "    tHrowreduced = GaussJordan(tHtest)\n",
    "\n",
    "    print(\"\\nRowReduced appliqué à transposée de H (comme çi on manipulait les colonnes de H):\\n\",tHrowreduced)\n",
    "    tHrowreducedT = np.transpose(tHrowreduced)\n",
    "    print(\"\\nPuis on transpose cette matrice:\\n\",tHrowreducedT)\n",
    "\n",
    "    rowreducedtHrowreducedT = GaussJordan(tHrowreducedT)\n",
    "\n",
    "    print(\"\\n et on lui applique GaussJordan :\\n\",rowreducedtHrowreducedT)\n",
    "\n",
    "    print(Htest.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fonction CodingMatrix "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CodingMatrix(MATRIX,use_sparse=1):\n",
    "\n",
    "    \"\"\" \n",
    "    CAUTION: RETURNS tG TRANSPOSED CODING MATRIX. \n",
    "    \n",
    "    Function Applies GaussJordan Algorithm on Columns and rows of MATRIX in order\n",
    "    to permute Basis Change matrix using Matrix Equivalence.\n",
    "\n",
    "    Let A be the treated Matrix. refAref the double row reduced echelon Matrix.\n",
    "\n",
    "    refAref has the form:\n",
    "\n",
    "    (e.g) : |1 0 0 0 0 0 ... 0 0 0 0|  \n",
    "            |0 1 0 0 0 0 ... 0 0 0 0| \n",
    "            |0 0 0 0 0 0 ... 0 0 0 0| \n",
    "            |0 0 0 1 0 0 ... 0 0 0 0| \n",
    "            |0 0 0 0 0 0 ... 0 0 0 0| \n",
    "            |0 0 0 0 0 0 ... 0 0 0 0| \n",
    "\n",
    "\n",
    "    First, let P1 Q1 invertible matrices: P1.A.Q1 = refAref\n",
    "\n",
    "    We would like to calculate:\n",
    "    P,Q are the square invertible matrices of the appropriate size so that:\n",
    "\n",
    "    P.A.Q = J.  Where J is the matrix of the form (having MATRIX's shape):\n",
    "\n",
    "    | I_p O | where p is MATRIX's rank and I_p Identity matrix of size p.\n",
    "    | 0   0 |\n",
    "\n",
    "    Therfore, we perform permuations of rows and columns in refAref (same changes\n",
    "    are applied to Q1 in order to get final Q matrix)\n",
    "\n",
    "\n",
    "    NOTE: P IS NOT RETURNED BECAUSE WE DO NOT NEED IT TO SOLVE H.G' = 0 \n",
    "    P IS INVERTIBLE, WE GET SIMPLY RID OF IT. \n",
    "    \n",
    "    Then\n",
    "    \n",
    "    solves: inv(P).J.inv(Q).G' = 0 (1) where inv(P) = P^(-1) and \n",
    "    P.H.Q = J. Help(PJQ) for more info.\n",
    "    \n",
    "    Let Y = inv(Q).G', equation becomes J.Y = 0 (2) whilst:\n",
    "    \n",
    "    J = | I_p O | where p is H's rank and I_p Identity matrix of size p.\n",
    "        | 0   0 |\n",
    "    \n",
    "    Knowing that G must have full rank, a solution of (2) is Y = |  0  | Where k = n-p. \n",
    "                                                                 | I-k |\n",
    "    \n",
    "    Because of rank-nullity theorem. \n",
    "    \n",
    "    -----------------\n",
    "    parameters:\n",
    "    \n",
    "    H: Parity check matrix. \n",
    "    use_sparse: (optional, default True): use scipy.sparse format to speed up calculations\n",
    "    ---------------\n",
    "    returns:\n",
    "    \n",
    "    tG: Transposed Coding Matrix. \n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    H = np.copy(MATRIX)\n",
    "    m,n = H.shape\n",
    "\n",
    "    if m > n: \n",
    "        raise ValueError('MATRIX must have more rows than columns (a parity check matrix)')\n",
    "    \n",
    "    if n > 500 and use_sparse:\n",
    "        sparse = 1\n",
    "    \n",
    "    else:\n",
    "        sparse = 0\n",
    "    ##### DOUBLE GAUSS-JORDAN:\n",
    "\n",
    "    Href_colonnes,tQ = GaussJordan(np.transpose(H),1)\n",
    "\n",
    "    Href_diag = GaussJordan(np.transpose(Href_colonnes))    \n",
    "\n",
    "    Q=np.transpose(tQ)\n",
    "    \n",
    "    k = n - sum(Href_diag.reshape(m*n))\n",
    "\n",
    "    \n",
    "    Y = np.zeros(shape=(n,k)).astype(int)\n",
    "    Y[n-k:,:] = np.identity(k)\n",
    "    \n",
    "    if sparse:\n",
    "        Q = csr_matrix(Q)\n",
    "        Y = csr_matrix(Y)\n",
    "\n",
    "    tG = BinaryProduct(Q,Y)\n",
    "\n",
    "    return tG\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def BinaryRank(MATRIX):\n",
    "    \"\"\" Computes rank of a binary Matrix using Gauss-Jordan algorithm\"\"\"\n",
    "    A = np.copy(MATRIX)\n",
    "    m,n = A.shape\n",
    "    \n",
    "    A = GaussJordan(A)\n",
    "    \n",
    "    return sum([a.any() for a in A])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Condition H.G' = 0 (1) Vérifiée ? True\n",
      "\n",
      "Contrainte de rang (1') vérifiée ? True\n"
     ]
    }
   ],
   "source": [
    "##### TEST DE LA FONCTION CODINGMATRIX:\n",
    "\n",
    "CodingMatrix_test=1\n",
    "\n",
    "if CodingMatrix_test:\n",
    "    n = 100\n",
    "    d_c = 4\n",
    "    d_v = 3\n",
    "    Htest = RegularH(n,d_v,d_c)\n",
    "    tGtest = CodingMatrix(Htest)\n",
    "    \n",
    "    print(\"Condition H.G' = 0 (1) Vérifiée ?\",(BinaryProduct(Htest,tGtest)==0).all())\n",
    "\n",
    "    print(\"\\nContrainte de rang (1') vérifiée ?\", BinaryRank(tGtest)==n -BinaryRank(Htest) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> <h3 align=\"center\" font=\"b\"> Solution 2: Matrice de Codage systématique </h3>\n",
    "\n",
    "### La théorie:\n",
    "\n",
    " <u><b> Lemme: H sous la forme systémtique et de rang maximal: </b></u>\n",
    " \n",
    " H s'écrit :\n",
    " \n",
    "<img src=\"Equations/Construction6.png\">\n",
    "\n",
    " Alors la matrice systématique :\n",
    " \n",
    "<img src=\"Equations/Construction7.png\">\n",
    " \n",
    " convient.\n",
    "\n",
    " <u><b> Cas particulier: Matrice de rang maximal K:</b></u> \n",
    "\n",
    " Si H est une matrice de parité de rang maximal m = n-k, alors, un simple Gauss-Jordan permet de l'écrire sous la forme: \n",
    "\n",
    "<img src=\"Equations/Construction8.png\">\n",
    "\n",
    "on note T la matrice inversible telle que:  refreduced(K) = TH  \n",
    "\n",
    "Or pour obtenir G sous la forme systématique, il faut que le bloc identité soit à droite, il faut donc appliquer une permutation aux colonnes de TH pour obtenir: \n",
    "\n",
    "<img src=\"Equations/Construction9.png\">\n",
    " \n",
    "Si l'on note par P la matrice de cette permutation:\n",
    "\n",
    "<img src=\"Equations/Construction10.png\">\n",
    "\n",
    "On rappelle que l'inverse de P est la transposée de P. \n",
    "\n",
    "Et donc\n",
    "\n",
    "<img src=\"Equations/Construction11.png\">\n",
    "\n",
    "Et comme H_s est systématique, d'après le lemme une solution GP' l'est aussi. \n",
    "Posons G_s = GP'. \n",
    "\n",
    "Il vient: \n",
    "\n",
    "<img src=\"Equations/Construction12.png\">\n",
    "\n",
    "Et donc on obtient finalement: \n",
    "\n",
    " -  Une matrice de parité H modifiée par l'inverse de la permutation P.\n",
    " -  Une matrice de codage G_s systématique. \n",
    "\n",
    " <u><b> Cas général: Matrice de parité de Gallager:</b></u> \n",
    "\n",
    "\n",
    "Comme H n'est pas forcément de rang maximal m, appliquer Gauss-Jordan à H donne une matrice de la forme:\n",
    "\n",
    "<img src=\"Equations/Construction13.png\">\n",
    "\n",
    "Il faut effectuer des permutations sur les colonnes de refreduced(H) = T.H pour l'écrire sous la forme: \n",
    "\n",
    "<img src=\"Equations/Construction14.png\">\n",
    "\n",
    "On note P_1 la première matrice relative à cette permutation donc\n",
    "\n",
    "<img src=\"Equations/Construction15.png\">\n",
    "\n",
    "Notez bien que H_{ss} s'écrit sous la forme du <b> cas particulier </b> avec des lignes nulles en plus.  \n",
    "\n",
    "Ces lignes ne changent en rien l'équation (1), et lui appliquant la permutation \n",
    "\n",
    "<img src=\"Equations/Construction10.png\">\n",
    "\n",
    "les permutations se combinent en une seule matrice:\n",
    "\n",
    "<img src=\"Equations/Construction16.png\">\n",
    "\n",
    "et comme H_s est systématique, d'après le lemme une solution GP' l'est aussi se construisant facilement à partir de H_s.\n",
    "\n",
    "Posons G_s = GP'. \n",
    "\n",
    "Il vient:\n",
    "\n",
    "<img src=\"Equations/Construction17.png\">\n",
    "\n",
    "Et donc on obtient finalement: \n",
    "\n",
    " -  Une matrice de parité H modifiée par l'inverse de la permutation P = P_2.P_1.\n",
    " -  Une matrice de codage $G_s$ systématique que l'on construit avec H_s.\n",
    " \n",
    " ### L'algorithme:\n",
    " \n",
    "La fonction CodingMatrix_systematic effectue toutes ces opérations: \n",
    "\n",
    "- Elle applique Gauss-Jordan à H. \n",
    "- Permute les dernières colonnes de la forme échelonée réduite de H pour obtenir le bloc identité d'abord à gauche.\n",
    "- Applique la permutation sigma pour obtenir le bloc identité à droite. \n",
    "- Récupère la matrice de permutation globale P et applique son inverse (transposée) à H pour renvoyer Hp.\n",
    "- Récupère le bloc A dans la forme systématique de H pour construire Gs. \n",
    "- Renvoie le couple Hp,Gs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CodingMatrix_systematic(MATRIX,use_sparse = 1):\n",
    "\n",
    "    \"\"\" \n",
    "    Description:\n",
    "\n",
    "    Solves H.G' = 0 and finds the coding matrix G in the systematic form : [I_k  A] by applying permutations on MATRIX.\n",
    "    \n",
    "    CAUTION: RETURNS TUPLE (Hp,tGS) WHERE Hp IS A MODIFIED VERSION OF THE GIVEN PARITY CHECK MATRIX, tGS THE TRANSPOSED \n",
    "    SYSTEMATIC CODING MATRIX ASSOCIATED TO Hp. YOU MUST USE THE RETURNED TUPLE IN CODING AND DECODING, RATHER THAN THE UNCHANGED \n",
    "    PARITY-CHECK MATRIX H. \n",
    "\n",
    "    -------------------------------------------------\n",
    "    Parameters: \n",
    "\n",
    "    MATRIX: 2D-Array. Parity-check matrix.\n",
    "    use_sparse: (optional, default True): use scipy.sparse matrices to speed up calculations if n>100.\n",
    "\n",
    "    ------------------------------------------------\n",
    "\n",
    "    >>> Returns Tuple of 2D-arrays (Hp,GS):\n",
    "        Hp: Modified H: permutation of columns (The code doesn't change)\n",
    "        tGS: Transposed Systematic Coding matrix associated to Hp.\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    H = np.copy(MATRIX)\n",
    "    m,n = H.shape\n",
    "    \n",
    "    if n>100 and use_sparse:\n",
    "        sparse = 1\n",
    "    else:\n",
    "        sparse = 0 \n",
    "        \n",
    "    P1 = np.identity(n,dtype=int)\n",
    "    \n",
    "    Hrowreduced = GaussJordan(H)\n",
    "    \n",
    "    k = n - sum([a.any() for a in Hrowreduced ])\n",
    "\n",
    "    ## After this loop, Hrowreduced will have the form H_ss : | I_(n-k)  A |\n",
    "    permut = np.array(list(range(n)))\n",
    "\n",
    "    while(True):\n",
    "        zeros = [i for i in range(min(m,n)) if not Hrowreduced[i,i]]\n",
    "        if len(zeros)==0:\n",
    "            break\n",
    "        indice_colonne_a = min(zeros)\n",
    "        list_ones = [j for j in range(indice_colonne_a+1,n) if Hrowreduced[indice_colonne_a,j] ]\n",
    "        if not len(list_ones):\n",
    "            break\n",
    "\n",
    "        indice_colonne_b = min(list_ones)\n",
    "        \n",
    "        aux = np.copy(Hrowreduced[:,indice_colonne_a])\n",
    "        Hrowreduced[:,indice_colonne_a] = Hrowreduced[:,indice_colonne_b]\n",
    "        Hrowreduced[:,indice_colonne_b] = aux \n",
    "        \n",
    "        aux = np.copy(P1[:,indice_colonne_a])\n",
    "        P1[:,indice_colonne_a] = P1[:,indice_colonne_b]\n",
    "        P1[:,indice_colonne_b] = aux\n",
    "        \n",
    "    ############ NOW, Hrowreduced has the form: | I_(n-k)  A | , the permutation above makes it look like : \n",
    "    ########### |A  I_(n-k)|\n",
    "    \n",
    "    P1 = P1.T\n",
    "    identity = list(range(n))\n",
    "    sigma = identity[n-k:]+identity[:n-k]\n",
    "    \n",
    "    P2 = np.zeros(shape=(n,n),dtype=int)\n",
    "    P2[identity,sigma] = np.ones(n)\n",
    "    \n",
    "    if sparse:\n",
    "        P1 = csr_matrix(P1)\n",
    "        P2 = csr_matrix(P2)\n",
    "        H = csr_matrix(H)\n",
    "\n",
    "    P = BinaryProduct(P2,P1)\n",
    "    \n",
    "    if sparse:\n",
    "        P = csr_matrix(P)\n",
    "        \n",
    "    Hp = BinaryProduct(H,np.transpose(P))\n",
    "\n",
    "    GS = np.zeros((k,n),dtype=int)\n",
    "    GS[:,:k] = np.identity(k)\n",
    "    GS[:,k:] = np.transpose(Hrowreduced[:n-k,n-k:])\n",
    "    \n",
    "    \n",
    "    return Hp,np.transpose(GS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Condition Hp.Gs' = 0 (1) Vérifiée ? True\n",
      "\n",
      "Contrainte de rang (1') vérifiée ? True\n",
      "\n",
      "Gs est-elle systématique ? True\n"
     ]
    }
   ],
   "source": [
    "Test_systematic = 1\n",
    "\n",
    "if Test_systematic:\n",
    "    n =100\n",
    "    d_v = 3\n",
    "    d_c = 4\n",
    "    H = RegularH(n,d_v,d_c)\n",
    "    t = time()\n",
    "    Hp,tGS = CodingMatrix_systematic(H)\n",
    "    k = BinaryRank(tGS)\n",
    "    \n",
    "    print(\"Condition Hp.Gs' = 0 (1) Vérifiée ?\",(BinaryProduct(Hp,tGS)==0).all())\n",
    "\n",
    "    print(\"\\nContrainte de rang (1') vérifiée ?\", k == n - BinaryRank(Hp) )\n",
    "    \n",
    "    print(\"\\nGs est-elle systématique ?\", (tGS[:k,:]==np.identity(k)).all() )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2- Codage et Décodage\n",
    "\n",
    "Fixons le triplet (n,d_v,d_c). \n",
    "On commence par générer une matrice de parité H régulière (n,d_v,d_c) de taille(m,n)  avec la méthode de Gallager. On obtient la matrice de codage G avec la fonction CodingMatrixG ci-dessus. On récupère k en accédant à la taille de G. \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1 CODAGE \n",
    "\n",
    "Prenons un message (suite de bits) aléatoire v de taille k. Le code transmis dans le canal est donc d=v.G . En appliquant la modulation BPSK, on travaille avec le vecteur émis de k-ième coordonnée x_k=(-1)^{d_k}. Une fois transmis dans le canal, on reçoit à l'autre bout un vecteur x + e où est un bruit blanc gaussien (si l'on suppose que le bruit du canal est additif). \n",
    "\n",
    "On définit le <i> Signal-Noise-Ratio (SNR)</i> par: \n",
    "Pour un bruit gaussien (toujours centré en 0) et de variance <img src=\"Equations/Construction18.png\">\n",
    "\n",
    "Le SNR est défini ci-dessus en décibels.\n",
    "\n",
    "On présente deux fonctions de codage: \n",
    "\n",
    "- La fonction Coding_random génère aléatoirement un message de taille k et le code. \n",
    "\n",
    "- La fonction Coding code un message v de taille k passé en argument. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Coding_random(tG,SNR):\n",
    "    \"\"\"\n",
    "    IMPORTANT: tG can be transposed coding matrix scipy.sparse.csr_matrix object to speed up calculations.\n",
    "\n",
    "    Randomly computes a k-bits message v, where G's shape is (k,n). And sends it\n",
    "    through the canal.\n",
    "    \n",
    "    Message v is passed to G: d = v,G. d is a n-vector turned into a BPSK modulated\n",
    "    vector x. Then Additive White Gaussian Noise (AWGN) is added. \n",
    "    \n",
    "    SNR is the Signal-Noise Ratio: SNR = 10log(1/variance) in decibels, where variance is the variance of the AWGN.\n",
    "    Remember: \n",
    "        1. d = v.G => tG.tv = td\n",
    "        2. x = BPSK(d) (or if you prefer the math: x = pow(-1,d) )\n",
    "        3. y = x + AWGN(0,snr) \n",
    "    \n",
    "    ----------------------------\n",
    "\n",
    "    Parameters:\n",
    "\n",
    "    tG: 2D-Array (OR scipy.sparse.csr_matrix) Transpose of Coding Matrix obtained from CodingMatrix functions.\n",
    "    \n",
    "    SNR: the Signal-Noise Ratio: SNR = 10log(1/variance) in decibels. \n",
    "        >> SNR = 10log(1/variance) \n",
    "    \n",
    "    -------------------------------\n",
    "\n",
    "    Returns\n",
    "\n",
    "    Tuple(v,y) (Returns v to keep track of the random message v)\n",
    "    \n",
    "    \"\"\"\n",
    "    n,k = tG.shape\n",
    "    \n",
    "    v=np.random.randint(2,size=k)\n",
    "\n",
    "    d = BinaryProduct(tG,v)\n",
    "    x=pow(-1,d)\n",
    "\n",
    "    sigma = 10**(-SNR/20)\n",
    "    \n",
    "    e = np.random.normal(0,sigma,size=n)\n",
    " \n",
    "    y=x+e\n",
    "\n",
    "    return(v,y)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "def Coding(tG,v,SNR):\n",
    "    \"\"\"\n",
    "    \n",
    "    IMPORTANT: if H is large, tG (transposed coding matrix) can be  scipy.sparse.csr_matrix object to speed up calculations.\n",
    "\n",
    "    Codes a message v with Coding Matrix G, and sends it through a noisy (default)\n",
    "    channel. \n",
    "    \n",
    "    G's shape is (k,n). \n",
    "\n",
    "    Message v is passed to tG: d = tG.tv d is a n-vector turned into a BPSK modulated\n",
    "    vector x. Then Additive White Gaussian Noise (AWGN) is added. \n",
    "    \n",
    "    SNR is the Signal-Noise Ratio: SNR = 10log(1/variance) in decibels, where variance is the variance of the AWGN.\n",
    "    Remember: \n",
    "    \n",
    "        1. d = v.G (or (td = tG.tv))\n",
    "        2. x = BPSK(d) (or if you prefer the math: x = pow(-1,d) )\n",
    "        3. y = x + AWGN(0,snr) \n",
    "\n",
    "    \n",
    "    Parameters:\n",
    "\n",
    "    tG: 2D-Array (OR scipy.sparse.csr_matrix)Transposed Coding Matrix obtained from CodingMatrix functions.\n",
    "    v: 1D-Array, k-vector (binary of course ..) \n",
    "    SNR: Signal-Noise-Ratio: SNR = 10log(1/variance) in decibels. \n",
    "\n",
    "    -------------------------------\n",
    "\n",
    "    Returns y\n",
    "\n",
    "    \"\"\"\n",
    "    n,k = tG.shape\n",
    "\n",
    "    if len(v)!= k:\n",
    "        raise ValueError(\" Size of message v must be equal to number of Coding Matrix G's rows \" )  \n",
    "    \n",
    "    d = BinaryProduct(tG,v)\n",
    "    x=pow(-1,d)\n",
    "\n",
    "    sigma = 10**(-SNR/20)\n",
    "    e = np.random.normal(0,sigma,size=n)\n",
    "\n",
    "\n",
    "    y=x+e\n",
    "\n",
    "    return y\n",
    "   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TEST DE CODAGE, SNR = 0 (Variance =1)\n",
      "\n",
      "CODAGE ALEATOIRE \n",
      "\n",
      "Le message aléatoire envoyé est:\n",
      " [0 0 0 0 1 0 1 1 1 0 1 1 1 1 0 1 0 1 1 0 0 0 1 0 1 0 0] \n",
      " \n",
      "     Le vecteur reçu y après transmission:\n",
      " [ 0.84356401  0.8485863  -0.67503432  1.13481934 -2.63565612 -1.26209385\n",
      "  0.95918682 -1.28426828  2.07781204  1.15412043  0.09061849 -0.08333776\n",
      " -1.07985257 -0.87670084  1.42841226 -2.85720088 -1.03686701 -0.24959313\n",
      " -0.96140121  0.63464732  3.06917118 -0.28327719  1.25665331 -2.14365771\n",
      " -1.98529827 -2.01421305  1.05257607 -0.83013099  1.26012478 -2.54079882\n",
      "  0.0842997   1.31318741 -0.0257579   0.73875371  0.90034267  1.63582258\n",
      "  0.05511692 -0.17850981  0.51879047  0.26590847  2.14331619  0.38223382\n",
      "  0.89601517  0.89529574 -0.937976    2.76909374 -2.01156082  1.34676893\n",
      "  0.82940958  0.85393581 -2.76727851 -0.51531004 -0.73242723  0.13610703\n",
      "  1.2730714   0.74358222  0.42102699 -0.98415196 -3.21547971 -0.66723172\n",
      "  0.8687968  -0.49842351 -0.26925289  2.17250146 -1.3305789   1.76958569\n",
      " -0.40448254  0.7521897  -0.15195302  0.92620834 -0.81147931 -1.37944994\n",
      " -1.83640086  1.33058825 -1.61559854  1.0477162  -2.38320571  1.77758644\n",
      " -0.38288842  0.92968981  1.09720278 -1.0920031  -1.16432399 -0.08001647\n",
      "  1.87162502 -1.73099136 -0.67870679  0.50720283 -0.64186245  1.13375053\n",
      " -1.80190257  2.25810483 -1.49102334  1.89363321 -2.27709189  0.99734289\n",
      " -1.6869807   0.22222977 -1.26043442 -0.82173542]\n",
      "\n",
      " _____________________________________________________________ \n",
      " CODAGE SPECIFIQUE:\n",
      "\n",
      "Le message à k bits vtest donné est:\n",
      " [1 0 1 1 0 1 1 1 1 0 1 0 0 1 1 1 0 0 0 0 0 1 0 0 1 1 1] \n",
      " \n",
      "     Le vecteur reçu y après transmission:\n",
      " [-1.62797564  0.77663291 -0.969397   -2.25085881  2.30780875 -0.50444616\n",
      " -1.67064666 -1.05894375 -1.99850641  1.25198001  0.25492153  1.94803927\n",
      " -0.32055496  0.25344396 -1.30489198 -1.77103313 -0.24000353  0.98390333\n",
      "  3.01806686  0.76248406  0.64885812 -0.43660328 -0.96654935  0.99626078\n",
      " -1.77851452 -1.75095913  0.0704118  -1.81998099  1.0312298  -1.20620389\n",
      "  0.45806882 -0.92250278 -0.41039815 -0.33524871  1.44756536 -0.12791318\n",
      " -1.45820269 -0.97669547 -0.57924878  1.90548343  0.56116202  0.32849929\n",
      " -1.27094006  0.25776131  0.21058773 -0.06037066  2.09717684  2.56702199\n",
      " -1.59391336  2.0390097  -0.33024849 -2.47163622  0.62478296 -2.15995569\n",
      "  0.41879649 -1.88357696 -1.46760289  0.46608449  1.43595817  1.24667171\n",
      "  0.6766702  -3.10733472 -1.22786443  0.63451299  0.77350162  1.3191263\n",
      " -0.77548778 -2.59114998  1.08816647 -1.02754142  2.72699767  0.02947916\n",
      "  3.09089658 -1.40044569  2.13254434 -0.44087587  0.0340262  -1.50965168\n",
      " -0.65461274  2.44363654  0.01269201 -0.06240705  1.02874002  1.04376966\n",
      "  1.68455902  0.05989818  1.34960075 -0.69838656 -0.5760999   0.57403642\n",
      " -0.63236633  1.65249184  0.51063698 -2.32600998  1.00674098  1.2174459\n",
      "  1.10740209  0.76458201 -1.6002973   1.09124614]\n"
     ]
    }
   ],
   "source": [
    "########### TEST DES FONCTIONS DE CODAGE:\n",
    "CODING_test = 1\n",
    "\n",
    "if CODING_test:\n",
    "    ntest =100\n",
    "    d_vtest = 3\n",
    "    d_ctest = 4\n",
    "    Htest = RegularH(ntest,d_vtest,d_ctest)\n",
    "    Htest,tGtest = CodingMatrix_systematic(Htest)\n",
    "    tGsparse = csr_matrix(tGtest)\n",
    "    \n",
    "    print(\"TEST DE CODAGE, SNR = 0 (Variance =1)\\n\")\n",
    "    print(\"CODAGE ALEATOIRE \\n\")\n",
    "    print(\"Le message aléatoire envoyé est:\\n {} \\n \\n \\\n",
    "    Le vecteur reçu y après transmission:\\n {}\".format(Coding_random(tGtest,0)[0],Coding_random(tGsparse,0)[1]))\n",
    "    print(\"\\n _____________________________________________________________ \\n CODAGE SPECIFIQUE:\\n\")\n",
    "    ktest = tGtest.shape[1]\n",
    "    vtest = np.random.randint(2,size=ktest)\n",
    "\n",
    "    print(\"Le message à k bits vtest donné est:\\n {} \\n \\n \\\n",
    "    Le vecteur reçu y après transmission:\\n {}\".format(vtest,Coding(tGsparse,vtest,0)))\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2-2 Décodage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 - 2 -1 Algorithme BP\n",
    "\n",
    "On introduit d'abord une série de fonctions qui serviront dans les fonctions de décodage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pi = math.pi\n",
    "def f1(y,sigma):\n",
    "    \"\"\" Normal Density N(1,sigma) \"\"\" \n",
    "    return(np.exp(-.5*pow((y-1)/sigma,2))/(sigma*math.sqrt(2*pi)))\n",
    "\n",
    "def fM1(y,sigma):\n",
    "    \"\"\" Normal Density N(-1,sigma) \"\"\" \n",
    "\n",
    "    return(np.exp(-.5*pow((y+1)/sigma,2))/(sigma*math.sqrt(2*pi)))\n",
    "\n",
    "def Bits2i(H,i):\n",
    "    \"\"\"\n",
    "    Computes list of elements of N(i)-j:\n",
    "    List of variables (bits) connected to Parity node i.\n",
    "\n",
    "    \"\"\"\n",
    "    if type(H)!=scipy.sparse.csr_matrix:\n",
    "        m,n=H.shape\n",
    "        return ([a for a in range(n) if H[i,a] ])\n",
    "    \n",
    "    indj = H.indptr\n",
    "    indi = H.indices\n",
    "    \n",
    "    return [indi[a] for a in range(indj[i],indj[i+1])]\n",
    "\n",
    "\n",
    "def Nodes2j(tH,j):\n",
    "    \n",
    "    \"\"\"\n",
    "    Computes list of elements of M(j):\n",
    "    List of nodes (PC equations) connecting variable j.\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    return Bits2i(tH,j) \n",
    "\n",
    "def InCode(H,x):\n",
    "\n",
    "    \"\"\" Computes Binary Product of H and x. If product is null, x is in the code.\n",
    "\n",
    "        Returns appartenance boolean. \n",
    "    \"\"\"\n",
    "        \n",
    "    return  (BinaryProduct(H,x)==0).all()\n",
    "\n",
    "def BitsAndNodes(H):\n",
    "    \n",
    "    m,n = H.shape\n",
    "    if type(H)==scipy.sparse.csr_matrix:\n",
    "        tH = scipy.sparse.csr_matrix(np.transpose(H.toarray()))\n",
    "    else:\n",
    "        tH = np.transpose(H)\n",
    "        \n",
    "    Bits = [Bits2i(H,i) for i in range(m)]\n",
    "    Nodes = [Nodes2j(tH,j)for j in range(n)]\n",
    "    \n",
    "    return Bits,Nodes\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### FONCTION DECODAGE - ALGORITHME BP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Decoding_BP(H,y,SNR,max_iter=1):\n",
    "\n",
    "    \"\"\" Decoding function using Belief Propagation algorithm.\n",
    "        \n",
    "        \n",
    "        IMPORTANT: H can be scipy.sparse.csr_matrix object to speed up calculations if n > 1000 highly recommanded. \n",
    "    -----------------------------------\n",
    "    Parameters:\n",
    "    \n",
    "    H: 2D-array (OR scipy.sparse.csr_matrix object) Parity check matrix, shape = (m,n) \n",
    "\n",
    "    y: n-vector recieved after transmission in the channel. (In general, returned \n",
    "    by Coding Function)\n",
    "\n",
    "\n",
    "    Signal-Noise Ratio: SNR = 10log(1/variance) in decibels of the AWGN used in coding.\n",
    "    \n",
    "    max_iter: (default = 1) max iterations of the main loop. Increase if decoding is not error-free.\n",
    "\n",
    "     \"\"\"\n",
    "        \n",
    "    m,n=H.shape\n",
    "    if not len(y)==n:\n",
    "        raise ValueError('Size of y must be equal to number of parity matrix\\'s columns n')\n",
    "\n",
    "    if m>=n:\n",
    "        raise ValueError('H must be of shape (m,n) with m < n')\n",
    "\n",
    "    sigma = 10**(-SNR/20)\n",
    "    \n",
    "    p0 = np.zeros(shape=n)\n",
    "    p0 = f1(y,sigma)/(f1(y,sigma) + fM1(y,sigma))\n",
    "    p1 = np.zeros(shape=n)\n",
    "    p1 = fM1(y,sigma)/(f1(y,sigma) + fM1(y,sigma))\n",
    "\n",
    "\n",
    "    #### ETAPE 0 : Initialization \n",
    "    q0 = np.zeros(shape=(m,n))\n",
    "    q0[:] = p0\n",
    "\n",
    "    q1 = np.zeros(shape=(m,n))\n",
    "    q1[:] = p1\n",
    "\n",
    "    r0 = np.zeros(shape=(m,n))\n",
    "    r1 = np.zeros(shape=(m,n))\n",
    "\n",
    "    count=0\n",
    "    prod = np.prod\n",
    "    \n",
    "    Bits,Nodes = BitsAndNodes(H)\n",
    "\n",
    "    \n",
    "    while(True):\n",
    "\n",
    "        count+=1 #Compteur qui empêche la boucle d'être infinie .. \n",
    "\n",
    "        #### ETAPE 1 : Horizontale\n",
    "\n",
    "        deltaq = q0 - q1\n",
    "        deltar = r0 - r1 \n",
    "        \n",
    "\n",
    "        for i in range(m):\n",
    "            Ni=Bits[i]\n",
    "            for j in Ni:\n",
    "                Nij = Ni.copy()\n",
    "\n",
    "                if j in Nij: Nij.remove(j)\n",
    "                deltar[i,j] = prod(deltaq[i,Nij])\n",
    "                \n",
    "\n",
    "        r0 = 0.5*(1+deltar)\n",
    "        r1 = 0.5*(1-deltar)\n",
    "\n",
    "\n",
    "        #### ETAPE 2 : Verticale\n",
    "\n",
    "        for j in range(n):\n",
    "            Mj = Nodes[j]\n",
    "            for i in Mj:\n",
    "                Mji = Mj.copy()\n",
    "                if i in Mji: Mji.remove(i)\n",
    "                    \n",
    "                q0[i,j] = p0[j]*prod(r0[Mji,j])\n",
    "                q1[i,j] = p1[j]*prod(r1[Mji,j])\n",
    "                \n",
    "                if q0[i,j] + q1[i,j]==0:\n",
    "                    q0[i,j]=0.5\n",
    "                    q1[i,j]=0.5\n",
    "              \n",
    "                else:\n",
    "                    alpha=1/(q0[i,j]+q1[i,j]) #Constante de normalisation alpha[i,j] \n",
    "\n",
    "                    q0[i,j]*= alpha\n",
    "                    q1[i,j]*= alpha # Maintenant q0[i,j] + q1[i,j] = 1\n",
    "\n",
    "\n",
    "        #### Calcul des probabilites à posteriori:\n",
    "        q0_post = np.zeros(n)\n",
    "        q1_post = np.zeros(n)\n",
    "        \n",
    "        for j in range(n):\n",
    "            Mj=Nodes[j]\n",
    "            q0_post[j] = p0[j]*prod(r0[Mj,j])\n",
    "            q1_post[j] = p1[j]*prod(r1[Mj,j])\n",
    "            \n",
    "            if q0_post[j] + q1_post[j]==0:\n",
    "                q0_post[j]=0.5\n",
    "                q1_post[j]=0.5\n",
    "                \n",
    "            alpha = 1/(q0_post[j]+q1_post[j])\n",
    "            \n",
    "            q0_post[j]*= alpha\n",
    "            q1_post[j]*= alpha\n",
    "        \n",
    "\n",
    "        x = np.array(q1_post > q0_post).astype(int)\n",
    "        \n",
    "        if InCode(H,x) or count >= max_iter:  \n",
    "            break\n",
    "  \n",
    "    return x\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2-2 Algorithme full-log BP "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def Decoding_logBP(H,y,SNR,max_iter=1):\n",
    "    \n",
    "    \"\"\" Decoding function using Belief Propagation algorithm (logarithmic version)\n",
    "\n",
    "    IMPORTANT: if H is large (n>1000), H should be scipy.sparse.csr_matrix object to speed up calculations\n",
    "    (highly recommanded. )\n",
    "    -----------------------------------\n",
    "    Parameters:\n",
    "    \n",
    "    H: 2D-array (OR scipy.sparse.csr_matrix object) Parity check matrix, shape = (m,n) \n",
    "\n",
    "    y: n-vector recieved after transmission in the channel. (In general, returned \n",
    "    by Coding Function)\n",
    "\n",
    "\n",
    "    Signal-Noise Ratio: SNR = 10log(1/variance) in decibels of the AWGN used in coding.\n",
    "    \n",
    "    max_iter: (default = 1) max iterations of the main loop. Increase if decoding is not error-free.\n",
    "\n",
    "     \"\"\"\n",
    "        \n",
    "    m,n=H.shape\n",
    "\n",
    "    if not len(y)==n:\n",
    "        raise ValueError('La taille de y doit correspondre au nombre de colonnes de H')\n",
    "\n",
    "    if m>=n:\n",
    "        raise ValueError('H doit avoir plus de colonnes que de lignes')\n",
    "    \n",
    "    \n",
    "    var = 10**(-SNR/10)\n",
    "\n",
    "    ### ETAPE 0: initialisation \n",
    "    \n",
    "    Lc = 2*y/var\n",
    "\n",
    "    Lq=np.zeros(shape=(m,n))\n",
    "\n",
    "    Lr = np.zeros(shape=(m,n))\n",
    "    \n",
    "    count=0\n",
    "    \n",
    "    prod=np.prod\n",
    "    tanh = np.tanh\n",
    "    log = np.log\n",
    "    \n",
    "    \t\n",
    "    Bits,Nodes = BitsAndNodes(H)\n",
    "    while(True):\n",
    "\n",
    "        count+=1 #Compteur qui empêche la boucle d'être infinie .. \n",
    "\n",
    "        #### ETAPE 1 : Horizontale\n",
    "        for i in range(m):\n",
    "            Ni = Bits[i]\n",
    "            for j in Ni:\n",
    "                Nij = Ni.copy()\n",
    "\n",
    "                if j in Nij: Nij.remove(j)\n",
    "            \n",
    "                if count==1:\n",
    "                    X = prod(tanh(0.5*Lc[Nij]))\n",
    "                else:\n",
    "                    X = prod(tanh(0.5*Lq[i,Nij]))\n",
    "                num = 1 + X\n",
    "                denom = 1 - X\n",
    "                if num == 0: \n",
    "                    Lr[i,j] = -1\n",
    "                elif denom  == 0:\n",
    "                    Lr[i,j] =  1\n",
    "                else: \n",
    "                    Lr[i,j] = log(num/denom)\n",
    "              \n",
    "\n",
    "        #### ETAPE 2 : Verticale\n",
    "        for j in range(n):\n",
    "            Mj = Nodes[j]\n",
    "            \n",
    "            for i in Mj:\n",
    "                Mji = Mj.copy()\n",
    "                if i in Mji: Mji.remove(i)\n",
    "\n",
    "                Lq[i,j] = Lc[j]+sum(Lr[Mji,j])\n",
    "        \n",
    " \n",
    "        #### LLR a posteriori:\n",
    "        L_posteriori = np.zeros(n)\n",
    "        for j in range(n):\n",
    "            Mj = Nodes[j]\n",
    "\n",
    "            L_posteriori[j] = Lc[j] + sum(Lr[Mj,j])\n",
    "\n",
    "        x = np.array(L_posteriori <= 0).astype(int)\n",
    "            \n",
    "        product = InCode(H,x)\n",
    "\n",
    "        if product or count >= max_iter:  \n",
    "            break\n",
    "        \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2-3 Réception du message à k bits après décodage\n",
    "Après décodage et correction des erreurs du vecteur à n bits transmis, pour conclure quant à la réussite de la transmission, on doit pouvoir retrouver le message à k bits codé avec la matrice G à partir du vecteur décodé x. \n",
    "On sait que pour un vecteur donné x appartenant au code, il existe v vecteur à k bits tel que :\n",
    "    \n",
    "                        x = v . G  \n",
    "Donc: \n",
    "\n",
    "                        G'.v' = x'\n",
    "\n",
    "On propose de résoudre ce système en v en lui appliquant l'algorithme du pivot de Gauss codé ci-dessous dans la fonction GaussElimination. On sait que ce système a une solution car G est de rang k. Il s'agit d'un système de rang k à k inconnues. La solution existe et est unique.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def GaussElimination(MATRIX,B):\n",
    "\n",
    "    \"\"\" Applies Gauss Elimination Algorithm to MATRIX in order to solve a linear system MATRIX.X = B. \n",
    "    MATRIX is transformed to row echelon form: \n",
    "\n",
    "         |1 * * * * * |\n",
    "         |0 1 * * * * |\n",
    "         |0 0 1 * * * |\n",
    "         |0 0 0 1 * * | \n",
    "         |0 0 0 0 1 * |\n",
    "         |0 0 0 0 0 1 |\n",
    "         |0 0 0 0 0 0 |\n",
    "         |0 0 0 0 0 0 |\n",
    "         |0 0 0 0 0 0 |\n",
    "\n",
    "\n",
    "    Same row operations are applied on 1-D Array vector B. Both arguments are sent back.\n",
    "    \n",
    "    --------------------------------------\n",
    "    \n",
    "    Parameters:\n",
    "    \n",
    "    MATRIX: 2D-array. \n",
    "    B:      1D-array. Size must equal number of rows of MATRIX.\n",
    "            \n",
    "    -----------------------------------\n",
    "    Returns:\n",
    "    \n",
    "    Modified arguments MATRIX, B as described above.\n",
    "    \n",
    "         \"\"\"\n",
    "\n",
    "    A = np.copy(MATRIX)\n",
    "    b = np.copy(B)\n",
    "    n,k = A.shape\n",
    "\n",
    "\n",
    "    if b.size != n:\n",
    "        raise ValueError('Size of B must match number of rows of MATRIX to solve MATRIX.X = B')\n",
    "\n",
    "    for j in range(min(k,n)):\n",
    "        listeDePivots = [i for i in range(j,n) if A[i,j]]\n",
    "        if len(listeDePivots)>0:\n",
    "            pivot = np.min(listeDePivots)\n",
    "        else:\n",
    "            continue\n",
    "        if pivot!=j:\n",
    "            aux = np.copy(A[j,:])\n",
    "            A[j,:] = A[pivot,:]\n",
    "            A[pivot,:] = aux\n",
    "\n",
    "            aux = np.copy(b[j])\n",
    "            b[j] = b[pivot]\n",
    "            b[pivot] = aux\n",
    "\n",
    "        for i in range(j+1,n):\n",
    "            if A[i,j]:\n",
    "                A[i,:] = abs(A[i,:]-A[j,:])\n",
    "                b[i] = abs(b[i]-b[j])\n",
    "\n",
    "    return A,b\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction DecodedMessage applique l'algorithme du pivot de Gauss au système linéaire précédent et retrouve les k bits du message v en remontant. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def DecodedMessage(tG,x):\n",
    "\n",
    "    \"\"\"\n",
    "    Let G be a coding matrix. tG its transposed matrix. x a n-vector received after decoding.\n",
    "    DecodedMessage Solves the equation on k-bits message v:  x = v.G => G'v'= x' by applying GaussElimination on G'.\n",
    "    \n",
    "    -------------------------------------\n",
    "    \n",
    "    Parameters:\n",
    "    \n",
    "    tG: Transposed Coding Matrix. Must have more rows than columns to solve the linear system. Must be full rank.\n",
    "    x: n-array. Must be in the Code (in Ker(H)). \n",
    "\n",
    "    \"\"\"\n",
    "    n,k = tG.shape \n",
    "    \n",
    "    if n < k:\n",
    "        raise ValueError('Coding matrix G must have more columns than rows to solve the linear system on v\\': G\\'v\\' = x\\'')\n",
    "    \n",
    "                         \n",
    "    rtG, rx = GaussElimination(tG,x)\n",
    "    \n",
    "    rank = sum([a.any() for a in rtG])\n",
    "    if rank!= k:\n",
    "        raise ValueError('Coding matrix G must have full rank = k to solve G\\'v\\' = x\\'')\n",
    "            \n",
    "    message=np.zeros(k).astype(int)\n",
    "\n",
    "    message[k-1]=rx[k-1]\n",
    "    for i in reversed(range(k-1)):\n",
    "        message[i]=abs(rx[i]-BinaryProduct(rtG[i,list(range(i+1,k))],message[list(range(i+1,k))]))\n",
    "\n",
    "    return message\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##  3 -Satistiques:  Test du décodage \n",
    "On commence par générer une matrice de parité régulière et sa matrice de codage après avoir fixé les paramètres.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "TestDecodage = 1\n",
    "if TestDecodage:   \n",
    "    n= 300\n",
    "    d_v = 3\n",
    "    d_c  = 4\n",
    "    H = RegularH(n,d_v,d_c)\n",
    "    H,tG = CodingMatrix_systematic(H)\n",
    "    tGs = sparse.csr_matrix(tG)\n",
    "    Hs = sparse.csr_matrix(H)\n",
    "k = tG.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "On teste le décodage en transmettant 10 messages aléatoires et comparant le message aléatoire à k bits envoyé v_envoye au vecteur v_recu obtenu après décodage complet. On récupère dans la liste \"erreurs\" le nombre de bits erronés pour chaque message. La variable \"good\" filtre cette liste pour n'en garder que les messages décodés sans erreurs et enregistre sa taille, c'est-à-dire le nombre de décodages réussis.\n",
    "\n",
    "On utilise ici une matrice assez large pour montrer l'efficacité du format csr."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3 - 1 Test du décodage BP\n",
    "Pour un canal avec un bruit modélisé par un bruit gaussien avec un SNR = 8db \n",
    "\n",
    "Matrices classiques:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 100%\n",
      "\n",
      " Temps de décodage classique: 1.1923930644989014\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    snr = 8\n",
    "    max_iter = 1\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(10): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_BP(H,y,snr,max_iter)\n",
    "        v_recu = DecodedMessage(tG,x)\n",
    "        erreurs.append(sum(abs(v_envoye-v_recu)))\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0])*10 #Nombre de messages bien décodés\n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de décodage classique:\",t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 100%\n",
      "\n",
      " Temps de décodage sparse: 0.935791015625\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    max_iter = 1\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(10): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_BP(Hs,y,snr,max_iter)\n",
    "        v_recu = DecodedMessage(tG,x)\n",
    "        erreurs.append(sum(abs(v_envoye-v_recu)))\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0])*10 #Nombre de messages bien décodés\n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de décodage sparse:\",t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour un canal avec un bruit modélisé par un bruit gaussien avec un SNR = 0db (Variance =1), therefore we need more iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 50%\n",
      "\n",
      " Temps de décodage classique: 2.242175817489624\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    snr = 0\n",
    "    max_iter = 5\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(10): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_BP(H,y,snr,max_iter)\n",
    "        v_recu = DecodedMessage(tG,x)\n",
    "        erreurs.append(sum(abs(v_envoye-v_recu)))\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0])*10 #Nombre de messages bien décodés\n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de décodage classique:\",t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 60%\n",
      "\n",
      " Temps de décodage sparse: 2.653632879257202\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    max_iter=5\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(10): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_BP(Hs,y,snr,max_iter)\n",
    "        v_recu = DecodedMessage(tG,x)\n",
    "        erreurs.append(sum(abs(v_envoye-v_recu)))\n",
    "\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0])*10 #Nombre de messages bien décodés\n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de décodage sparse:\",t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3_2 Test du décodage full-log BP\n",
    "\n",
    "Pour un bruit gaussien avec SNR = 8 db."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 100%\n",
      "\n",
      " Temps de codage/décodage classique: 0.970768928527832\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    snr = 8\n",
    "    max_iter = 1\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(10): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_logBP(H,y,snr,max_iter)\n",
    "        v_recu = DecodedMessage(tG,x)\n",
    "        erreurs.append(sum(abs(v_recu-v_envoye)))\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0])*10 #Nombre de messages bien décodés\n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de codage/décodage classique:\",t)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse matrices: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 100%\n",
      "\n",
      " Temps de codage/décodage sparse: 0.7703440189361572\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    max_iter=1\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(10): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_logBP(Hs,y,snr,max_iter)\n",
    "        v_recu = DecodedMessage(tG,x)\n",
    "        erreurs.append(sum(abs(v_envoye-v_recu)))\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0])*10 #Nombre de messages bien décodés\n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de codage/décodage sparse:\",t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour un canal avec un bruit modélisé par un bruit gaussien avec snr = 0, on a donc besoin de plus d'une itération:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 70%\n",
      "\n",
      " Temps de codage/décodage classique: 2.1588029861450195\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    snr = 0\n",
    "    max_iter = 6\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(10): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_logBP(H,y,snr,max_iter)\n",
    "        v_recu = DecodedMessage(tG,x)\n",
    "        erreurs.append(sum(abs(v_recu-v_envoye)))\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0])*10 #Nombre de messages bien décodés\n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de codage/décodage classique:\",t)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 80%\n",
      "\n",
      " Temps de codage/décodage sparse: 2.1334900856018066\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    max_iter=6\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(10): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_logBP(Hs,y,snr,max_iter)\n",
    "        v_recu = DecodedMessage(tG,x)\n",
    "        erreurs.append(sum(abs(v_envoye-v_recu)))\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0])*10 #Nombre de messages bien décodés\n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de codage/décodage sparse:\",t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "BP et Log_BP sont-ils équivalents ? En utilisant H et Hsparse ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Le taux de décodage réussi est: 100%\n",
      "\n",
      " Temps de codage/décodage sparse: 35.69674587249756\n"
     ]
    }
   ],
   "source": [
    "if 1:    \n",
    "    max_iter=6\n",
    "    erreurs=[]\n",
    "    t = time()\n",
    "    for i in range(100): \n",
    "        v_envoye,y = Coding_random(tG,snr) \n",
    "        x = Decoding_logBP(H,y,snr,max_iter)\n",
    "        xx = Decoding_BP(Hs,y,snr,max_iter)\n",
    "        erreurs.append(sum(abs(xx-x)))\n",
    "    t = time() - t\n",
    "    good = len([i for i in erreurs if i==0]) #Nombre de messages égaux \n",
    "    print(\"Le taux de décodage réussi est: {}%\".format(good))\n",
    "    print(\"\\n Temps de codage/décodage sparse:\",t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4 - Applications: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}