Table des matières

L'objectif de ce TD est d'introduire la visualisation avec Matplotlib et le calcul scientifique avec Numpy. Nous implémenterons une fonction simple avec différentes méthodes de programmation, mesurerons les performances de chacune et visualiserons les résultats.

Méthodes de programmation classiques

Récursion

Développez une fonction fibonacci_rec(n) qui calcule la valeur du n-ième terme de la suite de Fibonacci. Cette fonction doit être récursive. Les résultats attendus pour les premières valeurs sont : {$ 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181... $}
Remarque pour les étudiants avancés : ne pas utiliser la mémoïsation.

---

ANSWER ELEMENTS

Recall the definition of recursive functions.

def fibonacci_rec(n):
    if n == 0:
        return 0
    if n == 1:
        return 1
    else:
        return fibonacci_rec(n-1)+fibonacci_rec(n-2)

---

Analyse du temps d'exécution

La fonction timeit, définie dans le package homonyme, est utilisée pour mesurer le temps d'exécution d'un code Python. Dans l'exemple suivant, timeit renvoie le temps total de 100 exécutions du code Python pour lisser les résultats.

Les temps d'exécution peuvent varier d'une machine à l'autre ; vous devez vous attendre à un temps d'exécution de quelques millièmes de seconde. Essayez-le.

from timeit import timeit
timeit(lambda : fibonacci_rec(10), number=100)

Nous avons défini une fonction measure() qui détermine le temps d'exécution d'une fonction pour une séquence de valeurs.
Lorsque nous mesurons les performances de fibonacci() avec différentes valeurs de {$n$}, nous devrions observer un temps d'exécution élevé pour des valeurs au-delà de {$n = 20$}.

def measure(func, values, number=100):
    m = []
    for i in values :
        m.append(timeit(lambda : func(i), number=number))
    return m

measure(fibonacci_rec, [1, 2, 3, 5, 10, 20, 25])

Itération

Développez une fonction fibonacci_iter(n) qui calcule la valeur du n-ième terme de la suite de Fibonacci. Cette fonction doit être itérative.
Mesurez le temps d'exécution de fibonacci_iter(n). Il devrait être significativement plus rapide que la version récursive.

---

ANSWER ELEMENTS

def fibonacci_iter(n):
    if n == 0:
        return 0
    else:
        f = [0, 1]
        for i in range(2, n + 1):
            f.append(f[-1] + f[-2])
        return f[n]

---

Visualisation

Nous allons maintenant tracer ces valeurs. Nous utilisons le module pyplot de la bibliothèque matplotlib. Les instructions suivantes tracent une série de points définis par deux listes de coordonnées x et y.

import matplotlib.pyplot as plt
plt.plot(xvalues, yvalues, label="name for legend")
plt.legend()
plt.show()

Le graphique doit toujours contenir un titre et une légende, avoir des axes étiquettés, l'échelle de l'axe des y doit être logarithmique, l'échelle de l'axe des x doit commencer à 0. Voici quelques fonctions utiles.

plt.yscale("log")
plt.xlabel("Name")
plt.ylabel("Name")
plt.title("Name")

N'oubliez pas les daltoniens (environ 1 personne sur 20) et choisissez des styles de ligne qui peuvent être différenciés sans utiliser de couleurs. Les styles de ligne peuvent être définis à l'aide des codes suivants : ":", "-", "-." et "--".

plt.plot(xvalues, yvalues, "--", label="Dashed line")

Tout mettre en place

Ajoutez une section principale (main) à votre script Python comme suit :

if __name__ == "__main__":
    # measure execution times
    # plot performance curves

Ajoutez les instructions pour mesurer le temps d'exécution de fibonacci_rec() et fibonacci_iter(), et tracer les courbes de performance.

---

ANSWER ELEMENTS

import matplotlib.pyplot as plt

if __name__ == "__main__":
    l = list(range(1, 20, 1))
    frec = measure(fibonacci_rec, l)
    fiter = measure(fibonacci_iter, l)
    plt.yscale("log")
    plt.title("Execution time of fibonacci functions")
    plt.xlabel("Fibonacci number")
    plt.ylabel("Time in second")

    plt.plot(l, frec, "-.", label="fibonacci_rec")
    plt.plot(l, fiter, "--", label="fibonacci_iter")

    plt.legend()
    plt.show()

---

Méthode alternative

Une méthode alternative pour calculer la suite de Fibonacci consiste à utiliser des matrices. La suite de Fibonacci peut être représentée sous forme matricielle de la manière suivante :

{$$ \left( \begin{matrix} F_1 \\ F_2 \end{matrix} \right) = \left( \begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix} \right) * \left( \begin{matrix} F_0 \\ F_1 \end{matrix} \right) $$}

