Objectif

Ce document propose une méthode simple pour implémenter un affichage dynamique sous python. La documentation officielle propose de manière très claire des exemples qui sont rapidement exploitables. Le principe, assez simple, est résumé sur la figure suivante :

demarche_affichage_dynamique

Programme

On prend pour application le tracé dynamique de la réponse d'un système du second ordre quelconque, en réponse à une consigne de type échelon. On appelera \(s(t)\) la grandeur à tracer, \(E_0\) l'amplitude de l'échelon, \(K\) le gain statique de fonction de transfert, \(\omega_0\) sa pulsation propre et \(z\) le coefficient d'amortissement: $$S(p)=\frac{K.E_0}{(1+\dfrac{2.z}{\omega_0}.p+\dfrac{1}{\omega_0^2}.p^2).p}$$

On commence par importer les modules nécessaires:

2
3
import matplotlib.animation as animation
from matplotlib.pylab import *

Le module matplotlib.animation est utilisé pour, comme son nom l'indique, l'affichage dynamique (animé) et matplotlib.pylab pour la création des figures, tracés et fonctions mathématiques.
Suit la définition de la fonction à tracer:

 5
 6
 7
 8
 9
 10
 11
 12
 13
 14
 15
def s(t,w0,z,K,E0):
      if z>1:
            r1=w0*(-z-sqrt(z**2-1))
            r2=w0*(-z+sqrt(z**2-1))
            T1=-1/r1
            T2=-1/r2
            return K*E0/(T1-T2)*(T1*(1-exp(-t/T1))-T2*(1-exp(-t/T2)))
      elif z<1:
            a=sqrt(1-z**2)
            return K*E0*(1-1/a*exp(-z*w0*t)*sin(w0*a*t-atan(-a/z)))
      return K*E0*(1-(1+t*w0)*exp(-w0*t))

On affecte des valeurs aux paramètres caractéristiques de la fonction de transfert, ainsi que l'amplitude de l'échelon. Puis on définit l'instant initial et final pour le tracé, suivi du nombre de points de la courbe:

17
18
19
20
21
22
23
24
w0=1
z=0.1
K=1
E0=1

tmin=0
tmax=30*w0
n=200

On nomme lt la liste du temps (abscisses), et ls la liste des valeurs prises par la fonction s image des élements de lt:

26
27
lt=[tmin+(tmax-tmin)*i/(n-1) for i in range(n)]
ls=[s(i,w0,z,K,E0) for i in lt]

La création de la figure et la définition des tracés se fait grâce au code suivant:

30
31
32
33
34
35
36
37
38
39
fig=figure(num=0,figsize=(12,8))
fig.suptitle("Une animation d'un système du"+"$2^{nd}$"+" ordre.")
a0=subplot2grid((1,1),(0,0))
a0.set_ylim(0,2)
a0.set_xlim(0,tmax)
a0.grid(True)
a0.set_xlabel("Temps (s)")
a0.set_ylabel("Réponse")
p0,=a0.plot([],[],'b',label="Réponse")
a0.legend([p0],[p0.get_label()])

Nous allons ensuite définir une fonction d'initialisation, qui sera appelée au début de l'animation. Dans notre exemple, elle n'est pas nécessaire et est présentée pour illustrer son utilisation. Elle permet d'obtenir un "tracé" (p0) ne contenant aucune donnée (cette opération a déjà été faite lors de la création de p0). A noter que cette fonction doit renvoyer un tuple pour être exploitable par la suite, d'où la virgule à la ligne 43 : return p0, :

41
42
43
def initialisation():
      p0.set_data([],[])
      return p0,

On définit ensuite la fonction animer :

45
46
47
def animer(i):            
      p0.set_data(lt[:i],ls[:i])
      return p0,

Cette fonction sera appelée à intervalle de temps régulier (ou quasiment...). Elle renvoie le tracé p0 (toujours un 1-uplet), dont on a modifié les données via la méthode set_data. L'argument de cette fonction est i, un entier qui sera incrémenté de 1 en 1 automatiquement. A chaque appel, on créé donc un tracé en ajoutant à chaque fois un point sur la sous-figure. Cette opération est réalisée grâce à une copie superficielle partielle des listes lt et ls. La méthode set_data efface les données précédentes pour les remplacer par celle fournies en argument.
Voici enfin la dernière étape:

49
50
anim=animation.FuncAnimation(fig,animer,range(n),interval=1,blit=True)
show()

On créer un objet animation.FuncAnimation qui prend pour arguments:

Finalement, le code complet est donné ci-dessous:

 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
import matplotlib.animation as animation
from matplotlib.pylab import *

def s(t,w0,z,K,E0):
      if z>1:
            r1=w0*(-z-sqrt(z**2-1))
            r2=w0*(-z+sqrt(z**2-1))
            T1=-1/r1
            T2=-1/r2
            return K*E0/(T1-T2)*(T1*(1-exp(-t/T1))-T2*(1-exp(-t/T2)))
      elif z<1:
            a=sqrt(1-z**2)
            return K*E0*(1-1/a*exp(-z*w0*t)*sin(w0*a*t-atan(-a/z)))
      return K*E0*(1-(1+t*w0)*exp(-w0*t))

