path: root/PVCM/cama/fr/ma1 np05 Notation Einstein.ipynb

{
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Présentation de la notation d'Einstein\n",
    "\n",
    "Cette notation permet d'exprimer de facon concise un nombre extraordinaire d'opérations. Cette section n'est pas\n",
    "nécessaire pour la suite mais c'est un exercice intellectuel intéressant.\n",
    "\n",
    "L'idée est base est de sommer les termes d'une équation lorsqu'on y trouve 2 fois un même indice qui n'est pas défini par ailleurs. Ainsi :\n",
    "\n",
    "$$ A_{i,i} \\quad \\textrm{veut dire} \\quad \\sum_{i=1}^N A_{i,i} \\; \\textrm{(la trace de la matrice)}$$\n",
    "\n",
    "Si on regarde la multiplication matricielle $A\\, B$ on a pour tout indice (i,j) du résultat :\n",
    "\n",
    "$$ C_{i,j} = A_{i,k} \\, B_{k,j} \\quad \\textrm{c.a.d} \\quad  C_{i,j} = \\sum_{k=1}^N A_{i,k} \\, B_{k,j} $$ \n",
    "\n",
    "Le nom complet de la notation d'Einstein étant la convention de sommation d'Einstein, la fonction de Numpy est [`einsum`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html). Voici ce que cela donne pour nos 2 premiers exemples :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trace 'ii' :  12 \n",
      "\n",
      "Multiplication matricielle A A 'ij,jk->ik' :\n",
      "[[ 15  18  21]\n",
      " [ 42  54  66]\n",
      " [ 69  90 111]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.arange(9).reshape(3,3)\n",
    "\n",
    "print(\"Trace 'ii' : \", np.einsum('ii', A), '\\n')   # 0 + 4 + 8 = 12\n",
    "\n",
    "print(\"Multiplication matricielle A A 'ij,jk->ik' :\")\n",
    "print(np.einsum('ij,jk->ik', A, A))  # notez que j'ai nommé différement les indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 15,  18,  21],\n",
       "       [ 42,  54,  66],\n",
       "       [ 69,  90, 111]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.dot(A)  # on vérifie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "On note que les arguments de `einsum` sont :\n",
    "\n",
    "* la règle de sommation dans une chaine de caractères avec une virgule pour séparer chaque composant\n",
    "* les composants sur lesquels s'applique la règle\n",
    "\n",
    "On peut aller un peu plus loin avec Numpy. Voici l'ensemble des règles qu'utilise `einsum` :\n",
    "\n",
    "#### Règles de base et complémentaires\n",
    "\n",
    "1. un indice répété implique une sommation sur cet indice sauf si cet indice est cité dans le résultat<br/>(cf exemple diagonale de A ci-dessous pour l'exception)\n",
    "2. les indices qui se répètent d'un composant à un autre impliquent que les éléments référencés seront multipliés entre eux<br/> (cf exemple du produit matriciel)\n",
    "3. un lettre omise dans le résultat (après `->`) implique que l'on sommme les éléments sur cet indice<br/> (cf exemple de la somme des éléments d'un vecteur ci-dessous)\n",
    "4. si on ne met pas la flèche, einsum la met pour vous avec à droite tous les indices qui ne sont pas doublés rangés par ordre alphabétique<br/> (cf exemple de la transposée ci-dessous)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Voici une liste d'opérations prises sur le [blog du Dr Goulu](https://www.drgoulu.com/2016/01/17/einsum) :\n",
    "\n",
    "<table align=\"center\">\n",
    "<tbody>\n",
    "<tr>\n",
    "<th>Signature de einsum</th>\n",
    "<th>équivalent NumPy</th>\n",
    "<th>Description</th>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('i-&gt;', v)</code></td>\n",
    "<td><code>sum(v)</code></td>\n",
    "<td>somme des valeurs du vecteur v</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('i,i-&gt;i', u, v)</code></td>\n",
    "<td><code>u \\* v</code></td>\n",
    "<td>multiplication des vecteurs u et v élément par élément</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('i,i', u, v)</code></td>\n",
    "<td><code>inner(u, v)</code></td>\n",
    "<td>produit scalaire de u et v</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('i,j', u, v)</code></td>\n",
    "<td><code>outer(u, v)</code></td>\n",
    "<td>produit dyadique de u et v</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij', A)</code></td>\n",
    "<td><code>A</code></td>\n",
    "<td>renvoie la matrice A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ji', A)</code></td>\n",
    "<td><code>A.T</code></td>\n",
    "<td>transposée de A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ii-&gt;i', A)</code></td>\n",
    "<td><code>diag(A)</code></td>\n",
    "<td>diagonale de A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ii', A)</code></td>\n",
    "<td><code>trace(A)</code></td>\n",
    "<td>somme de la diagonale de A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td class=\"tg-031e\"><code>('ij-&gt;', A)</code></td>\n",
    "<td class=\"tg-031e\"><code>sum(A)</code></td>\n",
    "<td class=\"tg-031e\">somme des valeurs de A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij-&gt;j', A)</code></td>\n",
    "<td><code>sum(A, axis=0)</code></td>\n",
    "<td>somme des colonnes de A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij-&gt;i', A)</code></td>\n",
    "<td><code>sum(A, axis=1)</code></td>\n",
    "<td>somme des lignes de A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,ij-&gt;ij', A, B)</code></td>\n",
    "<td><code>A \\* B</code></td>\n",
    "<td>multiplication des matrices A et B élément par élément</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,ji-&gt;ij', A, B)</code></td>\n",
    "<td><code>A \\* B.T</code></td>\n",
    "<td>multiplication des matrices A et B transposée élément par élément</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,jk', A, B)</code></td>\n",
    "<td><code>dot(A, B)</code></td>\n",
    "<td>produit scalaire de A et B</td> </tr>\n",
    "<tr>\n",
    "<td><code>('ij,jk-&gt;ij', A, B)</code></td>\n",
    "<td><code>inner(A, B)</code></td>\n",
    "<td>produit intérieur de A et B</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,jk-&gt;ijk', A, B)</code></td>\n",
    "<td><code>A[:, None] \\* B</code></td>\n",
    "<td>chaque ligne de A multipliée par B</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,kl-&gt;ijkl', A, B)</code></td>\n",
    "<td><code>A[:, :, None, None] \\* B</code></td>\n",
    "<td>chaque valeur de A multipliée par B</td>\n",
    "</tr>\n",
    "</tbody>\n",
    "</table>\n",
    "\n",
    "Le `None` dans les deux dernières lignes est une facon de redimensionner un tableau. Ainsi `np.arange(6)` est un tableau de dimension 1, `np.arange(6)[:]` est le même tableau de dimension 1, alors que `np.arange(6)[:,None]` est un tableau de dimension 2 à savoir `6 x 1` quand `np.arange(6)[None,:,None]` a 3 dimensions : `1 x 6 x 1`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Mise en pratique\n",
    "\n",
    "On va comparer les performances d'`einsum` et des fonctions Numpy correspondante. Pour cela calculez\n",
    "\n",
    "* le cube de chaque élement d'un vecteur\n",
    "* d'une matrice carrée, $A^3$,\n",
    "\n",
    "avec `einsum` et sans. Comparez la vitesse d'exécution dans tous les cas avec 10000 éléments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a vous de jouer\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.5 µs ± 112 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      "15.2 µs ± 131 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "u = np.random.random(10000)\n",
    "\n",
    "%timeit u*u*u\n",
    "%timeit np.einsum('i,i,i->i', u,u,u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "138 µs ± 9.48 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "134 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "A = u.reshape(100,100)\n",
    "\n",
    "%timeit A.dot(A).dot(A)\n",
    "%timeit np.einsum('ij,jk,kl->il', A, A, A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "On constate qu'`einsum` est plus lent (\\*). Mais en regardant sur le web on constate que cela n'a pas toujours été le cas et que c'est lié à la version de la bibliothèque Numpy. Conclusion si on veut de la performance, il faut tester son code avant pour choisir la méthode la plus rapide.\n",
    "\n",
    "(\\*) un peu plus lent pour le calcul vectoriel mais 1000 fois plus lent sur mon portable pour `A.dot(A)` (tous mes processeurs tournent à 100% alors que le calcul par `einsum` n'est fait que sur 1 seul processeur). La version `A.dot(A)` est nettement plus rapide grace à la bibliothèque MKL qu'utilise Numpy sur ma machine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}