et en général :

{$$ \left( \begin{matrix} F_n \\ F_{n+1} \end{matrix} \right) = \left( \begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix} \right)^n * \left( \begin{matrix} F_0 \\ F_1 \end{matrix} \right) $$}

Nous utiliserons Numpy pour manipuler les matrices.

Voici quelques étapes simples pour vous aider.

Création de matrices et de vecteurs.

import numpy as np
mat1 = np.array([[5,6],[10,3]])
vec = np.array([[9],[2]])

Accès à un élément.

mat1[1,0]     # => 10

Quelques opérations terme à terme.

mat2 = mat1 + mat1  # => mat2[i,j] = mat1[i, j] + mat1[i,j]
mat3 = mat1 * 2     # => mat3[i,j] = mat1[i,j] * 2
 ...

Quelques opérations matricielles.

mat4 = np.matmul(mat1, mat1)      # matrix product
mat5 = mat1 @ mat1            # abbreviated product syntax
mat6 = np.linalg.matrix_power(mat1, 2)   # mat1**2

Prenez un moment pour vous familiariser avec les opérations matricielles numpy dans un shell Python.

ATTENTION : vous avez peut-être déjà utilisé np.matrix au lieu de np.array. Leur fonctionnement est légèrement différent. En particulier, des opérations comme M ** 2 ou M * vec sont disponibles pour np.matrix et cela semble assez pratique. CEPENDANT, la documentation de numpy est claire, l'objet np.matrix sera bientôt retiré de la bibliothèque, vous devez donc cesser de les utiliser : "Il n'est plus recommandé d'utiliser cette classe, même pour l'algèbre linéaire. Utilisez plutôt des tableaux normaux. La classe pourrait être supprimée dans le futur".

Matrice

Développez une fonction fibonacci_mat(n) permettant de calculer la valeur du n-ième terme de la suite de Fibonacci. Cette fonction doit utiliser des fonctions numpy.

Incluez cette fonction dans votre analyse de performances.

---

ANSWER ELEMENTS

import numpy as np
def fibonacci_mat(n):
    matrix = np.array([[0,1],[1,1]])
    vec = np.array([[0],[1]])
    return (np.linalg.matrix_power(matrix, n-1) @ vec)[1, 0] 

if __name__ == "__main__":
    l = list(range(1, 20, 1))
    frec = measure(fibonacci_rec, l)
    fiter = measure(fibonacci_iter, l)
    fmat = measure(fibonacci_mat, l)
    plt.yscale("log")
    plt.title("Execution time of fibonacci functions")
    plt.xlabel("Fibonacci number")
    plt.ylabel("Time in second")

    plt.plot(l, frec, "-.", label="fibonacci_rec")
    plt.plot(l, fiter, "--", label="fibonacci_iter")
    plt.plot(l, fmat, "-", label="fibonacci_mat")

    plt.legend()
    plt.show()

---

Grandes valeurs de {$n$}

Votre fonction récursive est probablement beaucoup plus lente que les deux autres variantes. Le temps d'exécution semble augmenter de façon exponentielle avec la valeur de {$n$}. Et c'est le cas. Testez vos fonctions avec des valeurs plus importantes de {$n$}. Commencez par range(1, 20, 2), puis range(1, 30, 3) et ainsi de suite jusqu'à ce que le temps d'exécution de fibonacci_rec devienne trop long.

Testez les trois fonctions.

---

ANSWER ELEMENTS

if __name__ == "__main__":
    l = list(range(0, 20, 1))
    l2 = list(range(0, 20, 1)) + list(range(20, 1000, 10))
    frec = measure(fibonacci_rec, l)
    fiter = measure(fibonacci_iter, l2)
    fmat = measure(fibonacci_mat, l2)
    plt.yscale("log")
    plt.title("Execution time of fibonacci functions")
    plt.xlabel("Fibbnacci number")
    plt.ylabel("Time in second")

    plt.plot(l, frec, "-.", label="fibonacci_rec")
    plt.plot(l2, fiter, "--", label="fibonacci_iter")
    plt.plot(l2, fmat, "-", label="fibonacci_mat")

    plt.legend()
    plt.show()

---

Le diable se cache dans les détails

Comparez les résultats de fibonacci_iter(93) et fibonacci_mat(93). D'où vient le problème ?

---

ANSWER ELEMENTS

En Python, les entiers n'ont pas de taille fixe, la seule limite est la mémoire disponible. L'interpréteur gère les débordements et, si nécessaire, un plus grand nombre de bits est utilisé.

Les valeurs entières dans Numpy ont une taille fixe. Par défaut, un entier est une valeur signée sur 64 bits. Les débordements ne sont pas détectés. Conclusion : Numpy est une bibliothèque précieuse mais il faut être prudent quand on l'utilise.

---

Bravo

Vous avez terminé la première partie de ce TD. Vous avez étudié les bases de la visualisation avec matplotlib et, pour la première fois, vous avez touché du doigt la simplicité des calculs matriciels avec numpy. S'il vous reste du temps, vous voudrez peut-être répondre aux questions suivantes.

Questions facultatives

Le problème de la fonction récursive vient du nombre exponentiellement grand d'appels récursifs. La fonction peut être appelée plusieurs fois avec la même valeur en argument. Il est possible de contourner ce problème grâce à la mémoïsation. Le principe de la mémoïsation est de stocker les résultats intermédiaires dans un dictionnaire et de ne faire un appel récursif que si le résultat de l'appel n'est pas encore dans le dictionnaire.

Mémoïsation

Modifiez l'implémentation de votre fonction récursive pour utiliser la mémoïsation et comparez ses performances avec celles des deux autres fonctions.

Vous pouvez en savoir plus sur la mémoïsation ici.

---

ANSWER ELEMENTS

Attention : lorsque l'on donne à un paramètre une valeur par défaut qui est un objet mutable (dictionnaire, liste, ...), cet objet est conservé entre les appels. Donc dans l'exemple suivant, le dictionnaire de mémoïsation n'est complété qu'une unique fois. Cela fausse les comparaisons car le module timeit fait une moyenne sur plusieurs appels.

def fibonacci_rec(n, memo={}):
    if n == 0:
        return 0
    if n == 1:
        return 1
    else:
        if n in memo:
            return memo[n]
        else :
            f = fibonacci_rec(n-1, memo=memo) + fibonacci_rec(n-2, memo=memo)
            memo[n] = f
            return f

Deux solutions sont possibles. Soit utiliser une fonction auxiliaire pour explicitement initialiser le dictionnaire. Soit fixer la valeur par défaut à None.

def fibonacci_rec(n, memo=None):
    if n == 0:
        return 0
    if n == 1:
        return 1
    else:
        if memo is None:
            memo = {}
        if n in memo:
            return memo[n]
        else :
            f = fibonacci_rec(n-1, memo=memo) + fibonacci_rec(n-2, memo=memo)
            memo[n] = f
            return f

---

Très grandes valeurs de {$n$}

Jusqu'à quelle valeur de {$n$} pouvez-vous calculer des résultats ? Cela dépend de votre machine. Augmentez la plage de valeurs de mesure, par exemple de {$0$} à {$20,000$} avec un pas de {$1,000$}. Une des fonctions devrait générer un message d'erreur pour des valeurs trop grandes de {$n$}.

A votre avis

1. Pourquoi fibonacci_mat(n) a un temps d'exécution qui augmente pour de petites valeurs de {$n$} et semble ensuite constant ?
2. Que signifie le message d'erreur fibonacci_rec(n) pour de grandes valeurs de {$n$} ?

---

ANSWER ELEMENTS

1. Pour calculer la puissance d'une matrice {$M^n$}, numpy utilise un algorithme d'exponentiation rapide. Le nombre d'opérations effectuées est proportionnel à {$\log_2(n)$}, pour observer cette augmentation il faudrait des valeurs de {$n$} beaucoup plus importantes. La notion associée à cette question est la complexité des algorithmes, nous y reviendrons lors du cours d'algorithmique ST2.
2. Les appels récursifs sont accumulés sur la pile du processeur. La pile a une taille de mémoire limitée. Lorsque la taille de la pile est dépassée, le processeur ne peut plus continuer son exécution normale, une exception est levée. Vous verrez cette notion de plus près si vous suivez la mention de troisième année "Architecture des systèmes informatiques".

---

Essayez le calcul matriciel à votre façon

Puisque le calcul matriciel avec Numpy génère des débordements, vous pouvez créer votre propre méthode en utilisant des listes d'entiers natifs de Python. Pour cela, vous avez besoin d'une fonction de produit matriciel (basée sur des listes de listes) et d'une fonction d'exponentiation rapide.

---

ANSWER ELEMENTS

def prod(m1, m2):
    result = [[sum(a * b for a, b in zip(X_row, Y_col)) for Y_col in zip(*m2)] for X_row in m1]
    return result

def fast_expo(m, e):
    result = [[1, 0], [0, 1]]
    while e > 0:
        if e % 2 == 1:
            result = prod(result, m)
        m = prod(m, m)
        e = e // 2
    return result

def fibonacci_mat(n):
    matrix = [[0, 1], [1, 1]]
    vec = [[0], [1]]
    e = fast_expo(matrix, n-1)
    return prod(e, vec)[1][0]

---