Bases de programmation scientifique avec le langage Python.
Frédéric Faure
Université Grenoble Alpes, France
(version: 15 avril 2022)

Table des matières

Introduction
Partie I Bases du langage python
Chapitre 1 Objets de bases
Chapitre 2 Objets plus élaborés
Chapitre 3 Lecture et écriture de fichiers
Chapitre 4 Résoudre une équation différentielle ordinaire (EDO) avec Scipy
Chapitre 5 Micro-projets
Partie II Compléments du langage python
Chapitre 6 Classes et Objets
Chapitre 7 Analyse de données avec numpy
Chapitre 8 Quelques modules utiles
Chapitre 9 Librairie PyRoot pour le graphisme et l'analyse de données
Chapitre 10 Librairie opencv pour l'analyse d'images
Chapter 11 Python virtual environments
Chapitre 12 Librairie Caffe pour le deep Learning (DL)
Chapitre 13 Deep learning
Chapitre 14 Convertir python en javascript

Attention: Pour lire ce document, utiliser Firefox qui donne un bon rendu des équations.

Introduction

0.1 Déroulement des 10 séances en Licence L3 de physique, 2018-2019


0.2 Des références générales sur python

Voici des références pour compléter ce didacticiel (qui est d'ailleurs construit à partir de ces références).
Sur internet:
Des MOOC
Livres (téléchargeables en version pdf):

0.3 Installation de Python sur votre ordinateur

On conseille d’installer python3 et l'environnement de programmation spyder.
ainsi que les librairies scientifiques.
Par exemple utiliser la distribution Anaconda

0.3.0.1 Sur Linux, Ubuntu (recommandé)

Commandes dans un terminal:
sudo apt install python3-pip python3-dev
sudo apt install spyder

0.3.0.2 Sur Windows

Voir WinPython

0.3.0.3 Sur Mac-Os

0.3.0.4 Sur Android

0.3.0.5 En général sur un navigateur

0.4 Utilisation de python avec le logiciel « spyder »

Lancement du logiciel spyder:
dans un terminal, lancer: spyder (ou spyder3)
Il apparaît un environnement de programmation pour le langage python appelé « spyder ».
premiers Tests à faire:
 
Programmes avec graphisme:
dans ce cas il est conseillé dans spyder dans le menu Outils/Préférences/ConsoleIPython/Graphiques/Sortie : de choisir le mode automatique.
(et relancer la fenetre IPython)

0.4.1 Lancement de votre programme python dans une fenêtre de commandes

Voici comment éxecuter votre programme temp.py ci-dessus dans une fenêtre de commandes.
  1. Ouvrir un terminal (fenêtre de commandes)
  2. Aller dans le répertoire qui contient votre programme temp.py
  3. Ecrire:
    python3 temp.py
 
Dans la suite de ce didacticiel, vous pourrez écrire les lignes de codes proposées en mode interactif ou en mode programme. Pour recopier le(s) code(s) proposé(s), vous pouvez faire « copier/coller ».

0.5 (supplément) Utilisation de python avec le logiciel « Jupyter »

Si vous avez choisit de travailler avec « spyder » comme ci-dessous, vous pouvez ignorer cette Section.
Références: sur datacamp
Lancement du logiciel Jupyter
  1. dans un terminal, lancer: anaconda-navigator
  2. Cliquer sur « Launch JupyterLab ». Il apparaît un environnement de programmation pour le langage python appelé « JupyterLab ».
  3. Cliquer sur l'icone « NoteBook/Python3 »
Mode commandes: si le bord est bleu.
Mode éditeur: si le bord est vert.

0.6 (supplément) Utilisation de python en ligne

avec le site Colaboratory de google:
https://colab.research.google.com/notebooks/welcome.ipynb

0.7 Utilisation de commandes dans un terminal

Si vous avez une question sur l'utilisation de ubuntu, par exemple comment ouvrir et utiliser un terminal, écrire dans google: “ubuntu terminal”.
Dans ubuntu, ouvrir une fenêtre de commandes (ou terminal). Dans cette fenêtre de commandes, vous pouvez écrire des commandes à l'ordinateur. Ces commandes sont dans le langage UNIX, système avec lequel fonctionne l'ordinateur. Cette section vous explique les principales commandes de UNIX qu'il faut connaître: les commandes qui permettent de gérer l'emplacement et le nom de vos fichiers (qui seront par exemples vos programmes en python).
Il est important de savoir que les fichiers sont classés dans l'ordinateur (sur le disque magnétique) selon une structure arborescente. Chaque noeud de l'arborescence s'appelle un répertoire ou dossier ou directory (en anglais). Dans ces répertoires se trouvent des fichiers (ou file en anglais).
Par exemple /h/moore/u4/phyliphrec/aaal/toto.py signifie que le fichier toto.py se trouve dans le répertoire aaal qui se trouve lui même dans le répertoire phyliphrec, qui se trouve lui même dans ...

0.7.1 Les répertoires

image: 18_home_faure_enseignement_informatique_c++_cours_c++_cours_1_terminal.png

0.7.2 Exercices sur les répertoires

A l'aide des commandes expliquées ci-dessus:
  1. Placez vous dans votre répertoire d'origine: cd
  2. Vérifiez que vous y êtes bien: pwd
  3. Regardez la liste des fichiers par: ls , puis avec: ls -als
  4. Créer un répertoire du nom: essai Vérifier son existence (avec ls).
  5. Aller dans ce répertoire: cd essai Vérifier que vous y êtes (avec pwd).
  6. Revenir dans le répertoire d'origine par cd .. Vérifiez que vous y etes. Détruire le répertoire essai

0.7.3 Les fichiers (en exercice)

0.7.4 La documentation sur les commandes unix

Partie I Bases du langage python

Chapitre 1 Objets de bases

Référence sur les types.

1.1 Nombres entiers et nombres à virgule. Entrée et sortie à l'écran.

Pour créer une variable, il suffit de lui affecter une valeur. Le texte après le symbole # est un commentaire et est ignoré par l'ordinateur. Ecrire:
x = 1 # met le chiffre 1 dans la variable x
print (x)
résultat:
1
On voit ici une différence importante avec d'autres langages : il n'est pas nécessaire de déclarer a priori le type de la variable x (entier 'int' (integer), nombre flottant 'float', chaine de caractère 'string',...). Pourtant ce type existe bien, si vous écrivez ensuite :
print (type(x))
résultat:
<type 'int'>
Python a donc automatiquement décidé que x était un entier (int = integer). Ré-essayez avec d'autres initialisations : x=1.0 (nombre à virgule flottante ou float), ou x = 'bonjour' (chaine de caractère ou string),...
On peut affecter plusieurs variables à la fois. Ecrire:
x,y = 1, 2
print (x,y)
résultat:
1 2
ou écrire:
x = y = 1
print (x)
print (y)
résultat:
1
1
Remarque: dans l'exemple précédent, pour que les deux résultats soient sur la même ligne (sans saut de ligne), il faut rajouter ,end=' ' à la fin de la première ligne. Ecrire:
x = y = 1
print (x,end=' ')
print (y)
résultat:
1 1
Pour demander à l'utilisateur le contenu d'une variable, écrire:
x = input('x=? ') # demande un nombre et le met dans la variable x
print (x)
print (type(x))
y = int(x) # convertit x qui est de type 'str' en type 'int'
print (y)
print (type(y))
résultat:
x=? 2
2
<type 'str'>
2
<type 'int'>
Remarque 1.1.1. On observe dans l'exemple précédent qu'une variable obtenue par l'instruction input() est de type 'string'. Si on veut la convertir en nombre entier il faut appliquer la fonction int().
Attention : majuscules et minuscules sont différenciées par python, i.e. ma_variable, MA_VARIABLE et Ma_Variable sont trois variables distinctes. Ecrire:
ma_variable = 1
Ma_variable = 2
print (ma_variable)
 
résultat:
1
Pour que la console oublie la valeur de x (et de toutes les variables) il faut écrire: %reset

1.2 Chaînes de caractères

Une chaine de caractère est un tableau contenant des caractères. On peut ainsi construire des mots et des phrases et les manipuler. La numérotation des « cases » commence à zéro.
Ecrire:
nom='Vincent'
print (nom)
print (nom[2:6])    # extrait les caractères 2 (inclus) à 6 (exclu),
Résultat:
Vincent
ncen
Remarque 1.2.1. - le premier caractère est donc nom[0] et contient 'V'.
Il est possible d'additionner des chaînes de caractères. Ecrire:
print ('Programmation' + ' Scientifique')
résultat:
Programmation Scientifique
On peut obtenir la longueur d'une chaîne de caractère avec la fonction len(). Ecrire:
len('toto')
résultat:
4

1.2.1 Formatage

Référence:
format.
Il est possible d'écrire une chaîne de caractère en insérant des codes spécifiques dans la chaîne (de la même manière qu'avec la fonction printf en langage C ): (%s pour une chaîne de caractères, %d pour un entier, %f pour un nombre à virgule, et optionellement le nombre de caractères total et après la virgule: %4.2f veut dire 'nombre représenté sur 4 caractères, avec 2 chiffres après la virgule'), suivie de % avec entre parenthèses toutes les variables dont les valeurs sont à insérer. Ecrire:
age=5
nom="Vincent"
taille=1.73
maChaine="%s a %d ans et mesure %4.2f m" % (nom,age,taille)
print (maChaine)
résultat:
Vincent a 5 ans et mesure 1.73 m
Pour avoir la liste des fonctions utilisables avec des chaînes de caractères, écrire:
help(str)

1.3 Conversion de variables entre différents types

On peut convertir une variable d'un type à un autre. Ecrire:
x=float(1) # nombre réel à partir d'un entier
print (type(x))
résultat:
<type 'float'>
Ecrire:
x=float('1.0') # nombre réel à partir d'une chaîne de caractère
print (type(x))
résultat:
<type 'float'>
Ecrire:
print (x)
résultat:
1.0
Exercice 1.3.1. Ecrire un programme qui demande le nom et l'âge d'une personne. En retour elle affiche l'année de naissance. Par exemple:
ton nom=? blaise
ton age=? 25
Salut blaise , tu es né en 1993 ?

1.4 Opérations mathématiques

1.4.1 opérations de base

1.4.1.1 Opération de division et puissance