w0=1
z=0.1
K=1
E0=1

tmin=0
tmax=30*w0
n=200

lt=[tmin+(tmax-tmin)*i/(n-1) for i in range(n)]
ls=[s(i,w0,z,K,E0) for i in lt]


fig=figure(num=0,figsize=(12,8))
fig.suptitle("Une animation d'un système du"+"$2^{nd}$"+" ordre.")
a0=subplot2grid((1,1),(0,0))
a0.set_ylim(0,2)
a0.set_xlim(0,tmax)
a0.grid(True)
a0.set_xlabel("Temps (s)")
a0.set_ylabel("Réponse")
p0,=a0.plot([],[],'b',label="Réponse")
a0.legend([p0],[p0.get_label()])

def initialisation():
      p0.set_data([],[])
      return p0,

def animer(i):            
      p0.set_data(lt[:i],ls[:i])
      return p0,

anim=animation.FuncAnimation(fig,animer,range(n),init_func=initialisation,interval=1,blit=True)
show()

Enfin, il se peut que vous ayez envie de sauvegarder votre animation pour faire un petit film. La fin du programme est modifiée en rajoutant une ligne :

49
50
51
anim=animation.FuncAnimation(fig,animer,range(n),interval=1,blit=True)
anim.save('enregistrement.mp4',fps=25);
show()

Attention, il se peut que vous ayez besoin de "ffmpeg" s'il n'est pas installé. Un très bon tutoriel pour windows est disponible. Pour les autres OS, des ressources doivent exister.
Finalement, pour notre application, voici ce que vous devez obtenir :

Pour aller plus loin

Enrichissons notre animation! Supposons que le modèle décrive un système classique masse avec ressort et amortisseur. Pour ne pas trop alourdir le code, le ressort et l'amortisseur seront modélisés par un câble extensible, dont l'épaisseur varie en fonction de l'alongement. La masse, sera un simple carré. Le support d'accroche du câble, un trait noir assez épais. Le tout sera dans une sous figure . On se propose pour finir de rajouter le protrait de phase dans une troisième sous figure.
Pour réaliser cette "application", on commence comme précédemment par importer les modules nécessaires. On définit également la fonction dérivée ds pour le portrait de phase:

 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import matplotlib.animation as animation
from matplotlib.pylab import *

def s(t,w0,z,K,E0):
      if z>1:
            r1=w0*(-z-sqrt(z**2-1))
            r2=w0*(-z+sqrt(z**2-1))
            T1=-1/r1
            T2=-1/r2
            return K*E0/(T1-T2)*(T1*(1-exp(-t/T1))-T2*(1-exp(-t/T2)))
      elif z<1:
            a=sqrt(1-z**2)
            return K*E0*(1-1/a*exp(-z*w0*t)*sin(w0*a*t-atan(-a/z)))
      return K*E0*(1-(1+t*w0)*exp(-w0*t))

def ds(t,w0,z,k,EO):
      if z>1:
            r1=w0*(-z-sqrt(z**2-1))
            r2=w0*(-z+sqrt(z**2-1))
            T1=-1/r1
            T2=-1/r2
            return K*E0/(T1-T2)*(exp(-t/T1)-exp(-t/T2))
      elif z<1:
            a=sqrt(1-z**2)
            return -K*E0/a*(-z*w0*exp(-z*w0*t)*sin(w0*a*t-atan(-a/z))+w0*a*exp(-z*w0*t)*cos(w0*a*t-atan(-a/z)))
      return K*E0*w0*t*exp(-w0*t)

Pour la modélisation de la masse, on choisit un carré, dont les côtés sont de longueur dim. Le barycentre du carré est, pour \(t=0\), situé à l'origine du repère associé au tracé correspondant. Au cours du temps, seule sa position verticale change, et l'ordonnée de son barycentre vaudra s\((t)\). La fonction masse prendre pour argument une liste y, qui correspond à la liste de la réponse, dim la longueur des côtés et i l'incrément de l'animation. Cette fonction renvoie un tuple de deux listes : les abscisses et les ordonnées des sommets du carré. Il y a 5 élements dans chaque liste pour assuré que le carré soit bien dessiné (rebouclage sur le premier sommet).

29
30
def masse(y,dim,i):
      return [-dim/2,dim/2,dim/2,-dim/2,-dim/2],[y[i]-dim/2,y[i]-dim/2,y[i]+dim/2,y[i]+dim/2,y[i]-dim/2]

Le support du câble est assez simple à réaliser : c'est un segment de longueur dim , situé à linit (longueur initiale du câble) du côté supérieur de la masse pour \(t=0\). \(\forall t>0\), le support se translate brutalement de \(E_0\) pour simuler l'échelon.

32
33
34
35
36
37
def support(i,dim,linit):
      if(i==0):
            h=dim/2+linit
      else:
            h=dim/2+linit+E0
      return [-dim/2,dim/2],[h,h]

Le câble est un simple segment vertical, lié d'une part au support, d'autre part au côté supérieur de la masse:

39
40
41
42
43
44
def cable(y,i,dim,linit):
      if (i==0):
            h=dim/2+linit
      else:
            h=dim/2+linit+E0
      return [0,0],[dim/2+y[i],h]

Enfin, on définit une fonction largeur_cable qui renvoie un flottant traduisant la largeur courante du câble. Pour cette exemple, on choisi une loi linéaire en \(s(t)\), mais d'autres modéles sont biensûr tout aussi valables:

46
47
def largeur_cable(i,y,larg_init,coeff):
      return larg_init+coeff*y[i]

Viennent ensuite les affectations des paramètres de notre modèle:

49
50
51
52
53
54
55
56
57
58
59
larg_init=0.2#largeur initiale du câble
coeff=3#coeff pour amplifier l'effet d'épaississement/amincissement du câble
w0=1
z=0.1
K=1
E0=1
tmin=0
tmax=40*w0
n=200
dim=1#longueur des côtés de la masse
linit=2#longueur initiale du câble.

Et la création des listes du temps, de la réponse et de sa dérivée:

61
62
63
lt=[tmin+(tmax-tmin)*i/(n-1) for i in range(n)]
ls=[s(i,w0,z,K,E0) for i in lt]
lds=[ds(i,w0,z,K,E0) for i in lt]

Comme précédemment, on crée une figure globale:

65
66
fig=figure(num=0,figsize=(12,8))
fig.suptitle("Une animation d'un système du"+"$2^{nd}$"+" ordre.")

Par contre, on désire implémenter 3 sous figures. Plusieurs dispositions sont possibles. On va, par exemple, placer sur la demi figure globale supérieure le tracé de la réponse. La partie inférieure gauche sera réservée au dessin du système MRA et la partie inférieure droite au portrait de phase. Dressons un croquis pour bien configurer les sous tracés:

structure_fig_mra

Une fois la structure mise en place, on peut définir les sous-figures avec leur tracés. Ce qui donne pour la sous-figure de la réponse:

67
68
69
70
71
72
73
74
75
#sous-figure de la réponse temporelle et son objet tracé
a0=subplot2grid((2,2),(0,0),colspan=2)
a0.set_ylim(0,2)
a0.set_xlim(0,tmax)
a0.grid(True)
a0.set_xlabel("Temps (s)")
a0.set_ylabel("Réponse")
p0,=a0.plot([],[],'b',label="Réponse")
a0.legend([p0],[p0.get_label()])

A noter que l'instruction a0=subplot2grid((2,2),(0,0),colspan=2) signifie a0 est une sous figure (de fig) dans une figure globale subdivisée en 2 "lignes" et 2 "colonnes". Sa position est en "(0,0)" et s'étend sur 2 "colonnes".
On créé ensuite la sous figure du système MRA avec:

76
77
78
79
80
81
82
83
84
85
86
87
#idem pour le système masse-ressort-amortisseur
a1=subplot2grid((2,2),(1,0))
a1.axis("equal")
a1.set_ylim(-1,4)
a1.set_xlim(-5,5)
a1.set_ylabel("Un système MRA")
a1.set_xticks([])
a1.set_yticks([])

pmasse,=a1.plot([],[],'r',linewidth=2)#tracé de la masse
psupport,=a1.plot([],[],'k',linewidth=2)#du support
pcable,=a1.plot([],[],'c',linewidth=1)#du câble

Quelques remarques s'imposent:

Pour le portrait de phase, c'est classique:

88
89
90
91
92
93
94
95
96
97
#idem pour le portrait de phase
a2=subplot2grid((2,2),(1,1))
a2.axis("equal")
a2.grid(True)
a2.set_ylim(-3,3)
a2.set_xlim(-0.5,3)
a2.set_xlabel("Réponse")
a2.set_ylabel("Sa dérivée")
p2,=a2.plot([],[],'g',label="Portrait")
a2.legend([p2],[p2.get_label()])

La fin du programme est aussi classique : définition de la fonction d'initialisation, de la fonction animer, et création de l'animation :

99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def initialisation():
      p0.set_data([],[])
      pmasse.set_data([],[])
      psupport.set_data([],[])
      pcable.set_data([],[])
      p2.set_data([],[])
      return p0,pmasse,psupport,pcable,p2

def animer(i):            
      p0.set_data(lt[:i],ls[:i])
      pmasse.set_data(masse(ls,dim,i))
      psupport.set_data(support(i,dim,linit))
      pcable.set_data(cable(ls,i,dim,linit))
      pcable.set_linewidth(largeur_cable(i,ls,larg_init,coeff))
      p2.set_data(ls[:i],lds[:i])
      return p0,pmasse,psupport,pcable,p2

anim=animation.FuncAnimation(fig,animer,range(n),interval=1,init_func=initialisation,blit=True)
show()

En executant le programme, vous devriez obtenir cette animation :

Pour terminer, voici le programme complet:

  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
import matplotlib.animation as animation
from matplotlib.pylab import *

def s(t,w0,z,K,E0):
      if z>1:
            r1=w0*(-z-sqrt(z**2-1))
            r2=w0*(-z+sqrt(z**2-1))
            T1=-1/r1
            T2=-1/r2
            return K*E0/(T1-T2)*(T1*(1-exp(-t/T1))-T2*(1-exp(-t/T2)))
      elif z<1:
            a=sqrt(1-z**2)
            return K*E0*(1-1/a*exp(-z*w0*t)*sin(w0*a*t-atan(-a/z)))
      return K*E0*(1-(1+t*w0)*exp(-w0*t))

def ds(t,w0,z,k,EO):
      if z>1:
            r1=w0*(-z-sqrt(z**2-1))
            r2=w0*(-z+sqrt(z**2-1))
            T1=-1/r1
            T2=-1/r2
            return K*E0/(T1-T2)*(exp(-t/T1)-exp(-t/T2))
      elif z<1:
            a=sqrt(1-z**2)
            return -K*E0/a*(-z*w0*exp(-z*w0*t)*sin(w0*a*t-atan(-a/z))+w0*a*exp(-z*w0*t)*cos(w0*a*t-atan(-a/z)))
      return K*E0*w0*t*exp(-w0*t)

def masse(y,dim,i):
      return [-dim/2,dim/2,dim/2,-dim/2,-dim/2],[y[i]-dim/2,y[i]-dim/2,y[i]+dim/2,y[i]+dim/2,y[i]-dim/2]

def support(i,dim,linit):
      if(i==0):
            h=dim/2+linit
      else:
            h=dim/2+linit+E0
      return [-dim/2,dim/2],[h,h]

def cable(y,i,dim,linit):
      if (i==0):
            h=dim/2+linit
      else:
            h=dim/2+linit+E0
      return [0,0],[dim/2+y[i],h]

def largeur_cable(i,y,larg_init,coeff):
      return larg_init+coeff*y[i]

larg_init=0.2#largeur initial du câble
coeff=3#coeff pour amplifier l'effet d'épaississement/amincissement du câble
w0=1
z=0.1
K=1
E0=1
tmin=0
tmax=40*w0
n=200
dim=1#longueur des côtés de la masse
linit=2#longueur initiale du câbl.

lt=[tmin+(tmax-tmin)*i/(n-1) for i in range(n)]
ls=[s(i,w0,z,K,E0) for i in lt]
lds=[ds(i,w0,z,K,E0) for i in lt]

fig=figure(num=0,figsize=(12,8))
fig.suptitle("Une animation d'un système du"+"$2^{nd}$"+" ordre.")
#sous-figure de la réponse temporelle et son objet tracé
a0=subplot2grid((2,2),(0,0),colspan=2)
a0.set_ylim(0,2)
a0.set_xlim(0,tmax)
a0.grid(True)
a0.set_xlabel("Temps (s)")
a0.set_ylabel("Réponse")
p0,=a0.plot([],[],'b',label="Réponse")
a0.legend([p0],[p0.get_label()])
#idem pour le système masse-ressort-amortisseur
a1=subplot2grid((2,2),(1,0))
a1.axis("equal")
a1.set_ylim(-1,4)
a1.set_xlim(-5,5)
a1.set_ylabel("Un système MRA")
a1.set_xticks([])
a1.set_yticks([])

pmasse,=a1.plot([],[],'r',linewidth=2)#tracé de la masse
psupport,=a1.plot([],[],'k',linewidth=2)#du support
pcable,=a1.plot([],[],'c',linewidth=1)#du câble
#idem pour le portrait de phase
a2=subplot2grid((2,2),(1,1))
a2.axis("equal")
a2.grid(True)
a2.set_ylim(-3,3)
a2.set_xlim(-0.5,3)
a2.set_xlabel("Réponse")
a2.set_ylabel("Sa dérivée")
p2,=a2.plot([],[],'g',label="Portrait")
a2.legend([p2],[p2.get_label()])

def initialisation():
      p0.set_data([],[])
      pmasse.set_data([],[])
      psupport.set_data([],[])
      pcable.set_data([],[])
      p2.set_data([],[])
      return p0,pmasse,psupport,pcable,p2

def animer(i):            
      p0.set_data(lt[:i],ls[:i])
      pmasse.set_data(masse(ls,dim,i))
      psupport.set_data(support(i,dim,linit))
      pcable.set_data(cable(ls,i,dim,linit))
      pcable.set_linewidth(largeur_cable(i,ls,larg_init,coeff))
      p2.set_data(ls[:i],lds[:i])
      return p0,pmasse,psupport,pcable,p2

anim=animation.FuncAnimation(fig,animer,range(n),interval=1,init_func=initialisation,blit=True)
show()

Et la 3D dans tout cela?

Il est bien évidemment possible d'animer des tracés en 3D (courbes, surfaces, nuages de points, etc.) Nous allons, pour clore ce document sur l'affichage dynamique, prendre un simple basé sur les réseaux de neuronnes : on cherche l'extremum d'une fonction \( f(x,y)\) sur un pavé, par exemple \([-3,3]\times[-3,3]\) avec $$f(x,y) = 3.{\left( {1 - x} \right)^2}.{e^{ - {x^2} - {{\left( {y + 1} \right)}^2}}} - 10.\left( {\frac{x}{5} - {x^3} - {y^5}} \right).{e^{ - \left( {{x^2} + {y^2}} \right)}} - \frac{1}{3}.{e^{ - \left( {{{\left( {x + 1} \right)}^2} + {y^2}} \right)}}$$ Cette fonction est reprise d'un exemple de base pour les tracés en 3D sous Python. L'objectif de ce qui va suivre n'est pas d'expliquer le principe d'un algorithme génétique (on trouve énormément de ressources sur le web à ce sujet), mais de montrer quelques instructions bien utiles à l'animation en 3D.
On notera enfin que la partie du code liée à l'algorithme génétique est largement optimisable. Cet aspect ne nous interesse pas ici.
Pour commencer, on importe les modules nécessaires, dont le dernier spécifiquement utile pour cet exemple:

2
3
4
from pylab import *
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D

Toute la partie suivante concerne l'algorithme génétique, elle est donnée pour simple information :

  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#fonction à optimiser
def f(x,y):
      return 3*(1 - x)**2 * exp(-x**2 - (y + 1)**2)\
             - 10*(x/5 - x**3 - y**5)*exp(-x**2 - y**2)\
             - 1./3*exp(-(x + 1)**2 - y**2)

#fonction créant une liste prenant comme argument une liste
#à deux dimensions (point d'abscisse x,y)
# et une fonction, renvoyant la liste des élements image de [x,y] par f
def f_lfech(lp,f):
      return [f(i[0],i[1]) for i in lp]

#fonction de décroissance pour la probabilité de sélection.
#Plus coeff est grand, plus la selection est sévère
def dec_proba(i,nech,coeff):
      return exp(-i/(nech-1)*coeff)
#ordonnance une liste de point [x,y] suivant le tri d'une liste lz
def tri(lp,lz,rev):
      lptri=[lp[i] for i in sorted(range(len(lz)), key=lz.__getitem__,reverse=rev)]
      return lptri

#renvoie un élément trié au hasard dans une liste lp
def tirage(lp):
      k=len(lp)
      if(k==1):
            return lp[0]
      a=randint(0,k-1)
      parentp=lp[a]
      lp.remove(lp[a])
      return parentp

#grille initiale pour l'échantillon
xmin=-3
xmax=3
ymin=-3
ymax=3
#liste des résultats
lres=[]
#valeur absolue de l'amplitude maximale introduite par la phase de mutation,
#sur chaque composant x,y d'un point
delta_gen_mut=0.1
#nombre d'itérations
niter=100
coeff=1
#racine carrée de la taille de l'échantillon. La taille
#doit être un multiple de 4
nech=4
#choix de la recherche de l'extremum de la fonction : False = recherche du minimum,
#True = recherche du maximum
rev=True
lxech = [xmin+(xmax-xmin)*i/(nech-1) for i in range(nech)]
lyech = [ymin+(ymax-ymin)*i/(nech-1) for i in range(nech)]
#listes des abscisses [x,y] des individu de l'échantillon
lpech = [[i,j] for i in lxech for j in lyech]
#ordonnées (z) de l'échantillon
lfech=f_lfech(lpech,f)

for i in range(niter):
      #ordonnacement de l'échantillon suivant le critère voulu
      lpech_sorted=tri(lpech,lfech,rev)
      #création d'une liste de nombre aléatoire dont la valeur décroit statistiquement
      lproba_selec=[rand()*dec_proba(i,nech,coeff) for i in range(nech**2)]
      #Selection de la moitié de l'échantillon.
      lpech_selec=tri(lpech_sorted,lproba_selec,False)[-int(nech**2/2):]
      #nouvel échantillon
      lpnew=[]
      for j in range(int(nech**2/4)):
            #tirage des parents
            parent1=tirage(lpech_selec)
            parent2=tirage(lpech_selec)
            lpnew.append(parent1)
            lpnew.append(parent2)
            #alpha nombre aléatoire donnant la proportion de gènes venant des parents
            alpha=rand()
            #génération des enfants : leurs gènes dépendent des parents, ainsi que
            #d'une part aléatoire (mutation)
            enfant1=[delta_gen_mut*(2*rand()-1)+(1-0)*(alpha*ii+(1-alpha)*jj) for ii,jj in zip(parent1,parent2)]
            enfant2=[delta_gen_mut*(2*rand()-1)+(1-0)*(alpha*ii+(1-alpha)*jj) for ii,jj in zip(parent2,parent1)]
            lpnew.append(enfant1)
            lpnew.append(enfant2)
            lpech=lpnew
            lfech=f_lfech(lpech,f)#image de l'échantillon
      #lres est une liste de listes. lres[i] contient les données à la ième iteration
      #lres[i][0] est la liste contenant les abscisses de l'échantillon
      #lres[i][1] est la liste contenant les ordonnées correspondantes
      lres.append([lpech,lfech])

#remise en forme pour affichage
lxech_res=[[lres[i][0][j][0] for j in range(nech**2)] for i in range(niter)]
lyech_res=[[lres[i][0][j][1] for j in range(nech**2)] for i in range(niter)]
ntrans=40
#contient une liste de liste des x (parents) :lxech_parent[i] est une liste des x parents à la ième iteration
lxech_parent=[]
#la suite est explicite
lxech_enfant=[]
lyech_parent=[]
lyech_enfant=[]
lzech_parent=[]
lzech_enfant=[]
for i in range(niter):
      lparentsx=lxech_res[i][:int(nech**2/2)]
      lparentsy=lyech_res[i][:int(nech**2/2)]
      lenfantsx=lxech_res[i][int(nech**2/2):]
      lenfantsy=lyech_res[i][int(nech**2/2):]
      lcoeffsx=[(ii-jj)/(ntrans-1) for ii,jj in zip(lenfantsx,lparentsx)]
      lcoeffsy=[(ii-jj)/(ntrans-1) for ii,jj in zip(lenfantsy,lparentsy)]
      lparentsz=[f(ii,jj) for ii,jj in zip(lparentsx,lparentsy)]   
      for j in range(ntrans):
            lxech_parent.append(lparentsx)
            lyech_parent.append(lparentsy)
            lzech_parent.append(lparentsz)
            s_lxech_enfant=[lparentsx[k]+lcoeffsx[k]*j for k in range(int(nech**2/2))]
            s_lyech_enfant=[lparentsy[k]+lcoeffsy[k]*j for k in range(int(nech**2/2))]
            lxech_enfant.append(s_lxech_enfant)
            lyech_enfant.append(s_lyech_enfant)
            lzech_enfant.append([f(ii,jj) for ii,jj in zip(s_lxech_enfant,s_lyech_enfant)])

C'est à partir d'ici que les fonctions utiles pour les tracés interviennent. Elles ont pour but de fournir un tracé de la surface représentative de \(f(x,y)\). Cette surface sera immobile au cours de l'animation. Par contre, l'évolution de l'échantillon sera tracé dynamiquement : A chaque itération, les parents seront représentés par des points rouges, et les enfants par des points verts. On utilisera les listes suivantes:

A partir de là, construisons notre figure:

126
127
128
129
130
#création d'une figure
fig = figure(figsize=(10,8))
fig.suptitle("A la montagne!")
#a0 est une sous figure en 3D
a0 = subplot2grid((1,1),(0,0), projection='3d')

Pour indiquer à Python qu'on va tracer une figure en 3D, on utilise projection='3d'.
On définit ensuite:

132
133
134
135
136
137
138
139
140
141
142
143
144
#nombre de points en x et en y pour la surface à tracer
nplot=150
#limites du tracé
xminplot=-3
xmaxplot=3
yminplot=-3
ymaxplot=3
#mise en forme pour la figure en 3D
lxplot = [xminplot+(xmaxplot-xminplot)*i/(nplot-1) for i in range(nplot)]
lyplot = [yminplot+(ymaxplot-yminplot)*i/(nplot-1) for i in range(nplot)]
lxplotgrid = [lxplot for i in range(nplot)]
lyplotgrid = [[lyplot[i]]*nplot for i in range(nplot)]
lzplotgrid=[[f(k,l) for k,l in zip(i,j)] for i,j in zip(lxplotgrid,lyplotgrid)]

Toutes ces lignes servent à définir les données utiles pour le tracé de la surface. A noté qu'en utilisant meshgrid de numpy, la syntaxe est bien plus légère, mais rien n'est impossible sans numpy! Ensuite, on définit les labels pour chaque direction, ainsi que le point de vue : c'est la direction suivant laquelle la "caméra" pointe vers la surface. Cette direction est classiquement définie par les angles d'azimut et d'élévation:

145
146
147
148
149
150
#mise en place des label pour a0
a0.set_xlabel("x")
a0.set_ylabel("y")
a0.set_zlabel("z")
#mise en place du point de vue pour la figure 3d
a0.view_init(elev=40, azim=-130)

On peut ensuite tracer la surface avec la méthode plot_surface:

152
a0.plot_surface(lxplotgrid, lyplotgrid, lzplotgrid, linewidth=0,rstride=1,cstride=1,cmap=cm.jet)

On définit ensuite les tracés : un pour les parents et un pour les enfants. Ils sont créés vides:

153
154
155
156
157
158
#creation des tracés
#pour les parents
p0,=a0.plot([],[],[])
#et les enfants
p1,=a0.plot([],[],[])
#initialisation des tracés

Comme pour les exemples précédents, on retrouve la fonction d'initialisation :

159
160
161
162
163
164
def initialisation():
      p0.set_data([],[])
      p0.set_3d_properties([])
      p1.set_data([],[])
      p1.set_3d_properties([])
      return p0,p1

On notera que pour associer des données aux tracés, on utilise set_data qui prendra comme arguments la grille des x et des y et set_3d_properties qui prendra comme argument la grille des z
On retrouve ensuite la fonction qui met à jour les données pour l'animation:

168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
def animer(i):
      #mise en place des donnees (x,y) de la ième iteraion
      p0.set_data(lxech_parent[i],lyech_parent[i])
      #et en ordonnées
      p0.set_3d_properties(lzech_parent[i])
      #styles
      p0.set_marker("o")
      p0.set_ls(".")
      p0.set_ms(10)
      p0.set_mfc([1,0,0])
      p0.set_alpha(.5)
      #idem pour les enfants
      p1.set_data(lxech_enfant[i],lyech_enfant[i])
      p1.set_3d_properties(lzech_enfant[i])
      p1.set_marker("o")
      p1.set_ls(".")
      p1.set_ms(10)
      p1.set_mfc([0,1,0])
      p1.set_alpha(.5)
      return p0,p1

On termine enfin le programme classiquement par:

189
190
191
ani = animation.FuncAnimation(fig, animer, range(niter*ntrans),init_func=initialisation,
                                interval=10, blit=True, repeat=True)
show()

Voici ce qu'on obtient pour la recherche du maximum de la fonction :

Pour information, voici le code complet :

  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
from pylab import *
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D




#fonction à optimiser
def f(x,y):
      return 3*(1 - x)**2 * exp(-x**2 - (y + 1)**2)\
             - 10*(x/5 - x**3 - y**5)*exp(-x**2 - y**2)\
             - 1./3*exp(-(x + 1)**2 - y**2)

#fonction créant une liste prenant comme argument une liste
#à deux dimensions (point d'abscisse x,y)
# et une fonction, renvoyant la liste des élements image de [x,y] par f
def f_lfech(lp,f):
      return [f(i[0],i[1]) for i in lp]

#fonction de décroissance pour la probabilité de sélection.
#Plus coeff est grand, plus la selection est sévère
def dec_proba(i,nech,coeff):
      return exp(-i/(nech-1)*coeff)
#ordonnance une liste de point [x,y] suivant le tri d'une liste lz
def tri(lp,lz,rev):
      lptri=[lp[i] for i in sorted(range(len(lz)), key=lz.__getitem__,reverse=rev)]
      return lptri

#renvoie un élément trié au hasard dans une liste lp
def tirage(lp):
      k=len(lp)
      if(k==1):
            return lp[0]
      a=randint(0,k-1)
      parentp=lp[a]
      lp.remove(lp[a])
      return parentp

#grille initiale pour l'échantillon
xmin=-3
xmax=3
ymin=-3
ymax=3
#liste des résultats
lres=[]
#valeur absolue de l'amplitude maximale introduite par la phase de mutation,
#sur chaque composant x,y d'un point
delta_gen_mut=0.1
#nombre d'itérations
niter=100
coeff=1
#racine carrée de la taille de l'échantillon. La taille
#doit être un multiple de 4
nech=4
#choix de la recherche de l'extremum de la fonction : False = recherche du minimum,
#True = recherche du maximum
rev=True
lxech = [xmin+(xmax-xmin)*i/(nech-1) for i in range(nech)]
lyech = [ymin+(ymax-ymin)*i/(nech-1) for i in range(nech)]
#listes des abscisses [x,y] des individu de l'échantillon
lpech = [[i,j] for i in lxech for j in lyech]
#ordonnées (z) de l'échantillon
lfech=f_lfech(lpech,f)

for i in range(niter):
      #ordonnacement de l'échantillon suivant le critère voulu
      lpech_sorted=tri(lpech,lfech,rev)
      #création d'une liste de nombre aléatoire dont la valeur décroit statistiquement
      lproba_selec=[rand()*dec_proba(i,nech,coeff) for i in range(nech**2)]
      #Selection de la moitié de l'échantillon.
      lpech_selec=tri(lpech_sorted,lproba_selec,False)[-int(nech**2/2):]
      #nouvel échantillon
      lpnew=[]
      for j in range(int(nech**2/4)):
            #tirage des parents
            parent1=tirage(lpech_selec)
            parent2=tirage(lpech_selec)
            lpnew.append(parent1)
            lpnew.append(parent2)
            #alpha nombre aléatoire donnant la proportion de gènes venant des parents
            alpha=rand()
            #génération des enfants : leurs gènes dépendent des parents, ainsi que
            #d'une part aléatoire (mutation)
            enfant1=[delta_gen_mut*(2*rand()-1)+(1-0)*(alpha*ii+(1-alpha)*jj) for ii,jj in zip(parent1,parent2)]
            enfant2=[delta_gen_mut*(2*rand()-1)+(1-0)*(alpha*ii+(1-alpha)*jj) for ii,jj in zip(parent2,parent1)]
            lpnew.append(enfant1)
            lpnew.append(enfant2)
            lpech=lpnew
            lfech=f_lfech(lpech,f)#image de l'échantillon
      #lres est une liste de listes. lres[i] contient les données à la ième iteration
      #lres[i][0] est la liste contenant les abscisses de l'échantillon
      #lres[i][1] est la liste contenant les ordonnées correspondantes
      lres.append([lpech,lfech])

#remise en forme pour affichage
lxech_res=[[lres[i][0][j][0] for j in range(nech**2)] for i in range(niter)]
lyech_res=[[lres[i][0][j][1] for j in range(nech**2)] for i in range(niter)]
ntrans=40
#contient une liste de liste des x (parents) :lxech_parent[i] est une liste des x parents à la ième iteration
lxech_parent=[]
#la suite est explicite
lxech_enfant=[]
lyech_parent=[]
lyech_enfant=[]
lzech_parent=[]
lzech_enfant=[]
for i in range(niter):
      lparentsx=lxech_res[i][:int(nech**2/2)]
      lparentsy=lyech_res[i][:int(nech**2/2)]
      lenfantsx=lxech_res[i][int(nech**2/2):]
      lenfantsy=lyech_res[i][int(nech**2/2):]
      lcoeffsx=[(ii-jj)/(ntrans-1) for ii,jj in zip(lenfantsx,lparentsx)]
      lcoeffsy=[(ii-jj)/(ntrans-1) for ii,jj in zip(lenfantsy,lparentsy)]
      lparentsz=[f(ii,jj) for ii,jj in zip(lparentsx,lparentsy)]   
      for j in range(ntrans):
            lxech_parent.append(lparentsx)
            lyech_parent.append(lparentsy)
            lzech_parent.append(lparentsz)
            s_lxech_enfant=[lparentsx[k]+lcoeffsx[k]*j for k in range(int(nech**2/2))]
            s_lyech_enfant=[lparentsy[k]+lcoeffsy[k]*j for k in range(int(nech**2/2))]
            lxech_enfant.append(s_lxech_enfant)
            lyech_enfant.append(s_lyech_enfant)
            lzech_enfant.append([f(ii,jj) for ii,jj in zip(s_lxech_enfant,s_lyech_enfant)])

#création d'une figure
fig = figure(figsize=(10,8))
fig.suptitle("A la montagne!")
#a0 est une sous figure en 3D
a0 = subplot2grid((1,1),(0,0), projection='3d')

#nombre de points en x et en y pour la surface à tracer
nplot=150
#limites du tracé
xminplot=-3
xmaxplot=3
yminplot=-3
ymaxplot=3
#mise en forme pour la figure en 3D
lxplot = [xminplot+(xmaxplot-xminplot)*i/(nplot-1) for i in range(nplot)]
lyplot = [yminplot+(ymaxplot-yminplot)*i/(nplot-1) for i in range(nplot)]
lxplotgrid = [lxplot for i in range(nplot)]
lyplotgrid = [[lyplot[i]]*nplot for i in range(nplot)]
lzplotgrid=[[f(k,l) for k,l in zip(i,j)] for i,j in zip(lxplotgrid,lyplotgrid)]
#mise en place des label pour a0
a0.set_xlabel("x")
a0.set_ylabel("y")
a0.set_zlabel("z")
#mise en place du point de vue pour la figure 3d
a0.view_init(elev=40, azim=-130)
#tracé de la courbe sur la sous figure
a0.plot_surface(lxplotgrid, lyplotgrid, lzplotgrid, linewidth=0,rstride=1,cstride=1,cmap=cm.jet)
#creation des tracés
#pour les parents
p0,=a0.plot([],[],[])
#et les enfants
p1,=a0.plot([],[],[])
#initialisation des tracés
def initialisation():
      p0.set_data([],[])
      p0.set_3d_properties([])
      p1.set_data([],[])
      p1.set_3d_properties([])
      return p0,p1

#fonction de rafraichissement

def animer(i):
      #mise en place des donnees (x,y) de la ième iteraion
      p0.set_data(lxech_parent[i],lyech_parent[i])
      #et en ordonnées
      p0.set_3d_properties(lzech_parent[i])
      #styles
      p0.set_marker("o")
      p0.set_ls(".")
      p0.set_ms(10)
      p0.set_mfc([1,0,0])
      p0.set_alpha(.5)
      #idem pour les enfants
      p1.set_data(lxech_enfant[i],lyech_enfant[i])
      p1.set_3d_properties(lzech_enfant[i])
      p1.set_marker("o")
      p1.set_ls(".")
      p1.set_ms(10)
      p1.set_mfc([0,1,0])
      p1.set_alpha(.5)
      return p0,p1

ani = animation.FuncAnimation(fig, animer, range(niter*ntrans),init_func=initialisation,
                                interval=10, blit=True, repeat=True)
show()