écrire:
print (7/2, 7//2, 3**2)
résultat:
3.5 3 9
Noter que la première opération / a donné la division décimale et la deuxième opération // a donné la division Euclidienne (le résultat est un entier).

1.4.1.2 Modulo

Le « modulo »: n%q est le reste de la division euclidienne de n par q :
n:
-2
-1
0
1
2
3
4
5
n%3
1
2
0
1
2
0
1
2
Ecrire:
print (4%3, -1%3) # 4 modulo 3, (-1) modulo 3
print (7.5%3.0, -2.5%3.0) # 7.5 modulo 3.0, (-2.5) modulo 3.0
résultat:
1 2
1.5 0.5

1.4.1.3 Fonctions mathématiques

Pour calculer une racine carrée: (le texte après le signe # est un commentaire)
from math import * # cela importe le module math
sqrt(2.)
résultat:
1.4142135623730951
Pour avoir la liste des fonctions mathématiques disponibles écrire
help()
math
Remarque : pour importer un module (tel que math), il y a deux manières de le faire, avec (a) “import math” ou (b) “from math import *”. Dans le cas (a), les fonctions du module math s'utiliseront en les faisant précéder de “math.” (exemple : “math.sin(math.pi/3.0)”). Dans le cas (b), il n'est pas nécessaire d'utiliser le préfixe, i.e. python comprendra directement “sin(pi/3.0)”. A priori, la première notation n'est à utiliser que si il y a un risque de conflit de noms entre fonctions, si par exemple vous souhaitez redéfinir par ailleurs la fonction sin().

1.4.2 Nombre aléatoire

Documentation sur random.
Pour choisir un nombre à virgule au hasard entre 0 et 1, écrire:
import random
x = random.random()
print (x)
résultat:
0.9565461152605507
Pour choisir au hasard un nombre entier entre 1 et 6 compris, écrire:
import random
x = random.randint(1, 6)
print (x)
Résultat:
3
Exercice 1.4.1. Ecrire un programme qui demande un nombre entier x à l'utilisateur et en retour affiche un entier choisit au hasard entre 0 et x. On souhaite le résultat:
x=? 3
2
Solution:
 
import random
x = input('x=? ')
print random.randint(0, x)

1.4.3 Nombres complexes

Il est possible de manipuler des nombres complexes en Python, le nombre complexe i étant noté 1j . Ecrire:
a=2+1j
type(a)
résultat
<type 'complex'>
Ecrire
b=3.5+4j
print (a+b)
résultat
(5.5+5j)
Ecrire
print (a*b)
résultat:
(3+11.5j)
Pour utiliser des fonctions mathématiques sur des complexes, il faut importer le module cmath en plus du module math. Par exemple pour calculer e i π qui est égal à 1 , écrire:
from math import *
from cmath import *
print (exp(1j*pi))
résultat:
(-1+1.2246467991473532e-16j)

1.4.4 (supplément) Nombre en base 2, binaire.

Les 2 symboles utilisés en base 2 sont: 0,1 . Pour indiquer qu'un nombre est écrit en base 2, on précède son écriture du symbole 0b.
Ecrire:
#... Declaration et Affichage en base deux:
A=0b1001 # declaration en base 2 (prefixe 0b)
print 'A=' , A , ' en binaire = ' , bin(A)
 
#.. decalages de bits
b= A<<4 # equivalent: b= A * (2**4)
print ' b =', b, ' = ', bin(b)
c = b>>3 # equivalent: c= b / (2**3)
print ' c =', c, ' = ', bin(c)
 
#.. operations bits a bits
B = 0b1;
d = A | B # ou inclusif bit a bit cad 1|1 donne 1
e = A ^ B # ou exclusif bit a bit cad 1^1 donne 0
f = A & B # et bit a bit
g = ~A # non A, inversion des bits
 
print ' B=', B, ' = ', bin(B)
print ' A|B = ', d, ' = ', bin(d)
print ' A^B = ', e, ' = ', bin(e)
print ' A&B = ', f, ' = ', bin(f)
print ' ~A = ', g, ' = ', bin(g)
Résultats:
A= 9 en binaire = 0b1001
b = 144 = 0b10010000
c = 18 = 0b10010
B= 1 = 0b1
A|B = 9 = 0b1001
A^B = 8 = 0b1000
A&B = 1 = 0b1
~A = -10 = -0b1010
Référence.

1.4.5 (supplément) Nombre en base 16 (hexadécimale)

Les 16 symboles utilisés en base 16 sont: 0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,f . Pour indiquer qu'un nombre est écrit en base 16, on précède son écriture du symbole 0x. Ecrire:
nombre = 0xb # déclaration en base 16
print (nombre) # affichage en base 10
Résultat (écriture en base 10):
11

1.5 Listes

Un type particulièrement intéressant est la liste qui permet de construire et manipuler une liste d'objets divers. Référence sur le type list.
Ecrire:
ma_liste = [1,4, 'ABC']
print (ma_liste)
résultat:
[1, 4, 'ABC']
Remarque 1.5.1. La numérotation des objets de la liste commence à 0 (i.e. si il y a n éléments dans ma_liste, ils sont accessibles par ma_liste[0]...ma_liste[n-1]).
Ecrire:
ma_liste = [1,4, 'ABC']
print (ma_liste[2])
résultat:
ABC
Ecrire:
ma_liste = [1,4, 'ABC']
print (type(ma_liste))
print (type(ma_liste[2]))
résultat:
<type 'list'>
<type 'str'>

1.5.1 Opérations sur les listes

Il est possible de modifier une liste existante. Ecrire (le texte après le signe # est un commentaire)
L = [1,3,7,5]
L[1]=42 # le 2ème élément de MaListe vaut maintenant 42
print ('L1=', L)
 
L.append(45) # ajoute l'objet “45” à la fin de la liste
print ('L2=', L)
 
L.insert(2, 33) # insère (rajoute) l'élément 33 à l'indice 2 dans MaListe
print ('L3=', L)
 
del L[1] # enleve element de la case 1
print ('L3a=', L)
 
res = L.pop(1) # enleve element de la case 1 et le donne
print ('L3b=', L, ' res=',res)
 
L.remove(45) # enleve le premiere element trouve egal a 45
print ('L3c=', L)
 
L.reverse() # inverse l'ordre des éléments de la liste
print ('L4=', L)
 
L.sort() # trie les éléments de la liste par ordre croissant
print ('L5=', L)
 
L2 = [] # cree une liste vide
print ('L6=', L2)
 
L3 = list(L) # copie dans une autre liste
L3[0] = 0
L4 = L # ce n'est pas une vraie copie: liste4 est le meme objet que MaListe
L4[0] = -1
print ('L7=', L3)
print ('L8=', L)
résultat:
L1= [1, 42, 7, 5]
L2= [1, 42, 7, 5, 45]
L3= [1, 42, 33, 7, 5, 45]
L3a= [1, 33, 7, 5, 45]
L3b= [1, 7, 5, 45] res= 33
L3c= [1, 7, 5]
L4= [5, 7, 1]
L5= [1, 5, 7]
L6= []
L7= [0, 5, 7]
L8= [-1, 5, 7]
On peut obtenir la longueur d'une liste (le nombre d'objets dans la liste) avec la fonction len():
MaListe = [1,3,7,5]
len(MaListe)
résultat:
4
Il est possible de concaténer deux listes avec l'opérateur + :
liste1 = [1,'ab',3.5]
liste2 = ['toto',57]
print (liste1 + liste2)
résultat:
[1, 'ab', 3.5, 'toto', 57]
On peut également supprimer un élément d'une liste avec la fonction del() :
MaListe=[1, 'ab', 3.5, 'toto', 57]
print (MaListe)
del(MaListe[1])
print (MaListe)
résultat:
[1, 'ab', 3.5, 'toto', 57]
[1, 3.5, 'toto', 57]
Pour créer des listes particulières de nombres:
range(2,7) # liste des entiers consécutifs de 2 (inclu) à 7 (exclu)
résultat:
[2, 3, 4, 5, 6]
Si on commence de 0 il est plus simple d'écrire:
range(7) # liste de 0 (inclu) à 7 (exclu)
résultat:
[0, 1, 2, 3, 4, 5, 6]
Les listes sont des objets particulièrement utiles. Nous verrons que les boucles utilisent souvent une itération sur des éléments d'une liste. La liste complète des fonctions utilisables pour une liste peut être obtenue avec l'aide help(list).

1.5.2 Listes et hasard

Pour choisir au hasard un élément parmi une liste donnée, écrire:
import random
x = random.choice(['a', 'b', 'k', 'p', 'i', 'w', 'z'])
print (x)
Résultat:
'k'
Pour mélanger au hasard une liste donnée écrire:
import random
liste = ['a', 'b', 'k', 'p', 'i', 'w', 'z']
random.shuffle(liste)
print (liste)
résultat:
['p', 'k', 'w', 'z', 'i', 'b', 'a']
Exercice 1.5.2. Ecrire un programme qui choisit au hasard un élément parmi une liste donnée, comme le programme ci-dessus, mais en utilisant les fonction len() et random.randint().
Solution:
 
import random
liste = ['a', 'b', 'k', 'p', 'i', 'w', 'z']
print (liste[random.randint(0,len(liste)-1)])

1.6 Boucles et opérateurs

1.6.1 opérateurs ==, and, or , not

L'opérateur == permet de comparer deux objets. Ecrire:
1==2     # est-ce que 1 est égal à 2 ?
1==1
résultat:
False
True
Les autres opérateurs de comparaison sont <, >, <=, >=,!= (ce dernier signifie différent). Par convention, l'entier 0 correspond à False et les entiers non nuls correspondent à True.
L'opérateur not(x) renvoie la négation (au sens booléen, vrai/faux) d'une expression :
not(True)
not(False)
not(0)
résultat:
False
True
True
and et or permettent de réaliser les opérations logiques usuelles (attention il ne s'agit pas d'opérations bit à bit, il ne faut donc utiliser and et or que pour des tests booléens vrai/faux):
1 and 0
1 and 1
1 or 0
résultat:
0
1
1

1.6.2 L'instruction if...elif..else

Remarque 1.6.1. “elif” signifie “else if”.
On peut exécuter une instruction conditionnelle avec if. Écrire:
x=1
print ('x=',x)
if x==2:                 # est-ce que x est égal à 2 ?
    print ('x est égal à 2')  # notez les espaces en début de ligne (indentation)
elif x==3:
    print ('x est égal à 3')
else:
    print ("x n'est pas égal à 2")
    print ('ni à 3')
résultat:
x=1
x n'est pas égal à 2
ni à 3  
Les espaces en début de ligne (indentation) jouent un rôle important : tant que l'espace reste le même en début de ligne, on reste dans le même “niveau” d'éxécution, de même que dans d'autres langages on utilise des accolades {} pour délimiter les instructions à réaliser. On peut ainsi réaliser des instructions imbriquées en augmentant le niveau d'indentation :
x=12
if (x%2)==0:            # est-ce que x modulo 2 est égal à 0 ?
    if (x%3)==0:
        if (x%7)==0:
            print ('x est multiple de 2,3 et 7 !')
        else:
            print ('x est multiple de 2 et 3, mais pas de 7...')
    else:
        print ('x est multiple de 2 mais pas de 3')
else:
    print ('x n'est pas multiple de 2')
résultat:
x est multiple de 2 et 3, mais pas de 7...

1.6.3 La boucle while

On peut réaliser une série d'instructions tant qu'une condition est réalisée :
x=7
while x>=1:
    print (x,'>=1 ! On continue...')
    x = x-1                 # on pourrait aussi écrire x -= 1
résultat:
7 >=1 ! On continue...
6 >=1 ! On continue...
5 >=1 ! On continue...
4 >=1 ! On continue...
3 >=1 ! On continue...
2 >=1 ! On continue...
1 >=1 ! On continue...

1.6.4 Les boucles for

Une boucle for s'écrit de la manière suivante :
for i in [2,3,7]:    # De manière générale : for variable in liste:
    print (i)
résultat:
2
3
7
Comme il est fastidieux d'écrire un grand nombre de valeurs, on peut utiliser la fonction range() qui renvoie une liste d'entiers consécutifs:
for i in range(2,5):
    print i
résultat:
2
3
4
Remarque: range(5) est équivalent de range(0,5)
Exercice 1.6.2. « La Suite de Syracuse ». Soit un nombre entier x 0 1 . Par récurrence, on définit x t+1 à partir de x t par: x t+1 ={ 1 2 x t  si  x t  est pair 1 2 ( 3 x t +1 )  si  x t  est impair
Observer que la valeur de départ x 0 =1 donne la suite infinie périodique 1,2,1,2,1,
La valeur de départ x 0 =7 donne la suite 7,11,17,26,13,20,10,5,8,4,2,1,2,1,2,1,
La fameuse conjecture de Syracuse non démontrée à ce jour est que partant de toute valeur x 0 la suite arrive forcément au bout d'une date T (qui dépend de x 0 ) à la suite périodique 1,2,1,2,
image: 19_home_faure_enseignement_Systemes_dynamiques_M1_Images_chap_intro_syracuse_10.png
Ecrire un programme qui demande x 0 à l'utilisateur et affiche les valeurs de la suite ( x t ) t jusqu'à ce que la valeur x t =1 soit atteinte. Exemple de résultat:
x0=?7
11 , 17 , 26 , 13 , 20 , 10 , 5 , 8 , 4 , 2 , 1 ,
Solution:
x0 = int(input('x0=?'))

while(x0 != 1):
    if(x0%2 == 0):
        x0=x0/2
    else:
        x0=(3*x0+1)/2
    print (x0,',',end=' ')

1.7 Mesurer le temps

Ecrire:
from numpy import *
import time
 
t0 = time.time() # mesure de la date actuelle
 
for p in range(100000): # un calcul quelconque
  y = cos(p)
 
t1 = time.time() # mesure de la date
 
print ('Duree: t1-t0 = %6.3f s' %(t1-t0)) # affiche (t1-t0) avec une précision de 3 chiffres après la virgule, voir section 1.2
résultat:
Duree: t1-t0 = 0.116 s

1.8 Fonctions

Nous avons déjà utilisé les fonctions print() et type(). De manière générale, une fonction appelée f et prenant des paramètres x,y s'utilise de la manière suivante : f(x,y) ou (si la fonction renvoie une valeur y) : y=f(x,y)
Pour définir une nouvelle fonction, il faut utiliser le mot-clef “def”.
Ecrire
def MaFonction(nom,age): #definition de la fonction
    print nom,'a',age,'ans'        #noter les espaces en début de ligne
 
MaFonction('Obélix',5)   #utilisation de la fonction
résultat:
Obélix a 5 ans 
nom et age sont les paramètres de la fonction. On peut appeler la même fonction ci-dessus en utilisant une chaine de caractères pour age. Ecrire
MaFonction('Obélix','cinq')   #utilisation de la fonction
résultat:
Obélix a cinq ans 
Il faut que les types des paramètres fournis à la fonction soient conforme à l'utilisation qui va être faite. Ecrire
def Addition(a,b):
    return a+b          # Renvoie une valeur avec “return”
print Addition(1,2)        #deux entiers
print Addition(1,3.0)      #un entier plus un flottant= un flottant
print Addition('tata','toto')      # deux chaînes de caractères
résultat:
3
4.0
tatatoto 
Ecrire
print Addition(1,'toto')   # entier+chaîne de caractères = erreur !
résultat:
Traceback (most recent call last):
  File '<stdin>', line 1, in ?
  File '<stdin>', line 2, in Addition
TypeError: unsupported operand type(s) for +: 'int' and 'str'
Remarque 1.8.1. Comme pour les boucles c'est l'indentation (le nombre d'espaces en début de ligne) qui définit quelles sont les instructions qui font partie de la fonction. Cela remplace l'utilisation d'accolades {} dans d'autres langages.

1.8.1 Portée des objets : variables locales et globales

En python, les variables sont locales (par défaut), ce qui signifie qu'elles n'existent qu'au niveau où elles ont été définies.
Ecrire
x=12.3              # Première variable “x”
 
def fonction1(): # déclaration d'une fonction
    x=1              # Deuxième variable “x”
    print 'Etape 1, x=',x
    
def fonction2(): # déclaration d'une fonction
    x='Toto'         # Troisième variable “x”
    print 'Etape 2, x=',x
 
print 'Etape 0, x=',x
fonction1()
fonction2()
résultat:
Etape 0, x= 12.3
Etape 1, x= 1
Etape 2, x= Toto
Dans l'exemple ci-dessus, les 3 variables ont le même nom mais sont complètement indépendantes. On dit qu'elles sont locales. Il est à noter qu'avec un niveau d'indentation on a accès aux variables du niveau supérieur en lecture, mais on ne peut pas les modifier (contrairement au C/C++ où on peut les lire et les modifier). Essayer par exemple :
compte=0
 
def compteur():
    print (compte)
    compte = compte + 1
    return compte
 
compteur()
résultat:
0
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  File "<stdin>", line 2, in compteur
UnboundLocalError: local variable 'compte' referenced before assignment
Pour pouvoir utiliser (modifier) une variable à différents endroits (autrement qu'en la passant comme paramètre à une fonction), il faut utiliser le mot-clef global :
compte=0
 
def compteur():
    global compte
    compte = compte + 1
    return compte
 
print (compteur())
résultat:
1
Note : en général, il vaut mieux éviter d'utiliser des variables globales, et plutôt passer des variables par paramètres aux fonctions.
Exercice 1.8.2. (suite de l'exercice 1.6.2 sur la suite de Syracuse).
  1. Ecrire une fonction S( x ) qui prend un entier x en entrée et renvoit 1 2 x si x est pair ou 1 2 ( 3x+1 ) sinon.
  2. Pour x 0 donné, on note T( x 0 ) la plus petite valeur de t N telle que x t =1 . Faire une fonction T( x ) qui prend un entier x en entrée et renvoit T( x ) .
  3. Afficher les valeurs de T pour les valeurs de x 0 { 1,2, 100 } . Le résultat sera de la forme:
    x= 1 T= 1
    x= 2 T= 2
    x= 3 T= 6
    x= 4 T= 3
    x= 5 T= 5
    x= 6 T= 7
    x= 7 T= 12
    x= 8 T= 4
    x= 9 T= 14
    x= 10 T= 6
    etc
Solution:
def S(x):
    if(x%2==0):
        return x/2
    else:
        return (3*x+1)/2

def T(x):
    T0=1
    while(x!=1):
        x=S(x)
        T0 = T0+1
    return T0

for x in range(1,100):
    print ('x=', x, 'T=',T(x))

1.8.2 Récursivité

La récursivité est lorsqu'une fonction s'appelle elle même. Regarder l'exemple suivant et essayer de comprendre ce qu'elle fait.
Ecrire:
def f(n):
if(n>1):
return n*f(n-1);
else:
return 1;
print(f(5))
Résultat:
120

Chapitre 2 Objets plus élaborés

2.1 Vecteurs, matrices et algèbre linéaire

2.1.1 Vecteurs et représentation graphique

Voici comment utiliser arange() pour créer une liste de nombres. Comme pour 'list', les indices commencent à 0. Ecrire:
from pylab import *
from numpy import *
X=arange(0, 6, 0.5) # crée une liste de nombres X_{i} de 0 à 6 (exclu) par pas de 0.2 (il y a donc 6/0.2 = 30 éléments)
print ('X = ', X)
Y=sin(X) # crée une liste où chaque élément est Y_{i}=\sin\left(X_{i}\right)
print ('Y = ', Y)
 
plot(Y) # dessin des valeurs de \left(Y_{i}\right)_{i} en fonction de i
show() # montre le dessin
 
plot(X,Y) # dessin des valeurs de Yi en fonction des Xi
xlabel('x') # rajoute titre aux axes
ylabel('y=sin(x)')
show() # montre le dessin
résultat:
X = [ 0. 0.5 1. 1.5 2. 2.5 3. 3.5 4. 4.5 5. 5.5]
Y = [ 0. 0.47942554 0.84147098 0.99749499 0.90929743 0.59847214
0.14112001 -0.35078323 -0.7568025 -0.97753012 -0.95892427 -0.70554033]
image: 20_home_faure_enseignement_informatique_python_images_sinus.png
image: 17_home_faure_enseignement_informatique_python_images_sinus_xy.png
Pour connaitre le type des tableaux et listes de l'exemple précédent, écrire ensuite:
print (type(X))
résultat:
<class 'numpy.ndarray'>
Remarque 2.1.1. Ce type 'numpy.ndarray' ne peut contenir que des nombres, et pas d'autres types comme le permet 'list'.
Il y a une autre possibilité pour créer une liste de nombres avec linspace(). écrire:
from pylab import *
from numpy import *
X=linspace(0,1,10) # 10 nombres X i équidistribués de 0 à 1 compris
Y=sin(X)
plot(Y, marker='o', linestyle='solid') # options pour le dessin
show()
résultat:
image: 21_home_faure_enseignement_informatique_python_images_courbe.png
Ecrire
from numpy import *
zeros(10) # liste de 10 nombres tous égaux à zéro
Résultat:
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
Documentation: liste des markers, plot
Conversion list->array avec la fonction asarray:
Il peut être utile de convertir une liste de nombres de type 'list' en une liste de type 'numpy.ndarray' afin de la dessiner par exemple.
Ecrire:
from pylab import *
from numpy import *
Y = [0,1,4,9,16]
print (type(Y))
Y = asarray(Y)
print(type(Y))
plot(Y)
Résultat:
<class 'list'>
<class 'numpy.ndarray'>
image: 22_home_faure_enseignement_informatique_python_images_courbe_squares.png
Copie de tableaux:
On a vu que pour copier une liste L1 dans un nouvel objet L2 il faut écrire L2=list(L1) . De la même façon pour copier un tableau T1 dans un nouveau tableau T2 il faut écrire: T2 = np.copy(T1). Ecrire:
import numpy as np
T1 =arange(0, 2, 0.5) # crée une liste de nombres X_{i} de 0 à 2 (exclu) par pas de 0.5
T2 = copy(T1) # copie dans objet distinct
T3 = T1 # attention, ce n'est pas une vraie copie: T3 sera seulement un autre nom pour T1
T2[0] = 10
T3[0] = 20
print ('T1 = ', T1)
print ('T2 = ', T2)
print ('T3 = ', T3)
Résultat:
T1 = [ 20. 0.5 1. 1.5]
T2 = [ 10. 0.5 1. 1.5]
T3 = [ 20. 0.5 1. 1.5]

2.1.2 Matrices et représentation graphique

2.1.2.1 Premier exemple

Remarque 2.1.2. Attention les indices des tableaux commencent à 0.
 
from pylab import *
from numpy import *
 
A=zeros([3,3]) # matrice 3 × 3 remplit de zeros
A[1,1]=2 # on modifie un élément
A[2,1]=3
 
print (A)
 
pcolormesh(A) # dessin des contours des valeurs
colorbar() # rajoute une barre des valeurs
show()
résultat:
 
[[0. 0. 0.]
[0. 2. 0.]
[0. 3. 0.]]
image: 23_home_faure_enseignement_informatique_python_images_im1.png

2.1.2.2 Autres représentations graphiques

Regarder ces exemples.
Surface avec ombres
Illuminating surface plots

2.1.2.3 Création de matrices

Matrice avec valeurs aléatoires:
Ecrire
from pylab import *
from numpy import *
random.uniform(-1, 1, size=(3,3)) # matrice 3*3 avec éléments aléatoires dans [ -1,1 ]
Résultat:
array([[ 0.40517391, -0.47365213, -0.5208182 ],
[-0.23394281, 0.0347132 , -0.92326085],
[-0.94118242, -0.58734857, -0.19587783]])
Matrice identité:
 
from pylab import *
from numpy import *
eye(3)
Résultat
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])

2.1.2.4 Matrice comme produit tensoriel de vecteurs

L'exemple suivant montre comment à partir de deux vecteurs ( x i ) i et ( y j ) j on peut construire une matrice ( z i,j ) i,j . On termine l'exemple par une représentation graphique.
 
from pylab import *
from numpy import *
 
xmin,xmax,ymin,ymax = -3,3, -3,3 # bornes du domaine
M = 30 # nombre de points en x et y
 
x = linspace(xmin,xmax,M) # vecteur avec M nombres entre xmin et xmax
y = linspace(ymin,ymax,M)[:,newaxis] # vecteur dans la 2eme dimension
z = exp( -x*x-y*y) # matrice
 
levels = linspace(0,1,10) # defini les lignes de niveaux
contourf(z,levels,extent=(xmin,xmax,ymin,ymax))
colorbar()
xlabel('$x$')
ylabel('$y$')
title('fonction $f(x,y)= \exp (-x^2 -y^2)$. Echantillons: $M='+ ('%4i'%M)+'$')
Résultat:
image: 24_home_faure_enseignement_informatique_python_images_gaussienne.png

2.1.3 Opérations d'algèbre linéaire

Voir scipy, linear algebra et scipy, linear algebra.

2.1.3.1 opérations élémentaires

Ecrire
import numpy as np
from scipy import linalg
A = np.array([[2,1],[1,1]]) #matrice
V = np.array([[1], [1]]) # vecteur colonne
print ('A=\n', A , '\nV=\n', V)
print ('A*V=\n', np.dot(A,V)) # produit matrice * vecteur
print ('A*A=\n', np.dot(A,A)) # produit matrice * matrice
d = linalg.det(A) #calcule le determinant
print ('det(A) = ', d)
Résultat:
A=
[[2 1]
[1 1]]
V=
[[1]
[1]]
A*V=
[[3]
[2]]
A*A=
[[5 3]
[3 2]]
det(A) = 1.0

2.1.3.2 Valeurs propres et vecteurs propres

On rappelle que si M est une matrice diagonalisable n × n alors sa diagonalisation donne M=PD P -1 D=( d 0 0 0 d n-1 ) est la matrice diagonale des valeurs propres et P=( ( U 0 ) 0 ( U n-1 ) 0 ( U 0 ) n-1 ( U n-1 ) n-1 ) est une matrice contenant en colonnes les composantes des vecteurs propres (à droite) U 0 , U n-1 , qui vérifient M U j = d j U j ,j{ 0, n-1 }
Ecrire:
from pylab import *
from numpy import *
M = array(([1,3,3],[1,4,3],[1,3,4])) # matrice
D, P = eig(M) # -> valeurs propres D et matrice des vecteurs propres P
R = dot( dot(P , diag(D)) , inv(P)) # R = P * D * P^(-1)
print ('M=\n' ,M, '\n P*D*P^(-1) =\n' , R)
Résultat:
M=
[[1 3 3]
[1 4 3]
[1 3 4]]
 
P*D*P^(-1) =
[[1. 3. 3.]
[1. 4. 3.]
[1. 3. 4.]]
Ecrire:
import pylab as py
import numpy as np
 
N = 100
A= 0.1
M = np.random.uniform(0, A, size=(N,N)) # matrice N*N avec éléments aléatoires dans [0,A]
EV = py.eigvals(M) # calcul des valeurs propres de M
 
# dessin des valeurs propres ------------------
plot(real(EV), imag(EV), linestyle = 'none', marker = 'o', c = 'lime',
markersize = 10)
a = 0.5*sqrt(N)*A
xmin,xmax,ymin,ymax = -a,a,-a,a
plt.axis([xmin,xmax,ymin,ymax]) # selectionne la vue
xlabel('$\Re(z)$')
ylabel('$\Im(z)$')
title('Eigenvalues of a random $' + ('%d'%N) +'*'+ ('%d'%N) +'$ matrix with iid elements in $[0,' + ('%4.2f'%A) +']$')
Résultat:
image: 25_home_faure_enseignement_informatique_python_images_spectre.png

2.1.3.3 Matrices creuses (sparse)

Dans de nombreux problèmes on doit utiliser des matrices de taille très grande, ( 1 0 4 ou plus) mais avec seulement un petit nombre d'éléments non nuls. Dans ce cas l'ordinateur stocke seulement la liste des éléments non nuls. Ce mode de stockage s'appelle « matrice creuse ».
Matrices sparses avec scipy en python.

2.2 Objets graphiques avec Matplotlib

Lire ce site ou ici. Site de matplotlib. Voir Samples with matplotlib. Exemples. Cliquer sur les images pour voir le code. Exemples avec mathplotlib
Programmes avec graphisme:
Afin de pouvoir interagir avec les fenetres graphiques, dans spyder il est conseillé, dans le menu Outils/Préférences/ConsoleIPython/Graphiques/Sortie : de choisir le mode automatique.
(et relancer la fenetre IPython)

2.2.1 Courbe 1Dim simple

Exemple ici.
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
 
t = np.arange(0.0, 2.0, 0.01) # liste de nombre de 0 a 2 par pas de 0.01
s = 1 + np.sin(2 * np.pi * t)
 
plt.plot(t, s) # points relies (bleus)
plt.plot(t+0.1, s, 'r.') # red (r) et points isoles (.)
plt.plot(t+0.2, s, 'gx') # green (g) et croix (x)
plt.xlabel('t(s):temps')
plt.ylabel('V(mV): Tension')
plt.title('Titre')
plt.grid()
plt.savefig('V-t.png') # sauvegarde fichier image au format png
plt.show()
Résultat:
image: 26_home_faure_enseignement_informatique_python_images_V-t.png
Documentations spécifiques:
sur figure(),

2.2.2 Des sous figures, titres et labels sur les axes

Référence: Titres et labels sur un graphe.
Ecrire:
import matplotlib.pyplot as plt
 
#-- preparation des figures
fig = plt.figure(1, figsize=(10, 5)) # taille X,Y en pouces (inchs)
plt.subplots_adjust(left=0.1, wspace=0.4, top= 0.8)
plt.suptitle('Titre general', fontsize=16)
 
#--- figure 1
xmin,xmax,ymin,ymax = -2,2,-2,2 # fenetre
f1 = plt.subplot(121) # Il y a 1*2 sous figures. f1 est la figure 1
plt.title('Titre 1')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
f1.plot([0,-1.5], [0,1.5]) # segment (0,0)->(-1.5,1.5)
f1.axis([xmin,xmax,ymin,ymax]) # selectionne la vue
 
#--figure 2
amin, amax, bmin, bmax = -3,3, -2, 2
f2 = plt.subplot(122) # Il y a 1*2 sous figures. f2 est la figure 2
plt.title('Titre 2')
plt.xlabel('a')
plt.ylabel('b')
plt.grid(True)
f2.plot([-2,2], [-1,1]) # segment (-2,-1)->(2,1)
f2.axis([amin,amax,bmin,bmax]) # selectionne la vue
 
plt.tight_layout() # pour que les differentes sous figures ne se superposent pas
Résultat:
image: 27_home_faure_enseignement_informatique_python_images_sous_figures.png

2.2.3 Des objets graphiques en dimension 2

Référence, exemples
Ecrire:
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.path as mpath
import matplotlib.lines as mlines
import matplotlib.patches as mpatches
from matplotlib.collections import PatchCollection
from math import *
 
fig, ax = plt.subplots()
 
#texte
x,y = 0.5, 0.3
plt.text(x, y, 'Du texte', ha="center", family='sans-serif', size=14)
 
# Un cercle
x,y = 0.5,0.5
r = 0.1
c = plt.Circle((x,y), r, facecolor= 'green', edgecolor= 'red')
ax.add_patch(c)
 
#Une ellipse
x,y = 0.2,0.3
rx,ry = 0.2, 0.1
angle = 30 # en degres
e = mpatches.Ellipse((x,y), rx, ry, angle, facecolor= 'yellow', edgecolor= 'red')
ax.add_patch(e)
 
# un rectangle
x,y = 0.,0. # coin inferieur gauche
dx,dy=0.3, 0.1 # taille
rec = plt.Rectangle((x,y), dx,dy, facecolor= 'blue', edgecolor= 'black')
ax.add_patch(rec)
 
#Une flèche
x,y = -0.2, 0. # départ
dx,dy=0.2, 0.3 # deplacement
width = 0.1
a = mpatches.Arrow(x, y, dx, dy, width, facecolor= 'pink', edgecolor= 'red')
ax.add_patch(a)
 
plt.axis('equal') # pour avoir meme echelle en x,y
#plt.axis('off') # n'affiche pas le cadre ni les axes
ax.autoscale_view() # pour que la figure s'ajuste aux objets
plt.show()
Résultat (selon les options):
image: 28_home_faure_enseignement_informatique_python_images_objets_2D.png image: 29_home_faure_enseignement_informatique_python_images_objets_2D_2.png

2.2.4 Dessin d'un champ de vecteur en dimension 2 avec streamplot() ou quiver()

Reference.
Ecrire:
import numpy as np
import matplotlib.pyplot as plt
 
#fonction qui definit le champ de vecteur F
def F(x, y):
  return x, -y
 
# reseau de points en x,y
xmin, xmax, dx = -2,2, 0.2
ymin, ymax, dy = -2,2, 0.2
X,Y = np.meshgrid( np.arange(xmin, xmax, dx), np.arange(ymin, ymax, dy) )
Fx, Fy = F(X,Y) # calcule le champ de vecteur sur le reseau
 
#----------------
plt.streamplot(X,Y, Fx,Fy)
plt.title('avec streamplot')
 
#----------------
plt.figure() # nouvelle figure
norme = np.sqrt(Fx**2+Fy**2) # norme du champ F
strm = plt.streamplot(X, Y, Fx, Fy, color = norme, linewidth = norme, cmap=plt.cm.inferno, density=[1, 2], arrowstyle='->', arrowsize=1.5)
plt.colorbar(strm.lines) # montre l'echelle
plt.title("avec streamplot et plus d'options")
 
#----------------
plt.figure() # nouvelle figure
plt.quiver(X,Y,Fx,Fy)
plt.title('avec quiver')
Résultat:
image: 30_home_faure_enseignement_informatique_python_images_champ_vecteur_1.png image: 31_home_faure_enseignement_informatique_python_images_champ_vecteur_2.png image: 32_home_faure_enseignement_informatique_python_images_champ_vecteur_3.png

2.2.4.1 Dessin des lignes de champ passant par des points spécifiés avec streamplot()

Ecrire:
import numpy as np
import matplotlib.pyplot as plt
 
#definit le champ de vecteur F
def F(x, y):
return x, -y
 
# reseau de points en x,y
nx, ny = 64, 64 # nombre de points
x = np.linspace(-2, 2, nx) # 1dim array
y = np.linspace(-2, 2, ny)
X, Y = np.meshgrid(x, y) # 2dim array
Fx, Fy = F(X,Y) # calcule le champ de vecteur sur le reseau
 
# Liste de points (x,y)
points = np.array([[1, 1, -1, -1, 1], [0, 1, 1, -1, -1]])
plt.streamplot(x, y, Fx, Fy, start_points = points.T)
plt.plot(points[0], points[1], linestyle ='none', marker='o') # dessin des points
plt.title('Lignes de champ passant par des points spécifiques')
Résultat:
image: 33_home_faure_enseignement_informatique_python_images_champ_vecteur_4.png

2.2.5 Créer une animation

Reference.
Dans l'exemple suivant on considère N nombres complexes z j ( x ) = exp ( i2 π ( j+x ) /N ) avec les indices j=0,1, N-1 et un déphasage fixé x[ 0,1 ] . Ces N nombres sont équidistribués sur le cercle unité | z j ( x ) | =1 et on remarque que z j ( 1 ) = z j+1 ( 0 ) . Ainsi en faisant varier x=01 on verra les N points tourner. Remarquer que X j ( x ) = Re ( z j ( x ) ) = cos ( 2 π ( j+x ) /N ) et Y j ( x ) = Im ( z j ( x ) ) = sin ( 2 π ( j+x ) /N ) .
Ecrire:
import matplotlib.pyplot as plt
import numpy as np
import subprocess
from math import *
 
N, P = 5, 10 # nombre de points et nombre d'images,
xmin,xmax,ymin,ymax = -2,2,-2,2
 
for p in range(P): # boucle p = 0,1..,P-1
  x = p/P
  j = np.arange(N) # liste 0->(N-1)
  X = np.cos(2*pi*(j+x) /N) # liste des X_j(x)
  Y = np.sin(2*pi*(j+x) /N) # liste des Y_j(x)
  plt.cla() # efface le graphisme precedent
  plt.plot(X,Y, linestyle='none', marker='o')
  plt.axis([xmin,xmax,ymin,ymax]) # selectionne la vue
  plt.pause(0.01) # montre la figure et attend 0.01 sec.
 
plt.show() #laisse la figure à la fin et attend qu'on la referme
Pour voir le résultat, ouvrir un terminal (fenêtre de commandes) et lancer le programme par: python3 programme.py
Résultat:
image: 34_home_faure_enseignement_informatique_python_images_animation.gif

2.2.6 Créer une animation et crée un fichier gif animé avec convert ou gifsicle

Voici le même programme que ci-dessus, mais avec deux instructions supplémentaires (savefig() qui sauve chaque image dans la boucle et convert... qui fabrique la vidéo finale), afin de fabriquer à la fin un fichier gif animé. Il faut avoir installé le logiciel imagemagick.
import matplotlib.pyplot as plt
import numpy as np
import subprocess
from math import *
 
N, P = 5, 10 # nombre de points et nombre d'images,
xmin,xmax,ymin,ymax = -2,2,-2,2
 
for p in range(P): # boucle p = 0,1..,P-1
  x = p/P
  j = np.arange(N) # liste 0->(N-1)
  X = np.cos(2*pi*(j+x) /N) # liste des X_j(x)
  Y = np.sin(2*pi*(j+x) /N) # liste des Y_j(x)
  plt.cla() # efface le dessin precedent
  plt.plot(X,Y, linestyle='none', marker='o')
  plt.axis([xmin,xmax,ymin,ymax]) # selectionne la vue
  plt.savefig('image_' + str(p).zfill(4) + '.png', dpi=50) # sauve fichier image (pour ensuite gif anime). dpi: precision de l'image
  plt.pause(0.0O1) # montre la figure et attend 0.001 sec.
 
#--- cree un fichier gif anime. delay en 1/100s.-----------------------------
subprocess.getoutput('convert -delay 10 -loop 0 image_*.png animation.gif') #fabrique fichier gif final
subprocess.getoutput('rm image_*.png') # efface les fichiers images temporaires
print ('fichiers animation.gif cree. Veuillez le lire avec firefox par exemple.')
Pour voir le résultat, ouvrir un terminal (fenêtre de commandes) et lancer le programme par: python3 programme.py
Résultat:
fichier animation.gif cree. Veuillez le lire avec firefox par exemple.
image: 34_home_faure_enseignement_informatique_python_images_animation.gif
Pour visualiser le fichier produit, il faut l'ouvrir avec un navigateur (firefox par exemple). Ce fichier se trouve dans votre répertoire de travail.

2.2.6.1 Avec Gifsicle:

On peut aussi utiliser gifsicle pour fabriquer le fichier gif animé final. L'utilisation de gifsicle est conseillée, car les fichiers sont plus compressés.
Il faut remplacer la ligne subprocess.getoutput('convert -delay 10 -loop 0 image_*.png animation.gif')
par:
subprocess.getoutput('convert image_*.png GIF:- | gifsicle --delay=10 --loop --optimize=2 --colors=256 --multifile - > animation.gif')

2.2.7 Capturer les évènements de la souris (clic, position) et du clavier

Documentation, voir cet exemple et ces autres exemples.
Avec le programme on crée un point rouge en cliquant sur la fenetre. On déplace le point rouge avec les 4 flèches du clavier.
Ecrire:
import numpy as np
import matplotlib.pyplot as plt

x, y = 0,0 # point


#---- si evenement souris (clik)
def on_click(event):
    global x,y
    x, y = event.xdata, event.ydata # coord souris en x,y
    print('click en x=',x,' y=',y)
    plt.plot(x,y, 'r.')
    plt.pause(0.001) # montre le dessin sans bloquer            
    

#---- si evenement clavier
def on_key(event):
    print('key = ', event.key)
    global x,y
   # print('you pressed', event.key, event.xdata, event.ydata)
    if event.key == 'right':
            x=x+1
    elif event.key == 'left':
            x=x-1
    elif event.key == 'up':
            y=y+1
    elif event.key == 'down':
            y=y-1

    plt.plot(x,y, 'r.')
    plt.pause(0.001) # montre le dessin sans bloquer
    
#----------------------------    
plt.axis([-10,10,-10,10]) # selectionne la vue
#plt.axis('equal') # pour avoir meme echelle en x,y
plt.connect('button_press_event', on_click) # associe la fonction aux evenements  souris
plt.connect('key_press_event', on_key)  # associe la fonction aux evenements clavier
plt.show() # bloque le programme ici pour gérer les evenements

2.3 (Supplément) Dictionnaires (dict)

Un dictionnaire est également un type de liste d'objet, mais au lieu de repérer les objets par un indice commençant à 0, les différents objets sont accessible à l'aide d'une clef, qui peut être un objet quelconque.
Ecrire:
fra_deu={'un':'ein','trois':'drei','deux':'zwei'}    # Création
fra_deu['trois']            # Accès à un élément du dictionnaire
résultat:
'drei'
Ecrire:
fra_deu['quatre']              # Clef n'existant pas !
résultat:
Traceback (most recent call last):
 File '<stdin>', line 1, in ?
KeyError: 'quatre'
Ecrire:
fra_deu['quatre']='vier'    # Ajout de l'entrée manquante
fra_deu['quatre']
résultat:
'vier'
Vous pouvez obtenir les différentes fonctions disponibles pour un dictionnaire avec la commande help(dict).

Chapitre 3 Lecture et écriture de fichiers

On crée un objet fichier Python avec la fonction open() :
MonFichier = open('NomDuFichier','r')
Le mode peut être 'r' (read, en lecture seule), 'w' (en écriture seule), 'a' pour rajouter à la fin du fichier, D'autres modes (mixte lexture/écriture, etc) existent. La lecture des fichiers texte (i.e. non binaires) se fait avec les fonctions read(), readline() et readlines():
MonFichier.read()     # Lit la totalité du fichier sous forme d'un chaîne de caractères
MonFichier.readline()    # Lit la ligne suivante du fichier, sous forme d'une chaîne de caractères
MonFichier.readlines()    # Lit toutes les lignes du fichier sous forme d'une liste de chaînes de caractères
Monfichier.write('toto') # écrit
Pour lire une liste de nombres contenu dans une ligne, il suffit de lire cette ligne avec la fonction readline(), puis de séparer les nombres avec la fonction split(). Par exemple :
f=open('tmp0','r')
ligne=f.readline()       # lecture d'une ligne
print (ligne)      # Ecriture de la ligne sous forme d'une chaîne de caractères
print ligne.split()    # On sépare les éléments -> liste
print float(ligne.split()[1]) # Conversion en nombre réel 0.2
f.close()      # referme le fichier
résultat:
Toto 1.8 9.2
['Toto', '0.2', '9.2']
help(file) permet de lister l'aide sur l'utilisation des fichiers.
Ecriture/Lecture de données formatées:
On utilise le module pickle. Cela permet d'écrire des données et lire des données en conservant leur type (nombre, liste, chaine de caractère etc).
Pour écrire:
import pickle
x=2
f=open('fichier','w')
pickle.dump(x,f) #écrit x
f.close()
Pour relire:
f=open('fichier','r')
x=pickle.load(f) #lit x
f.close()
print x
résultat:
2

Chapitre 4 Résoudre une équation différentielle ordinaire (EDO) avec Scipy

Généralités sur les équations différentielles ordinaires (EDO)

Une équation différentielle ordinaire (EDO) est une équation de la forme dy dt =f( y,t ) où l'inconnue est la fonction y à n composantes: y:t R y( t ) R n , alors que la fonction f:( y,t ) R n × R R n et la condition initiale y( 0 ) R n sont supposées connues.

4.1.1 Transformer une équation différentielle d'ordre d en EDO

Prenons un exemple. L'équation du deuxième ordre en temps: d 2 y d t 2 = dy dt + sin ( y ) peut se ré-écrire comme une équation du premier ordre dY dt =F( Y,t ) à deux variables en posant Y=( y 1 , y 2 ) et y 1 =y , y 2 = dy dt et F( Y,t ) =( y 2 , y 2 + sin ( y 1 ) ) .
Plus généralement une équation différentielle d'ordre d à 1 variable y( t ) : d d y d t d =f( y( t ) , dy dt , d d-1 y d t d-1 ,t ) peut se réécrire comme une équation d'ordre 1 avec d variables Y=( y 1 ( t ) , , y d ( t ) ) . Il suffit de poser y 1 ( t ) =y( t ) , y 2 ( t ) = dy dt , , y d ( t ) = d d-1 y d t d-1 , c'est à dire d y 1 dt ( t ) = y 2 ( t ) , , d y d-1 ( t ) dt = y d ( t ) , d y d ( t ) dt =f( y 1 , y 2 , , y d-1 ,t ) . Ainsi l'équation se ré-écrit: dY dt =F( Y,t ) avec F( Y,t ) =( y 2 , y 3 , y d ,f( y 1 , y d-1 ,t ) ) .

Résoudre numériquement une EDO avec le langage python

Référence de odeint ou reference de ode.
On va utiliser une fonction appelée odeint(). Le principe est de donner les conditions initiales y( 0 ) , d'écrire la fonction f( y,t ) , et de donner un tableau de valeurs du temps ( t j ) j=1N . En retour la fonction odeint() renvoit un tableau de valeurs ( y( t j ) ) j=1N tel que y( t ) est la solution de l'ODE dy dt =f( y,t ) avec les conditions initiales spécifiées.
On peut passer des paramètres extérieurs à la fonction.

4.2.1 Exemple 1

On considère l'équation différentielle ordinaire à n=1 composante suivante où l'inconnue est la fonction y:t R y( t ) R : dy dt =y( t ) sin ( t )
Remarque 4.2.1. L'équation s'écrit dy dt =f( y,t ) avec la fonction est f( y,t ) =y sin ( t ) . La solution est y( t ) =y( 0 ) e 1- cos ( t ) .
Ecrire:
 
%matplotlib inline
from scipy.integrate import odeint
from pylab import *
from numpy import *
 
y_in = 0.1 # condition initiale
tmax= 50 # intervalle de temps
N = 501 # nombre de points en temps
 
def f(y,t):
  return y*sin(t)
 
t = linspace(0, tmax, N) # tableau de t = 0->tmax avec N valeurs
y = odeint(f, y_in,t)
 
plot(t,y)
xlabel('$t$')
ylabel('$y(t)$')
title("Solution de $y'(t)=y sin(t)$")
Résultat:
image: 35_home_faure_enseignement_informatique_python_images_dydt.png

4.2.2 Exemple 2

(4.2.2).On considère l'équation différentielle ordinaire à n=3 composantes suivante appelées équations de Lorenz. L'inconnue est la fonction y:t R y( t ) =( y 0 ( t ) , y 1 ( t ) , y 2 ( t ) ) R 3 :
{ d y 0 dt =- σ ( y 0 ( t ) - y 1 ( t ) ) d y 1 dt = ρ y 0 ( t ) - y 1 ( t ) - y 0 ( t ) y 2 ( t ) d y 2 dt =- β y 2 ( t ) + y 0 ( t ) y 1 ( t ) On connait les conditions initiales y( 0 ) et les paramètres σ =10, β =8/3, ρ =28 .
Ecrire:
#%matplotlib inline
from numpy import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
 
sigma, beta, rho = 10, 8/3., 28 # parametres
y_in = [0, 1, 1.05] # Conditions initiales
tmax= 100 # intervalle de temps
N = 10000 # nombre de points en temps
 
def f(y,t,sigma, beta, rho):
  """The Lorenz equations."""
  f0 = -sigma*(y[0] - y[1])
  f1 = rho*y[0] - y[1] - y[0]*y[2]
  f2 = -beta*y[2] + y[0]*y[1]
  return f0,f1,f2
 
t = linspace(0, tmax, N) # tableau de t = 0->tmax avec N valeurs
y = odeint(f, y_in, t, args=(sigma, beta, rho))
 
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(y.T[0], y.T[1], y.T[2])
#plt.show()
 
# rotate the axes and update
for angle in range(0, 360):
  ax.view_init(30, angle)
  plt.draw()
  plt.pause(.001) # montre image
Résultat: (on a fabriqué ce fichier gif animé d'après la Section 2.2.6 avec gifsicle)
image: 36_home_faure_enseignement_informatique_python_images_lorenz_anime.gif

Chapitre 5 Micro-projets

Dans ce chapitre, nous vous proposons de réaliser des micro-projets. Ce sont des programmes rapides à réaliser (quelques lignes de code), afin de vous montrer l'aspect ludique et créatif de la programmation, et aussi pour résumer différentes notions vues dans les chapitres précédents.
Nous demandons de produire un compte rendu sous forme pdf et html, que l'on rédigera avec Lyx (qui permet de créer des documents scientifiques) et exportée en html (avec LyxHtml) et en pdf: pour apprendre Lyx, voici un document en pdf qui sert d'exemple et d'exercice. Voici le même document en html, qui vous servira éventuellement pour copier/coller. Dans ce compte rendu qui sert d'exercice, il devra y avoir au moins une équation en ligne, une équation dans le texte, une figure, et une table des matières.
La Section suivante donne quelques recommandations générales et prélables pour effectuer un projet de programmation.

5.1 Introduction aux micro-projets

5.1.1 Que faire avant de programmer?

Avant de se retrouver devant un ordinateur pour écrire un programme, il faut avoir établi clairement ce qu'on désire lui faire faire. Voici l'ensemble des étapes à suivre, indispensables, à faire au préalable devant une feuille de papier.
  1. Définir clairement l'objectif du programme (surtout s'il y a plusieurs personnes à y collaborer).
  2. Définir l'ensemble des tâches qui doivent être accomplies dans un ordre logique pour atteindre cet objectif. Pour un projet scientifique il faut que toutes les formules mathématiques soient bien écrites.
  3. Tracer sur un papier un organigramme illustrant le fonctionnement de votre programme. Vérifiez que l'ordre d'exécution des diverses tâches va effectivement faire ce que vous attendez.
    En fait, votre organigramme se traduira directement dans la partie principale de votre programme (fonction int main()). A chaque tâche correspondra l'appel d'une sous partie (une fonction).
  4. Déterminez la méthode (algorithme) permettant de remplir chaque tâche. Il est possible qu'une sous-tâche soit commune à plusieurs fonctions. Dans ce cas, faites-en une nouvelle fonction. Cependant, pour éviter une multiplication de ces fonctions, ne le faites que pour celles faisant plusieurs lignes de programme. L'art de la programmation consiste à trouver un compromis entre la simplicité, la longueur, et la rapidité d'éxécution de votre programme.
  5. Si le programme effectue des taches complexes et “imprévisibles”, imaginer des situations simples où l'on connait la solution et qui permettront de tester le programme.

5.1.2 Quelques notions de qualité en informatique

La qualité d'un programme ne se traduit pas simplement en termes d'efficacité mais également en termes de sécurité (fiabilité) et d'exportabilité. Autrement dit, vous devez écrire votre programme de telle sorte (1) que n'importe qui puisse le lire et comprendre ce qu'il fait et (2) qu'il puisse être aisément testé. Voici quelques indications:

5.2 La fractale du dragon

Difficulté: moyenne.
Si on prend un ruban de papier, et que l'on replie le coté gauche sur le côté droit, puis ainsi N fois en rabattant chaque fois la gauche sur la droite, que l'on déplie le ruban et impose à chaque pli de faire un angle droit (90°) on obtient les formes suivantes assez surprenantes, ici pour N=1,2,3,4,5 :
image: 37_home_faure_enseignement_informatique_python_micro-projets_dragon_pli.JPG image: 38_home_faure_enseignement_informatique_python_micro-projets_dragon_im_2.JPG image: 39_home_faure_enseignement_informatique_python_micro-projets_dragon_im_3.JPG
image: 40_home_faure_enseignement_informatique_python_micro-projets_dragon_im_4.JPG image: 41_home_faure_enseignement_informatique_python_micro-projets_dragon_im_5.JPG
L'objectif de ce micro-projet est de dessiner ces figures pour N quelconque, de façon automatique et d'observer l'apparition de la fractale du dragon. Voir aussi suite de pliage.
image: 42_home_faure_enseignement_informatique_python_images_dragon_N_1.png image: 43_home_faure_enseignement_informatique_python_images_dragon_N_2.png image: 44_home_faure_enseignement_informatique_python_images_dragon_N_3.png image: 45_home_faure_enseignement_informatique_python_images_dragon_N_4.png image: 46_home_faure_enseignement_informatique_python_images_dragon_N_5.png image: 47_home_faure_enseignement_informatique_python_images_dragon_N_6.png image: 48_home_faure_enseignement_informatique_python_images_dragon_N_10.png image: 49_home_faure_enseignement_informatique_python_images_dragon_N_16.png
image: 50_home_faure_enseignement_informatique_python_images_animation_dragon.gif
Exercice 5.2.1. (Mathématiques) On peut coder la forme à N fixé, comme une suite S N de 2 N -1 valeurs +1,-1 , où +1 signifie que le pli est dans le sens direct, et -1 dans le sens indirect. D'après les images ci-dessus, S 1 ={ 1 } , S 2 ={ 1,1,-1 } , S 3 ={ 1,1,-1,1,1,-1,-1 } . Montrer que la suite S N se déduit facilement de la suite S N-1 par une formule simple.
Solution 5.2.2. On a S N = S N-1 { -1 } S N-1 ¯ - S N-1 ¯ est la suite S N-1 écrite à l'envers avec l'opposé des éléments.
Exercice 5.2.3. « Dragon, étape 1 ».(Programmation sans graphisme) Ecrire un algorithme sur papier puis un programme qui construit la suite S N . Afficher la liste pour vérifier. Aide: on peut utiliser le type 'list' et chercher (sur internet ou dans le cours) comment copier une liste, inverser une liste et concaténer des listes.
Solution 5.2.4. dragon_etape_1.py
 
Exercice 5.2.5. « Dragon, étape 2 »(Programmation avec graphisme) Ensuite, pour la partie graphique, écrire un algoritme sur papier puis un programme qui dessine la fractale du dragon. Il suffit de lire la liste S N et dessiner des segments à la suite.
dragon.py

5.3 Algorithme de Monte-Carlo pour le modèle d'Ising

Difficulté: moyenne à difficile.
image: 51_home_faure_enseignement_informatique_python_images_animation_ising_03.gif image: 52_home_faure_enseignement_informatique_python_images_animation_ising_04.gif image: 53_home_faure_enseignement_informatique_python_images_animation_ising_05.gif
Figure 5.1: Résultat de l'algorithme de Monte-Carlo pour le modèle d'aimantation de Ising, aux températures β =1/( k b T ) =0.3 (désordre), 0.4 (transition de phase) et 0.5 (ordre magnétique). On observe une « transition de phase ».
(Explication physique: voir Cours de Master 1: système dynamiques, chapitre “Dynamique probabiliste sur les graphes”).
On considère un réseau N × N périodique dont les sites sont notés X=( x,y ) et dont les variables sont f X = ± 1 et modélisent des spins (up/down). Si deux sites sont voisins (chaque site a 4 voisins) on note XY . Une configuration des spins est un champ donné: f= ( f X ) X . Son énergie ferromagnétique est: E( f ) := X Y,XY ( - f X . f Y )
Remarque 5.3.1. Il y a N 2 sites et 2 N 2 configurations possibles.
La condition de périodicité du réseau signifie par exemple que le voisin de gauche du point ( 0,j ) est ( N-1,j ) .
Les configurations d'énergie minimale sont celles donnant ( - f X . f Y ) minimal soit f X . f Y =1 pour tous X,Y . Il y a donc deux configurations d'énergie minimale qui sont: f up ={ f X =1 pour tous site X } et f down ={ f X =-1 pour tous site X }.
Les configurations d'énergie maximale sont celles donnant ( - f X . f Y ) maximal soit f X . f Y =-1 pour tous X,Y . Cela est possible si N est pair avec des spins dont les valeurs sont alternées (il y a deux configurations).

5.3.1 Algorithme de Monte-Carlo

L'algorithme d'évolution suivant est appelé algorithme de montecarlo . β 0 est un paramètre fixé qui est l'inverse de la température: β =1/( k b T ) .
  1. A l'instant n=0 , on part d'une configuration f choisie au hasard.
  2. On choisit un site X 0 au hasard, et on note f' la configuration identique à f sauf au site X 0 où le spin est opposé: f ' X 0 =- f X 0 (spin opposé).
  3. Choix de la configuration à l'instant n+1 :
    1. Si E( f' ) <E( f ) on choisit la configuration f' .
    2. Si E( f' ) E( f ) on choisit la configuration f' avec la probabilité P ff' = exp ( - β .( E( f' ) -E( f ) ) ) (et on choisit f avec la probabilité complémentaire).
  4. On revient en (2) pour poursuivre l'évolution du champ de spins.
Remarque 5.3.2. Cet algorithme correspond à une matrice stochastique réversible pour la mesure de Boltzmann: u( f ) = 1 Z exp ( - β E( f ) ) image: 54_home_faure_enseignement_informatique_python_images_ising_energie_2.png
Exercice 5.3.3.  
  1. Montrer que la condition ( E( f' ) -E( f ) ) s'écrit simplement en fonction du site X et de ses voisins.
  2. Ecrire un programme qui partant d'un champ se spins initial aléatoire f , le fait évoluer selon l'algorithme de Monte-Carlo et le représente après chaque N 2 itérations.
  3. L'aimantation de la configuration f est M( f ) := X f X . Représenter l'aimantation au cours du temps. Quand l'équilibre statistique est atteind, calculer la moyenne M et la variance V= ( M- M ) 2 . Représenter M ( β ) et V( β ) en fonction de β .
Solution 5.3.4. ising.py

5.4 Le pendule (système dynamique, ODE)

L'objectif est de dessiner (créer une animation) dans l'espace réel et dans l'espace des phases, le mouvement d'oscillation d'un pendule simple de masse m , longeur l , fixées, soumis à la force de pesanteur P=mg avec g=9.8m/ s 2 et à une force de frottement de coefficient γ .
image: 55_home_faure_enseignement_informatique_c++_cours_c++_cours_1_pendule.png
image: 56_home_faure_enseignement_informatique_python_images_animation_pendule_euler.gif

5.4.1 Equations physiques

si a( t ) est la position angulaire et b( t ) = da dt la vitesse angulaire au temps t on a da dt = F a ( a,b ) :=b (1) db dt = F b ( a,b ) :=- g l sin ( a ) - γ m b F a ( a,b ) , F b ( a,b ) sont des fonctions de a,b . L'espace des phases est la plan ( a,b ) .
La position du pendule dans le plan vertical est:
x=l sin ( a ) ,z=-l cos ( a )
Exercice 5.4.1. (Mathématiques sur le papier).
  1. Vérifier les formules (1) en appliquant la loi de Newton sur la droite verte tangente au cercle.
  2. Pour les petites oscillations, sans frottement, c'est à dire a,b1 , et γ =0 , linéariser les équations du mouvement et montrer que la trajectoire du point a( t ) ,b( t ) R 2 est une ellipse parcourue dans le sens indirect avec une période T=2 π l g . (2)
Solutions:
  1. Sur l'axe tangent au cercle, la loi de Newton s'écrit: -mg sin a- γ d( la ) dt =m d 2 ( la ) d t 2 -mg sin a- γ lb=ml db dt db dt =- g l sin a- γ m b
  2. Pour γ =0 , les équations linéarisées sont sin aa donc
da dt =b (3) db dt =- g l a On pose z( t ) =a( t ) +i l g b( t ) C alors dz dt = da dt +i l g db dt =b-i g l a =-i g l ( a+i l g b ) =-i g l z( t ) dont la solution est z( t ) =z( 0 ) e -i g l t =z( 0 ) e -i2 π t/T (mouvement circulaire indirect) avec la période T=2 π l g .

5.4.2 Méthode de résolution numérique de Euler

On a da dt = F a . Si dt1 , alors au premier ordre on peut écrire, d'après la formule de Taylor: a( t+dt ) =a( t ) +( da dt ) ( t ) dt+O( d t 2 ) =a( t ) + F a dt+O( d t 2 ) et de même pour b .
La méthode de Euler, consiste à considèrer un petit intervalle de temps dt1 et de négliger le terme O( d t 2 ) . Ainsi on écrit: a( t+dt ) =a( t ) + F a dt (4) b( t+dt ) =b( t ) + F b dt
Voici comment on utilise ces formules. On choisit un pas de temps dt très petit. Par exemple dt=0.01. On discrétise le temps avec un pas de temps dt : la date de départ est t 0 . ensuite il y a t 1 = t 0 +dt , t 2 = t 1 +dt= t 0 +2dt , etc..
On suppose que a 0 =a( t 0 ) , b 0 =b( t 0 ) sont les conditions initiales connues. Avec les formules de Euler on calcule a( t 1 ) ,b( t 1 ) , puis a( t 2 ) ,b( t 2 ) , etc

5.4.3 Programmation

A chaque étape suggérée ci-dessous, tester votre programme.
Exercice 5.4.2. « Etape 1, sans graphisme ».
  1. Déclarer les paramètres suivants comme variables globales: g=1,m=1,l=1, γ =0, le pas de temps dt=1 0 -2 , conditions initiales ( a,b ) =( 1 0 -1 ,0 ) . Calculer et afficher T . Tester ce programme.
  2. Définir une fonction F( a,b ) qui en entrée prend a,b et renvoit da dt , db dt d'après les lois du mouvement (1). Tester cette fonction avec quelques valeurs.
  3. Faire une boucle sur t=0T par pas de temps dt et itèrer les formules de Euler (4). Afficher les valeurs de ( t,a,b ) à chaque instant. Vérifier à la lecture des résultats que le point ( a( t ) ,b( t ) ) a bien parcouru un tour dans le sens indirect.
Solution 5.4.3. pendule_etape_1.py
 
Exercice 5.4.4. « Etape 2, avec graphisme ». Le but est d'obtenir les animations ci-dessous. Modifier et compléter le programme précédent pour
  1. Dans la boucle précédente de la méthode de Euler, rajouter quelques lignes de code qui font le dessin du pendule dans le plan vertical x,z . Pour cela on utilisera la fonction plot() pour tracer un segment qui modélise la barre du pendule. Remarque: on fait le dessin du pendule à l'instant « présent » dans la boucle, il n'est donc pas utile d'utiliser une liste pour mémoriser les états passés.
  2. Rajouter le dessin du pendule dans l'espace des phases a,b au cours du temps.
  3. Sauvegarder sous forme movie.mp4 et movie.gif (gif animé).
  4. Essayer avec γ =0 et dt=0.01 . Quel phénomène lié à l'approximation de Euler observez vous?
image: 56_home_faure_enseignement_informatique_python_images_animation_pendule_euler.gif
image: 57_home_faure_enseignement_informatique_python_images_animation_pendule_euler_2.gif
Solution 5.4.5. pendule.py

5.4.4 Méthode de résolution numérique avec la fonction odeint()

Afin d'avoir une meilleur précision (et rapidité) numérique on va remplacer l'approximation de Euler précédente par une approximation appelée « Runge Kutta » et qui est déjà implémentée dans la fonction odeint() expliquée en Section 4.
Solution 5.4.6. pendule_odeint.py

5.5 Equation de la chaleur

Difficulté: moyenne.
On souhaite visualiser la progression du champ de température dans un matériau lorsque la température est imposée et fixée au bords. Dans le premier exemple ci-dessous en dimension 2, la température initiale est 0 (bleu) partout et on impose un gradient de température au bord (chaud (rouge) en haut). Dans le deuxieme exemple, on impose de plus un point froid (à 0) au centre du matériau:
image: 58_home_faure_enseignement_informatique_python_images_animation_equation_chaleur.gif image: 59_home_faure_enseignement_informatique_python_images_animation_equation_chaleur_2.gif
Notons x=( x 1 , x n ) les variables d'espace et t la variable de temps. Joseph Fourier a découvert en 1822 que la propagation de la température T( x,t ) dans un matériau est régit par l'équation d'évolution T( x,t ) t =( Δ T ) ( x,t ) appelée équation de la chaleur, avec l'opérateur Laplacien: Δ T:= 2 T x 1 2 + + 2 T x n 2
La valeur de la température est imposée sur les bords du domaine.
On va modéliser cette évolution en dimension n=2 sur le domaine carré x 1 , x 2 [ 0,1 ] 2 . On discrétise l'espace avec un pas h1 très petit. On note donc x 1 =ih , x 2 =jh avec i,j=0 N-1 . Le champ de températures à l'instant t est modélisé par une matrice T( i,j ) de taille N × N . On discrétise le temps avec un pas ε 1 très petit et on note T' le champ de température à l'instant t+ ε .
Exercice 5.5.1. (Mathématiques) On remplaçant les dérivées par des différences finies montrer que l'équation de la chaleur ci-dessus devient (au premier ordre en h, ε 1 ):
T'( i,j ) =T( i,j ) +K( T( i+1,j ) +T( i-1,j ) +T( i,j+1 ) +T( i,j-1 ) -4T( i,j ) ) . (5) avec K= ε h 2 . Sur les 4 bords, c'est à dire en ( 0, x 2 ) ,( x 1 ,0 ) ,( 1, x 2 ) ,( x 1 ,1 ) , on souhaite imposer un gradient de « température » T( x 1 , x 2 ) = x 1 δ T avec δ T fixé. Ecrire ce que cela impose pour la matrice T( i,j ) . Montrer que le processus est convergent si K= ε h 2 1 . Quelle est la solution attendue de l'équation stationnaire Δ T=0 ?.
Solution:
D'après la formule de Taylor au 1er ordre: T( x 1 +h, x 2 ) T( x 1 , x 2 ) + T x ( x 1 , x 2 ) h donc T x ( x 1 , x 2 ) 1 h ( T( x 1 +h, x 2 ) -T( x 1 , x 2 ) ) De même T x ( x 1 -h, x 2 ) 1 h ( T( x 1 , x 2 ) -T( x 1 -h, x 2 ) ) x ( T x ) ( x 1 -h, x 2 ) 1 h ( T x ( x 1 , x 2 ) - T x ( x 1 -h, x 2 ) ) = 1 h 2 ( T( x 1 +h, x 2 ) -2T( x 1 , x 2 ) +T( x 1 -h, x 2 ) ) et de même pour y ( T y ) . On déduit Δ T( x 1 , x 2 ) 1 h 2 ( T( x 1 +h, x 2 ) +T( x 1 -h, x 2 ) +T( x 1 , x 2 +h ) +T( x 1 , x 2 -h ) -4T( x 1 , x 2 ) ) On a de même T t ( x,t ) 1 ε ( T( x,t+ ε ) -T( x,t ) ) . Donc l'équation de la chaleur T( x,t ) t =( Δ T ) ( x,t ) donne 1 ε ( T( x,t+ ε ) -T( x,t ) ) = 1 h 2 ( T( x 1 +h, x 2 ) +T( x 1 -h, x 2 ) +T( x 1 , x 2 +h ) +T( x 1 , x 2 -h ) -4T( x 1 , x 2 ) ) T( x,t+ ε ) =T( x,t ) + ε h 2 ( T( x 1 +h, x 2 ) +T( x 1 -h, x 2 ) +T( x 1 , x 2 +h ) +T( x 1 , x 2 -h ) -4T( x 1 , x 2 ) ) , ce qui donne (5).
Conditions aux bords T( x 1 , x 2 ) = x 1 δ T . On a x 1 =i/( N-1 ) et x 2 =j/( N-1 ) avec i,j=0N-1 . Sur les bords ( x 1 ,0 ) et ( x 1 ,1 ) cette condition s'écrit T[ i,0 ] =i δ T/( N-1 ) ,T[ i,N-1 ] =i δ T/( N-1 ) . Sur le bord ( 0, x 2 ) cette condition s'écrit T[ 0,j ] =0 Sur le bord ( 1, x 2 ) cette condition s'écrit T[ N-1,j ] = δ T Le processus est convergent si K= ε / h 2 1 .
Exercice 5.5.2. (Programmation)
  1. Ecrire un programme qui définit N=4 (pour tester), h=1/N , δ T=10 et initialise la matrice T avec les valeurs aux bords et 0 à l'intérieur. Afficher la matrice et dessiner pour vérifier (Aide: Section 2.1.2.1). Le résultat devrait ressembler à:
    N= 4
    T=
    [[ 0. 0. 0. 0. ]
    [ 3.33333333 0. 0. 3.33333333]
    [ 6.66666667 0. 0. 6.66666667]
    [ 10. 10. 10. 10. ]]
    image: 60_home_faure_enseignement_informatique_python_images_equation_chaleur_1.png
  2. Écrire une fonction def Evolue(): qui fait évoluer tous les éléments du tableau T sur un pas de temps, d'après (5). Vérifier dans le cas N=4 et dessiner le tableau. Incorporer l'appel de cette fonction dans une boucle temporelle qui ne s'arrête que si le tableau n'évolue plus à 1 0 -4 près. Pour évaluer cela il faut mesure la norme T'-T L 1 = i,j | T'( i,j ) -T( i,j ) | (que l'on calcule au fur et à mesure dans la boucle sur i,j ) et que la fonction Evolue() renvoie cette norme. Vérifier que vous obtenez la solution stationnaire attendue. Le résultat devrait ressembler à cela:
    T=
    [[ 0. 0. 0. 0. ]
    [ 3.33333333 3.32288355 3.32288355 3.33333333]
    [ 6.66666667 6.65621662 6.65621662 6.66666667]
    [ 10. 10. 10. 10. ]]
    norme= 9.67313700277e-05
    image: 61_home_faure_enseignement_informatique_python_images_equation_chaleur_2.png
  3. Modifier votre programme pour voir une animation. Avec N=40 , cela devrait ressembler aux animations présentées au début. Essayer différentes géométries.

5.6 Zeros d'un polynome aléatoire

On considère un polynome P( x ) de degré N1 et à coefficients P k aléatoires et indépendants, selon une loi uniforme P k [ -1,1 ] : P( x ) = k=0 N P k x k
On sait que P a N zéros ( z i ) i=0N dans le plan complexe. Dessiner les N zéros du polynome et observez. Quelle conjecture avez vous? Consulter cet article de Terence Tao et al.. On peut aussi considérer des coefficients complexes aléatoires: P k = R k +i I k avec R k , I k [ -1,1 ] aléatoires et indépendants.
Aide: les zéros du polynome sont les valeurs propres de la “matrice compagnon” du polynome. Construire cette matrice et utiliser une fonction pour trouver les valeurs propres.
image: 62_home_faure_enseignement_informatique_c++_cours_c++_cours_1_zeros.png

Partie II Compléments du langage python

Chapitre 6 Classes et Objets

Programmer une classe consiste à définir un nouveau type d'objets avec des opérations précises que l'on pourra faire avec. Vous avez utilisé par exemple la classe des nombres complexes déjà existante (écrite par d'autres informaticiens avant vous). Avec cette classe, on peut déclarer des variables complexes et effectuer directement des opérations sur les nombres complexes comme z3 = z1*z2*exp(z4),...sans devoir décomposer l'opération sur les parties réelles et imaginaires.

Exemple

Construisons par exemple une classe d'objets appelée Vecteur, et décrivant un déplacement dans le plan x,y .
import math
import matplotlib.pyplot as plt
 
class Vecteur:   # Déclaration de la classe
    x,y = 4,3        # des variables membre de la classe, pour stocker les coordonnées
 
    def Affiche(self): # fonction membre qui affiche les coordonnées
        print ('x=', self.x, 'y=', self.y)
 
    def Dessin(self): # fonction membre qui dessine le vecteur
        plt.plot([0, self.x], [0, self.y], marker='o')
        plt.show()
 
    def Longueur(self): # fonction membre qui renvoie la longueur
         self.L = math.sqrt(self.x**2 + self.y**2)    # rajout d'une variable membre
         return self.L
 
v = Vecteur()                 # Création d'un objet de type Vecteur
v.Affiche()                 # affiche ses coordonnées
v.Dessin()                 # dessin du vecteur
print ('L=', v.Longueur())      # affiche sa longueur
v.x = 0                  # modification de x
v.Affiche()                 # affiche ses coordonnées
print ('L=', v.Longueur())
 
résultat:
x=4 y=3
image: 63_home_faure_enseignement_informatique_python_images_line.png
L=5.0
x=0 y=3
L=3.0
Remarque 6.1.1. Le mot-clef “self” : désigne l'objet lui-même dans le code de la classe. Par exemple, “self.x” veut dire “la variable x dans l'objet Vecteur où le code est exécuté” (si il s'agit du Vecteur v, cela sera v.x).
 
Remarque 6.1.2. On peut ajouter une (des) variable(s) membre après avoir déclaré une classe. La variable membre est alors accessible dans tous les objets, même ceux créé avant. A la suite du code ci-dessus, écrire
Vecteur.z = 1   # Ajout de la variable membre 'z'
w = Vecteur()
print (w.z, w.x, v.z, v.x)
 
résultat
1 4 1 0

6.1.1 Fonctions membres spéciales

Certaines fonctions membre ont des noms prédéfinis qui commence et finit par deux traits soulignés. Ces fonctions sont appelées dans des situations spécifiques.

6.1.1.1 Constructeur __init__

Cette fonction appelée « constructeur » est appelée à l'initialisation d'un nouvel objet. Ecrire:
import math
class Vecteur:
    def __init__(self, x0, y0):   # fonction constructeur
        self.x, self.y = x0, y0
 
    def Affiche(self): # fonction membre qui affiche les coordonnées
        print ('x=', self.x, 'y=', self.y)
 
w = Vecteur(4, 3)                 # Création d'un objet avec le constructeur
w.Affiche()                 # affiche ses coordonnées
 
résultat:
x=4 y=3

6.1.1.2 Destructeur __del__

Le destructeur est la fonction qui est appelée lors de la destruction de l'objet. Il ne prend qu'un argument (self). Il peut servir à 'nettoyer' des données créées par l'objet.

6.1.1.3 Opérateurs __add__, __sub__, __mul__, __div__, __pow__

Au moment de leur utilisation ces fonctions sont appelées par les symboles respectivement (+,  -,  *,  /,  **). Ces fonctions prennent deux arguments, le premier qui est 'self' et un deuxième désignant le second membre de l'opération. Voici un exemple avec la fonction __add__()
import math
class Vecteur:
    def __init__(self, x0, y0):   # fonction constructeur
        self.x, self.y = x0, y0
 
    def Affiche(self): # fonction membre qui affiche les coordonnées
        print ('x=', self.x, 'y=', self.y)
 
    def __add__(self, w):  # addition de deux vecteurs
        return Vecteur( self.x + w.x, self.y + w.y)
 
v1 = Vecteur(1,2)
v2 = Vecteur(3,4)
v3 = v1 + v2     # cela appelle la fonction __add__(v1,v2)
v3.Affiche()
résultat
x=4, y=6

6.1.2 Héritage de classes

L'héritage permet de créer des classes qui héritent des propriétés de la (les) classe(s) parente(s).On peut ajouter des propriétés spécifiques à la classe « enfant ». L'héritage peut être très utile lorsque l'on génère de multiples classes ayant des propriétés similaires. Cela peut être aussi utile pour ré-utiliser du code, provenant par exemple d'une librairie développée par quelqu'un d'autre, en ajoutant quelques fonctions à une classe pré-existante.
Voici un exemple où la classe Segment hérite de Vecteur et rajoute l'information du point de départ.
import math
import matplotlib.pyplot as plt
 
class Vecteur:   # Déclaration de la classe
    def __init__(self, x0, y0):   # fonction constructeur
        self.x, self.y = x0, y0
 
    def Dessin(self): # fonction membre qui dessine le vecteur
        plt.plot([0, self.x], [0, self.y], marker='o')
        plt.show()
 
class Segment(Vecteur):        # Hérite de la classe Vecteur
    def __init__(self, xO, yO, x1, y1, col0):
        self.xO, self.yO = xO,yO # origine
        self.x, self.y = x1, y1 # deplacement
        self.col = col0  #couleur
 
    def Dessin(self): # fonction membre qui dessine le segment
        plt.plot([self.xO, (self.xO + self.x)], [self.yO, (self.yO + self.y)], self.col, marker='o')
        plt.show()
 
s1 = Segment(0,0, 2, 1, 'blue')
s1.Dessin()
s2 = Segment(2, 1, 2, -1 ,'red')
s2.Dessin()
résultat:
image: 64_home_faure_enseignement_informatique_python_images_segment1.png image: 65_home_faure_enseignement_informatique_python_images_segment2.png

6.1.3 Copie d'objets

Ecrire:
import copy
a = Vecteur(2, 1)
b = a # copie de l'adresse
c = copy.deepcopy(a) # copie de l'objet
a.x = 3
print('a.x=', a.x , ' b.x=', b.x, 'c.x=', c.x)
print( a==b, a==c)
résultat:
a.x= 3 b.x= 3 c.x= 2
True False
Remarque 6.1.3. Lorsque l'on écrit b=a, alors b et a deviennent deux noms différents du même objet (même adresse mémoire). Lorsque l'on écrit c = copy.deepcopy(a), alors c et a sont des objets différents mais contenant les même données. La commande c = copy.copy(a) effectue aussi une copie mais non récursive.

Chapitre 7 Analyse de données avec numpy

Curve fit

plot_curve_fit.html

Chapitre 8 Quelques modules utiles

8.1 Commandes du systèmes d'exploitation

voir module os.

8.1.1 Exemple

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Programme python pour creation/destruction de repertoires.
"""
Instructions
===========
  - sauvegarder ce programme dans un repertoire, par exemple dans  /home/nadia/python

   -Ouvrir un terminal.

   - Se placer dans le repertoire voulu, où tu veux créer/détruire des répertoires,  appelé "repertoire courant", par la commande:

         cd repertoire

   - Lancer le programme avec des informations supplementaires par exemple:
            python3 /home/nadia/python/creation_repertoires.py NOM 2 5

   - Sortie:
      il cree les repertoires NOM_0002 NOM_0003 ... NOM_0005
      dans le repertoire courant.

   - option pour detruire ces repertoires (si ils existent) rajouter l'option "destroy" a la fin par exemple:

         python3 /home/nadia/python/creation_repertoires.py NOM 2 5 destroy

         Attention: il detruira les repertoires  NOM_0002 ... NOM_0005, meme si ils contiennent des donnees.

 Remarques:
   - on peut modifier ces programmes pour changer des options, par exemple le nombre de zeros, etc..




References:
----------
sur le module os:
 https://docs.python.org/fr/3.6/library/os.html

sur le formattage des entiers:
   https://pyformat.info/

Informations utiles sur les commandes dans un terminal:

  https://doc.ubuntu-fr.org/tutoriel/console_commandes_de_base
"""

#--- importation des librairies que l'on utilise
import os
import sys

#--- en entree sys.argv est une liste avec les parametres donnes 
# print(sys.argv)

if (len(sys.argv) < 4): # teste la longueur de la liste de parametres
   print('il manque des arguments. Recommencez.')
   exit()
elif (len(sys.argv) > 5):
   print('il y a trop d arguments. Recommencez.')
   exit()

#--- tout semble ok on recupere les arguments
nom = sys.argv[1]
N1 = int(sys.argv[2])
N2 = int(sys.argv[3])

option='create'

if (len(sys.argv) == 5):
    if (sys.argv[4] != 'destroy'):
        print('l argument 4 devrait etre: destroy? Recommencez.')
        exit()
    else:   
        option = 'destroy'

# -- recapitulatif.
if (option == 'create'):
    print("Je vais creer les repertoires: "+ nom +  '_{:04d}'.format(N1)  + " ... " + nom +  '_{:04d}'.format(N2) )
elif (option == 'destroy'):
    print("Je vais detruire les repertoires: "+ nom +  '_{:04d}'.format(N1)  + " ... " + nom +  '_{:04d}'.format(N2) )

# ---- Effectue le travail --------------------
for i in range(N1, N2+1):
    nom_rep = nom +  '_{:04d}'.format(i)
    if(option == 'create'):
        try:
            os.mkdir(nom_rep)
        except:
            print('Le repertoire '+ nom_rep +' existe deja.')
    elif (option == 'destroy'):
        try:
            os.rmdir(nom_rep)
        except:
            print('Le repertoire '+ nom_rep +' ne peut etre detruit car il est non vide.')

print('fini.')

8.2 Envoyer un mail

 
import smtplib
server = smtplib.SMTP('smtp.free.fr')
expediteur='toto@free.fr'
destinataire='toto@free.fr'
Message='''\
From: (Regis) regis@free.fr
To: (Toto) toto@free.fr
Subject: Profil
 
  J'aime bien ton profil!'''
# Message (il faut laisser une ligne blanche avant le contenu message # et faire attention a la marge
try:
  server.sendmail(expediteur,destinataire,Message)
  server.quit()
except SMTPException as why:
  print (why)

8.3 Notation musicale

Abjad
Voir cette liste
JythonMusic

8.4 Analyse de signaux audio (sur fichiers, temps différé)

Référecnes mathématiques:
Reference informatiques:
Some examples:

8.4.1 Play/record sound

List of python modules for audio

8.4.2 Wavelett transform with scipy

Morlet2 wavelet:
Reference on morlet2.
For x R ,
ψ s,w ( x ) := 1 s π 1/4 e i w s x e - 1 2 ( x s ) 2 Sampling. Let
x j =j- M-1 2 So that if M is even, x j =- M-1 2 , - 1 2 ,+ 1 2 , + M-1 2 and if M is even, x j =- M-1 2 , -1,0,1, + M-1 2 .
We get a list ( ψ j ) j=0M-1 by Ψ M,s,w ( j ) := ψ s,w ( x j )
Mathematical wave-packet
ϕ t, ω ( t' ) =a e i ω t' e - 1 2 ( t'-t σ ) 2 with a= ( π σ 2 ) - 1 4 so that ϕ t, ω L 2 ( R ,dt' ) =1 .
Sampling: for M N , δ >0 , t 0 = M-1 2 δ , we sample Φ M, δ ,t, ω , σ ( j ) = ϕ t, ω ( t j ) with j=0M-1 and t j = δ j- t 0 so that t j =- t 0 t 0 with steps δ and
Relation
We have Φ M, δ ,t, ω , σ ( j ) =C Ψ M,s,w ( j ) if t j = δ x j t=0,s= σ δ ,w= σ ω ,C= s σ conversely take δ =1 .
Preuve. We have ω t j = w s x j ω δ = w s and x j s = t j σ 1 s = δ σ s= σ δ hence w=s ω δ = σ ω
and C s π 1/4 = ( π σ 2 ) - 1 4 C= s σ
image: 66_home_faure_enseignement_informatique_python_images_wave_packet.png
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Jan  2 09:19:29 2021

@author: faure
"""
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

#===================
#wave packet with morlet2

f = 300 # frequency
sigma = 0.01 # width 
omega = 2 *np.pi * f



delta = 1/44100.  # time step
t0 = 3*sigma # time interval
M = int(2*t0/delta) # size of list



s = sigma/delta 
w =   sigma*omega
wavelet = signal.morlet2(M, s, w)


Lt = np.linspace(-t0, t0, M)
plt.plot(Lt, np.real(wavelet))
plt.xlabel('$t$')
plt.ylabel('$\phi$')
plt.show()

8.4.3 Plot a file in wav format

From this wav file that was recorded with audacity.
image: 67_home_faure_enseignement_informatique_python_images_signal.png
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 23 17:03:59 2020

@author: faure

Plot a wav file

"""

import math
import wave
import numpy as np
import matplotlib.pyplot as plt

from scipy.io import wavfile


#=============================================
# using the module wave
# input: name of the file
# output: array with the data
def file_wav_to_array(name_file):
        
    global n_chan, data_size, rate, N
    # Read wav file 
    f = wave.open(name_file)
    #f = wave.open("offrandes_oubliees.wav")
    
    
    
    #informations on the sound
    N = f.getnframes()  # number of data
    n_chan = f.getnchannels()   # 1: mono, 2 : stereo
    data_size = f.getsampwidth()  # in octet: ex: 2
    rate = f.getframerate() # ex: rate = 44100 Hz
    
    
    # transfer the file to a list
    L_u = f.readframes(N)
    
    
    # Convert buffer to float32 using NumPy                                                                                 
    buf_as_np_int16 = np.frombuffer(L_u, dtype = np.int16)
    buf_as_np_float32 = buf_as_np_int16.astype(np.float32)
    
    # Normalise float32 array so that values are between -1.0 and +1.0                                                      
    max_int16 = 2**15
    return buf_as_np_float32 / max_int16, N, rate


#======================================
# plot the signal
# input: array with the data
def plot_signal(L_u):    
    N = len(L_u)
    t_max =  N /rate
    
    t = np.linspace(0, t_max, N)  # N values equidistributed in [0,tmax]
    

    plt.figure() # open a new window    
    plt.cla() # efface le graphisme precedent
    plt.plot(t, L_u) 
    #plt.plot( L_u)  # plot with data in abcisee
    plt.xlabel('t(s):temps')
    plt.ylabel('u(t)')
    plt.title('Signal u')
    plt.grid() # montre grille
    plt.axis([0, 0.01, -1,1]) # selectionne la vue
    plt.show() # montre le dessin
    

#========================================

#---- read wavefile using the module wave
#L_u, N, rate = file_wav_to_array("audio.wav") # L_u: list
#L_u, N, rate = file_wav_to_array("trumpet.wav") # L_u: list

#---- read wavefile using the module wavfile
rate, L_u = wavfile.read("audio.wav")



plot_signal(L_u)


D = 2 # duration
F= 44100 # frequency of sampling in Hz
dt = 1/F # time step
L_t = np.arange(0, D, dt) # list of discrete times
f0= 440
L_u = 0.5* np.cos(2.0 * np.pi *   f0 * L_t + np.pi/4) 

#---- write  wavefile using the module wavfile. The maximale amplitude is 1
wavfile.write('audio.wav', F, L_u)

8.4.4 Detect the pitch from a file in wav format

From this wav file that was recorded with audacity.
image: 67_home_faure_enseignement_informatique_python_images_signal.png image: 68_home_faure_enseignement_informatique_python_images_pitch.png
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 23 17:03:59 2020

@author: faure

Input: wav file

Ouput: pitch detect

"""


import math
import wave
import numpy as np
import matplotlib.pyplot as plt
from numpy import linalg as LA



f_A = 440 # freq of A in Hz

#=============================
def freq_to_pitch(f):
    return  12/math.log(2) * np.log(f/f_A) + 69 


#=============================
def pitch_to_freq(x):
    return  f_A * np.exp(math.log(2.)/12. *  (x - 69))
 

#=============================================
# using the module wave
# input: name of the file
# output: array with the data
def file_wav_to_array(name_file):
        
    global n_chan, data_size, rate, N
    # Read wav file 
    f = wave.open(name_file)
    #f = wave.open("offrandes_oubliees.wav")
    
    
    
    #informations on the sound
    N = f.getnframes()  # number of data
    n_chan = f.getnchannels()   # 1: mono, 2 : stereo
    data_size = f.getsampwidth()  # in octet: ex: 2
    rate = f.getframerate() # ex: rate = 44100 Hz
    
    
    # transfer the file to a list
    buf = f.readframes(N)
    
    
    # Convert buffer to float32 using NumPy                                                                                 
    buf_as_np_int16 = np.frombuffer(buf, dtype = np.int16)
    buf_as_np_float32 = buf_as_np_int16.astype(np.float32)
    
    # Normalise float32 array so that values are between -1.0 and +1.0                                                      
    max_int16 = 2**15
    return buf_as_np_float32 / max_int16, N, rate



#======================================
# plot the signal
# input: array with the data
def plot_signal(buf):    
    N = len(buf)
    t_max =  N /rate
    
    t = np.linspace(0, t_max, N)  # N values equidistributed in [0,tmax]
    

    plt.figure() # open a new window    
    plt.cla() # efface le graphisme precedent
    plt.plot(t, buf) 
    #plt.plot( buf)  # plot with data in abcisee
    plt.xlabel('t(s):temps')
    plt.ylabel('u(t)')
    plt.title('Signal u')
    plt.grid()
    plt.show() # montre le dessin
    
    
#=========================================
#Compute pitch
def Compute_pitch():
    itau_max = 300 # size max for the shift
    fmin =  rate/itau_max
    print("Lower frequency is :", fmin , " Lower pitch is:",   freq_to_pitch(fmin))
    N2 = int(N/itau_max) 
    C_seuil = 0.4 #  seuil pour C
    n2_seuil = itau_max * 1e-4 # seuil for norme n2
    L_f = np.zeros(N2)
    
    t_s = 1.8 # time selection for display 
    it2_s = int(t_s*rate /itau_max)
    L_C_s = np.zeros(itau_max) # for display only
    
    for it2 in range(0, N2 -2): # interval start
        #print(N2-it2)
        it = it2*itau_max
        L1 = buf[it : it + itau_max] #extract
        L3 = buf[it : it + 2*itau_max] #extract
        n2 = LA.norm(L3)
        C = 4;
        itau_min = 0;
        
        for itau in range(10, itau_max):  # shift
            L2 = buf[it+itau : it + itau + itau_max] #extract
            if(n2 > n2_seuil):
                C2 = LA.norm(L2-L1) / n2
                if(it2 == it2_s): # for display
                    L_C_s[itau] = C2
                 
                if(C2 < C):
                    C=C2
                    itau_min = itau
                    
             
    
        if(C < C_seuil):
            freq = rate/itau_min
            L_f[it2] =  freq_to_pitch(freq)
        
        
        
        
    #====== plot C
            
    plt.figure() # open a new window
    plt.cla() # efface le graphisme precedent
    plt.plot(L_C_s,'b.')
    plt.xlabel('tau(s):time shift')
    plt.ylabel('C(tau)')
    plt.title('Difference')
    
    
    #====== plot the pitch view
            
    plt.figure() # open a new window
    Lt2 = np.linspace(0, t_max, N2) # N2 values equidistributed in [0,tmax]
    plt.cla() # efface le graphisme precedent
    plt.plot(Lt2, L_f, 'r.') # red and points
    plt.xlabel('t(s):temps')
    plt.ylabel('pitch')
    plt.title('Pitch of the signal')
    
    
    plt.grid() # show a grid
    
    
 
    xmin = 60
    xmax = 73
    plt.axis([0, t_max, xmin, xmax]) # select the view
    
    plt.yticks(numpy.arange(xmin, xmax, 1)) # impose the list of labels
    ay = plt.gca() #  get current axis
    a = ay.get_yticks().tolist() # get list of current axis
    a = ['C5', 'C#', 'D', 'D#', 'E', 'F', 'F#', 'G', 'G#', 'A', 'A#', 'B', 'C6']
    ay.set_yticklabels(a)
    
     
    plt.show() # montre le dessin
    
        
#=============================================

buf = file_wav_to_array("flute.wav")

plot_signal(buf)

Compute_pitch()

8.4.5 Here is the code to draw a waveform and a frequency spectrum of a wavefile

Reference about specgram.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 18 17:57:04 2020

@author: faure
"""

import wave
import numpy as np
import matplotlib.pyplot as plt

signal_wave = wave.open('flute.wav', 'r')
sample_rate = 44100   # en Hz
sig = np.frombuffer(signal_wave.readframes(sample_rate), dtype=np.int32)

sig = sig[:] #For the whole segment of the wave file

#sig = sig[25000:32000] #For partial segment of the wave file

#left, right = data[0::2], data[1::2] #Separating stereo channels

plt.figure(1)

plot_a = plt.subplot(211)
plot_a.plot(sig)
plot_a.set_xlabel('sample rate * time')
plot_a.set_ylabel('energy')

plot_b = plt.subplot(212)
plot_b.specgram(sig, NFFT=1024, Fs=sample_rate, noverlap=900)
plot_b.set_xlabel('Time')
plot_b.set_ylabel('Frequency')

plt.show()

8.5 Audio en temps réel

Voir PyAudio, doc

8.5.1 Exemple d'enregistrement en temps réel.

import pyaudio
import numpy as np

CHUNK = 2**10
RATE = 44100

p=pyaudio.PyAudio()
stream=p.open(format=pyaudio.paInt16,channels=1,rate=RATE,input=True, frames_per_buffer=CHUNK)

for i in range(int(10*RATE/CHUNK)): #go for a few seconds
    data = np.fromstring(stream.read(CHUNK),dtype=np.int16)
    peak=np.average(np.abs(data))*2
    bars="#"*int(50*peak/2**16)
    print("%04d %05d %s"%(i,peak,bars))

stream.stop_stream()
stream.close()
p.terminate()
Avec dessin:
import pyaudio
import numpy as np
import matplotlib.pyplot as plt

CHUNK = 2**10 # 1024 buffer size
RATE = 44100  # frequency of sampling

p = pyaudio.PyAudio() #  nouvel objet p de la class PyAudio
stream = p.open(format = pyaudio.paInt16, channels=1, rate=RATE, input=True, frames_per_buffer = CHUNK)

for i in range(1000):
    plt.cla() # efface le graphisme precedent
    T = np.fromstring(stream.read(CHUNK), dtype=np.int16)
    plt.plot(T)
    plt.axis([0, 1024, -32000, 32000]) # selectionne la vue
    plt.pause(0.001) # montre la figure et attend 0.001 sec.

stream.stop_stream()
stream.close()
p.terminate()
Utilisant La capture audio dans un thread séparé:
#May 2020
import pyaudio
import numpy as np
import matplotlib.pyplot as plt
import time

CHUNK = 2**11 # 2**11 = 2048  buffer size
RATE = 44100  # frequency of sampling



p = pyaudio.PyAudio() #  nouvel objet p de la class PyAudio

"""
 this function is called on a separate thread for each new buffer on audio input
 input: in_data: the input data
         frame_count: taille du tableau = CHUNK
"""
def audio_event(in_data, frame_count, time_info, status):
    global T
    T = np.fromstring(in_data, dtype=np.int16) # transforme en tableau de chiffres
    return (in_data, pyaudio.paContinue)
  

#----------------------------------------------
stream = p.open(format = pyaudio.paInt16, channels=1, rate=RATE, input=True, frames_per_buffer = CHUNK, stream_callback = audio_event)


stream.start_stream() # lance la fonction audio_event() dans un autre thread


while stream.is_active():
    plt.cla() # efface le graphisme precedent
    time.sleep(0.1)
    plt.plot(T)
    plt.axis([0, 1024, -32000, 32000]) # selectionne la vue
    plt.pause(0.001) # montre la figure et attend 0.001 sec.
   
    
stream.stop_stream()
stream.close()
p.terminate()

8.6 Acquisition avec PyVisa

PyVISA: Control your instruments with Python
PyTektronixScope

Chapitre 9 Librairie PyRoot pour le graphisme et l'analyse de données

Références sur PyRoot

Chapitre 10 Librairie opencv pour l'analyse d'images

Voir opencv en python. Tutorial
Remarque 10.0.1. Pour vérifier que opencv est bien installé, écrire (en python 2):
import cv2 as cv
print(cv.__version__)
résultat:
2.4.9.1

10.1 Lecture et ecriture de fichiers images

Les commentaires expliquent les fonctions. Il faut avoir le fichier Image.jpg dans le répertoire du programme.
Ecrire:
 
import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
  
#..... Lecture d'une image -> img
img = cv.imread('Image.jpg')
 
#..... dessin de l'image img en utilisant opencv
cv.imshow('Titre',img)
 
#..... dessin de l'image img en utilisant matplotlib
plt.imshow(img, cmap = 'gray', interpolation = 'bicubic')
plt.xticks([]), plt.yticks([]) # to hide tick values on X and Y axis
plt.show()
 
#...... Ecriture d'une image: img -> fichier Image2.jpg
cv.imwrite('Image2.jpg',img)
  
#..... Lecture d'un pixel
x=img.shape[0]/2 # selectionne point (x,y) au centre de l'image
y=img.shape[1]/2
p= img[x,y] # valeur du pixel
print('La valeur du pixel en x=',x,', y=',y,' est (B,G,R): (',p,')')
 
#..... Attend qu'une touche soit appuyee.
c = cv.waitKey(0) # rem: montre le dessin. Renvoit le code de la touche. On peut mettre une duree t en ms. ex: 2000
Resultat:

10.2 Video

Tutorial.
Ecrire:
import numpy as np
import cv2
 
cap = cv2.VideoCapture(0)
 
while(True):
  # Capture frame-by-frame
  ret, frame = cap.read()
  # Our operations on the frame come here
  gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
  # Display the resulting frame
  cv2.imshow('frame',gray)
  if cv2.waitKey(1) & 0xFF == ord('q'):
    break
 
# When everything done, release the capture
cap.release()
cv2.destroyAllWindows()

10.3 Références sur Opencv with python and deep-learning

2019: On suit essentiellement les recommandations de www.pyimagesearch.com, voir aussi guide.

10.3.1 Installations

10.3.2 Exemples

10.3.3 Object detection with pretrained models

10.3.3.1 with YOLO algo

10.3.3.2 with opencv::dnn module and SSD (Single Shot Detectors) and MobileNets algos.

10.3.4 Quelques articles d'algorithmes de deep learning

10.3.5 DataBase

10.3.6 Utiliser Coral de Google

10.3.7 Différents modèles

10.3.8 Create a GUI for opencv

10.3.9 Sites web

Chapter 11 Python virtual environments

permet de travailler avec des versions spécifiques pour un projet donné.
Creation d'un envirronement appelé cv:
mkvirtualenv cv -p python3
Lancer l'environement:
workon cv
Se deconnecter:
deactivate
On peut installer des librairies, par ex:
pip install numpy

Chapitre 12 Librairie Caffe pour le deep Learning (DL)

12.1 Deep learning

Ref: article.

12.1.1 The model

Let σ : R R the function « sigmoid »: σ ( x ) := 1 1+ e -x so that σ ( - ) =0 , σ ( ) =1 .
We extend σ to σ : R n R n (keeping the same notation) by σ ( x ) :=( σ ( x 1 ) , , σ ( x n ) ) .
Let A: R m R n an affine map, i.e. Ax=Wx+b. with b R n and W: R m R n linear.
Définition 12.1.1. A « Layer » is a map L: R m R n of the form L( x ) = σ ( Ax ) . A d- layer map is F: R n 0 R n d-1 given by F= L d-1 L 0 with layer maps L j : R n j R n j+1 .

12.1.1.1 Training datas

Suppose that we have a set of k « data points » D:={ ( x 1 , y 1 ) , ,( x k , y k ) } with x j R n 0 in the initial space and y j R n d-1 in the final space, called the « answer ». Usually y j { A,B } can have two values.

12.1.1.2 Cost function

Définition 12.1.2. For a given map F as above and training datas D , we associate the « cost » or « objective » C( F,D ) := j=1 k y j -F( x j ) l 2 2 C j ( F ) R + . (6)

12.1.1.3 Training problem or optimization problem

For given datas D , find the map F min (in a connected component, i.e. only parameters of F ) that minimizes C( F,D ) : F min :={ F  s.t.  C( F min ,D ) = min F C( F,D ) } .

12.1.2 Algorithm of gradient flow

A problem is to compute F min .
Remarque 12.1.3. in matlab, there is the function lsqnonlin
Gradient flow is the flow φ t :MM generated by the gradient vector field X=-( C F ) . C on the space M={ F } of all F maps, equation of motion is dF dt =-( C F ) ( F ) =- j=1 k ( C j F ) ( F )

12.1.2.1 Computation of ( C j F ) ( F )

We have C j ( F ) = y j -F( x j ) l 2 2 with F= L d-1 L 0 and L α ( x ) = σ ( A α x ) , A α x= W α x+ b α . Parameters of F are W α , b α with α =0 d-1 .
Hence C j W α ( F ) =2( y j -F( x j ) ) .( - F W α ) F W α = L d-1 L 0 with D L α ( x ) = σ 'D A α

12.1.2.2 « stochastic gradient » method,

Numerically it is usual to apply the « stochastic gradient » method, which consists in:
  1. discretize time in small steps dt .
  2. Repeat many times:
    1. choose some j{ 1, k } randomly and compute dF=-dt( C j F ) ( F )
    2. apply the modification FF+dF.

Chapitre 13 Deep learning

fish detection
fish_detection
fish detection

13.1 Using Coral USB Accelerator with tensorflow Lite

Chapitre 14 Convertir python en javascript

Documents: