Analyse scientifique avec Python

Yannick Copin

Ce cours est une initiation au langage Python et à son utilisation dans le domaine de l'analyse scientifique (manipulation de données numériques). Il a avant tout été développé comme support pour des enseignements en présentiel de 20h environ, et ne se veut donc pas complet : beaucoup d'éléments sont développés à l'oral, notamment les solutions des exercices.

Pour réagir au contenu de ce tutoriel, un espace de dialogue vous est proposé sur le forum. Commentez

- Titre : Analyse scientifique avec Python

- Auteur : Yannick Copin

- Parution : 18 septembre 2017

- Licence : Copyright ® 2017 Yannick Copin. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.

Sommaire

- I. Introduction

- I-A. Pourquoi un module d'analyse scientifique ?

- I-B. Pourquoi Python ?

- I-C. Index

- II. Initiation à Python

- II-A. Introduction

- II-A-1. Installation

- II-A-2. Notions d'Unix

- II-A-3. L'invite de commande

- II-B. Types de base

- II-C. Structures de programmation

- II-D. Les chaînes de caractères

- II-D-1. Indexation

- II-D-2. Sous-liste (slice)

- II-D-3. Méthodes

- II-D-4. Formatage

- II-E. Objets itérables

- II-F. Fonctions

- II-G. Bibliothèques et scripts

- II-G-1. Bibliothèques externes

- II-G-2. Bibliothèques personnelles et scripts

- II-H. Exceptions

- II-I. Classes

- II-J. Entrées-sorties

- II-J-1. Intéractif

- II-J-2. Fichiers

- II-K. Éléments passés sous silence

- II-K-1. Python 3.x

- III. Bibliothèque standard

- III-A. Gestion des arguments/options de la ligne de commande

- III-B. [c]Pickle : sérialisation des données

- III-C. Batteries included

- III-D. Text/Graphical User Interfaces

- IV. Bibliothèques numériques de base

- IV-A. Numpy

- IV-A-1. Tableaux

- IV-A-1-a. Création de tableaux

- IV-A-1-b. Manipulations sur les tableaux

- IV-A-1-c. Opérations de base

- IV-A-2. Tableaux évolués

- IV-A-3. Entrées/sorties

- IV-A-4. Sous-modules

- IV-A-5. Performances

- IV-B. Scipy

- IV-B-1. Tour d'horizon

- IV-B-2. Quelques exemples complets

- IV-C. Matplotlib

- IV-C-1. pylab vs. pyplot

- IV-C-2. Figure et axes

- IV-C-3. Sauvegarde et affichage interactif

- IV-C-4. Anatomie d'une figure

- IV-C-5. Visualisation 3D

- V. Bibliothèques scientifiques avancées

- V-A. Pandas et xarray

- V-A-1. Structures

- V-A-2. Accès aux données

- V-A-3. Manipulation des données

- V-A-4. Regroupement et agrégation de données

- V-A-5. Visualisations

- V-A-6. Xarray

- V-B. Astropy

- V-B-1. Tour d'horizon

- V-B-2. Démonstration

- V-B-2-a. Fichiers FITS

- V-B-2-a-i. Lire un fichier FITS

- V-B-2-a-ii. World Coordinate System

- V-B-2-a-iii. Visualisation dans matplotlib

- V-B-2-a-iv. Tables

- V-B-2-a-v. Quantités et unités

- V-B-2-a-vi. Calculs cosmologiques

- V-C. Autres bibliothèques scientifiques

- VI. Développer en Python

- VI-A. Le zen du Python

- VI-A-1. Us et coutumes

- VI-A-2. Principes de conception logicielle

- VI-B. Développement piloté par les tests

- VI-C. Outils de développement

- VI-C-1. Integrated Development Environment

- VI-C-2. Vérification du code

- VI-C-3. Débogage

- VI-C-4. Profilage et optimisation

- VI-C-5. Documentation

- VI-C-6. Python packages

- VI-C-7. Système de gestion de versions

- VI-C-8. Intégration continue

- VI-D. Python 2 vs. Python 3

- VII. Références supplémentaires

- VII-A. Documentation générale

- VII-B. Listes de liens

- VII-C. Livres libres

- VII-D. Cours en ligne

- VII-D-1. Python

- VII-D-2. Scientifique

- VII-D-3. Snippets

- VIII. Exemples

- VIII-A. Mean power (fonction, argparse)

- VIII-B. Formes (POO)

- VIII-C. Cercle circonscrit (POO, argparse)

- VIII-D. Matplotlib

- VIII-D-1. Figure (relativement) simple

- VIII-D-2. Filtres du 2d ordre

- IX. Exercices

- IX-A. Introduction

- IX-A-1. Intégration : méthode des rectangles ★

- IX-A-2. Fizz Buzz ★

- IX-A-3. PGCD : algorithme d'Euclide ★★

- IX-A-4. Tables de multiplication ★

- IX-B. Manipulation de listes

- IX-B-1. Crible d'Ératosthène ★

- IX-B-2. Carré magique ★★

- IX-C. Programmation

- IX-C-1. Suite de Syracuse (fonction) ★

- IX-C-2. Flocon de Koch (programmation récursive) ★★★

- IX-C-3. Jeu du plus ou moins (exceptions) ★

- IX-C-4. Animaux (POO/TDD) ★

- IX-C-5. Jeu de la vie (POO) ★★

- IX-D. Manipulation de tableaux (arrays)

- IX-D-1. Inversion de matrice ★

- IX-D-2. Median Absolute Deviation ★

- IX-D-3. Distribution du pull ★★★

- IX-E. Méthodes numériques

- IX-E-1. Quadrature et zéro d'une fonction ★

- IX-E-2. Schéma de Romberg ★★

- IX-E-3. Méthode de Runge-Kutta ★★

- IX-F. Visualisation (matplotlib)

- IX-F-1. Quartet d'Anscombe ★

- IX-F-2. Diagramme de bifurcation : la suite logistique ★★

- IX-F-3. Ensemble de Julia ★★

- IX-G. Mise en œuvre de l'ensemble des connaissances acquises

- IX-G-1. Équation différentielle ★

- IX-H. Exercices en vrac

- X. Solutions aux exercices

- X-A. Méthode des rectangles

- X-B. Fizz Buzz

- X-C. Algorithme d'Euclide

- X-D. Crible d'Ératosthène

- X-E. Carré magique

- X-F. Suite de Syracuse

- X-G. Flocon de Koch

- X-H. Jeu du plus ou moins

- X-I. Animaux

- X-J. Particules

- X-K. Jeu de la vie

- X-L. Median Absolute Deviation

- X-M. Distribution du pull

- X-N. Calculs numériques (numerique.ipynb)

- X-N-1. Quadrature

- X-N-2. Zéro d'une fonction

- X-O. Quartet d'Anscombe

- X-P. Suite logistique

- X-Q. Ensemble de Julia

- X-R. Trajectoire d'un boulet de canon (canon.ipynb)

- XI. Note de la rédaction de Developpez.com

I - Introduction

I-A - Pourquoi un module d'analyse scientifique ?

- Pour générer ses données, p.ex. simulations numériques, contrôle d'expériences.

- Pour traiter ses données, i.e. supprimer les artefacts observationnels.

- Pour analyser ses données, i.e. extraire les quantités physiques pertinentes, p.ex. en ajustant un modèle.

- Pour visualiser ses données, et appréhender leur richesse multidimensionnelle.

- Pour présenter ses données, p.ex. générer des figures prêtes à publier.

Ce module s'adresse donc avant tout aux futurs expérimentateurs, phénoménologistes ou théoriciens voulant se frotter à la réalité des observations.

I-B - Pourquoi Python ?

Les principales caractéristiques du langage Python :

- syntaxe simple et lisible : langage pédagogique et facile à apprendre et à utiliser ;

- langage interprété : utilisation interactive ou script exécuté ligne à ligne, pas de processus de compilation ;

- haut niveau : typage dynamique, gestion active de la mémoire, pour une plus grande facilité d'emploi ;

- multiparadigme : langage impératif et/ou orienté objet, selon les besoins et les capacités de chacun ;

- logiciel libre et ouvert, largement répandu (multiplateforme) et utilisé (forte communauté) ;

- riche bibliothèque standard : Batteries included ;

- riche bibliothèque externe : de nombreuses bibliothèques de qualité, dans divers domaines (y compris scientifiques), sont déjà disponibles.

L'objectif est bien d'apprendre un seul langage de haut niveau, permettant tout aussi bien des analyses rapides dans la vie de tous les jours - quelques lignes de code en interactif - que des programmes les plus complexes (projets de plus de 100 000 lignes).

Liens :

- Getting Started

- Python Advocacy

I-C - Index

- Index

- Page de recherche

II - Initiation à Python

II-A - Introduction

II-A-1 - Installation

Cette introduction repose essentiellement sur les outils suivants :

- Python 2.7 (inclus l'interpréteur de base et la bibliothèque standard) ;

- les bibliothèques scientifiques Numpy et Scipy ;

- la bibliothèque graphique Matplotlib ;

- un interpréteur évolué, p.ex. ipython ;

- un éditeur de texte évolué, p.ex. emacs, vi, gedit ou Atom.

Ces logiciels peuvent être installés indépendamment, de préférence sous Linux (p.ex. Ubuntu ou votre distribution préférée), ou sous Windows ou MacOS. Il existe également des distributions « clés en main » :

- Anaconda (multiplateforme) ;

- Python(x,y) (Windows) ;

- Enthought Canopy (Windows, MacOS, Linux, gratuite pour les étudiants du supérieur).

II-A-2 - Notions d'Unix

Les concepts suivants sont supposés connus :

- ligne de commande : exécutables et options ;

- arborescence : chemin relatif (./]...) et absolu (/...), navigation (cd) ;

- gestion des fichiers (ls, rm, mv) et répertoires (mkdir) ;

- gestion des exécutables : $PATH, chmod +x ;

- gestion des processus : &, Control-c, Control-z + bg ;

- variables d'environnement : export, .bashrc .

Liens :

- Quelques notions et commandes d'UNIX

- Introduction to Unix Study Guide

II-A-3 - L'invite de commande

Il existe principalement deux interpréteurs intéractifs de commandes Python :

- Python : interpréteur de base :

1. $ python

2. Python 2.7.6 (default, Mar 22 2014, 22:59:56)

3. [GCC 4.8.2] on linux2

4. Type "help", "copyright", "credits" or "license" for more information.

5. >>>

- Control-d pour sortir ;

- help(commande) pour obtenir l'aide d'une commande ;

- A priori, pas d'historique des commandes ni de complétion automatique.

L'interpréteur de base permet également d'interpréter un « script », c.-à-d. un ensemble de commandes regroupées dans un fichier texte (généralement avec une extension .py) : python mon_script.py.

- ipython : interpréteur évolué (avec historique et complétion automatique des commandes) :

1. $ ipython

2. Python 2.7.12 (default, Nov 19 2016, 06:48:10)

3. Type "copyright", "credits" or "license" for more information.

5. IPython 5.3.0 -- An enhanced Interactive Python.

6. ?         -> Introduction and overview of IPython's features.

7. %quickref -> Quick reference.

8. help      -> Python's own help system.

9. object?   -> Details about 'object', use 'object??' for extra details.

11. In [1]:

- Control-d pour sortir ;

- Tab pour la complétion automatique ;

- Haut et Bas pour le rappel des commandes ;

- Aide ipython : object? pour une aide sur un objet, object?? pour une aide plus complète (au niveau source) ;

- Commandes magic (voir %magic) :

- %run mon_script.py pour exécuter un script dans l'interpréteur,

- %debug pour lancer le mode débogage interactif postmortem,

- %cpaste pour coller et exécuter un code préformaté.

Liens :

- Tutorial

- IPython Tips & Tricks

II-B - Types de base

- None (rien)

- Chaînes de caractères : str

- Entre (simples ou triples) apostrophes ' ou guillemets " : 'Calvin', "Calvin'n'Hobbes", '''Deux\nlignes''', """'Pourquoi?' demanda-t-il."""

- Conversion : str(3.2)

- Types numériques :

- Booléens bool (vrai/faux) : True, False, bool(3)

- Entiers int (pas de valeur limite explicite, correspond au moins au long du C) : -2, int(2.1), int("4")

- Réels float (entre ±1.7e±308, correspond au double du C) : 2., 3.5e-6, float(3)

- Complexes complex : 1+2j, 5.1j, complex(-3.14), complex('j')

1. >>> 5 / 2       # !!! Division euclidienne par défaut dans Python 2.x !!!

2. 2

3. >>> float(5)/2  # ou ajouter 'from __future__ import division' en début de script

4. 2.5

5. >>> 6 // 2.5    # Division euclidienne explicite

6. 2.0

7. >>> 6 % 2.5     # Reste de la division euclidienne

8. 1.0

9. >>> (1 + 2j)**-0.5  # Puissance entière, réelle ou complexe

10. (0.5688644810057831-0.3515775842541429j)

- Objets itérables :

- Listes list : ['a', 3, [1, 2], 'a']

- Listes immuables tuple : (2, 3.1, [1+0j], 'a') (selon les conditions d'utilisation, les parenthèses ne sont pas toujours nécessaires)

- Listes à clés dict : {'a':1, 'b':[1, 2], 3:'c'}

- Ensembles non ordonnés d'éléments uniques set : {1, 2, 3, 2}

1. >>> l = ['a', True] # Définition d'une liste

2. >>> x, y = 1, 2.5   # Affectations multiples via tuples (les parenthèses ne sont pas nécessaires)

3. >>> range(5)        # Liste de 5 entiers commençant par 0

4. [0, 1, 2, 3, 4]

5. >>> l + [x, y]      # Concaténation de listes

6. ['a', True, 1, 2.5]

7. >>> {2, 1, 3} | {1, 2, 'a'}  # Union d'ensembles

8. {'a', 1, 2, 3}

- type(obj) retourne le type de l'objet, isinstance(obj, type) teste le type de l'objet.

1. >>> type(l)

2. <type 'list'>

3. >>> isinstance(l, tuple)

4. False

Liens :

- The Floating Point Guide

- What Every Computer Scientist Should Know About Floating-Point Arithmetic

II-C - Structures de programmation

- Les blocs sont définis par l'indentation (en général par pas de quatre espaces)(1).

Évitez autant que possible les caractères de tabulation, source de confusion. Configurez votre éditeur de texte pour qu'il n'utilise que des espaces.

- Une instruction par ligne en général (ou instructions séparées par ;).

- Les commentaires commencent par #, et s'étendent jusqu'à la fin de la ligne.

- Expression booléenne : une condition est une expression s'évaluant à True ou False :

- False : test logique faux (p.ex. 3 == 4), valeur nulle, chaîne vide (''), liste vide ([]), etc.

- True : test logique vrai (p.ex. 2 + 2 == 4), toute valeur ou objet non nul (et donc s'évaluant par défaut à True sauf exception)

- Tests logiques : ==, !=, >, >=, etc.

Ne pas confondre « = » (affectation d'une variable) et « == » (test logique d'égalité).

- Opérateurs logiques : and, or, not

>>> (5 >= 6) or (not 3 > 4)
True

- Opérateur ternaire (PEP 308) : value if condition else altvalue, p.ex.

>>> y = x**0.5 if (x > 0) else 0  # Retourne sqrt(max(x, 0))

- Expression conditionnelle : if condition: ... [elif condition2: ...] [else: ...], p.ex. :

1. if (i > 0):    # Condition principale

2. print "positif"

3. elif (i < 0):  # Condition secondaire (si nécessaire)

4. print "négatif"

5. else:          # Cas final (si nécessaire)

6. print "nul"

- Boucle for : for element in iterable , s'exécute sur chacun des éléments d'un objet itérable :

1. >>> for val in ['un', (2, 3), 4]:  # Itération sur une liste de 3 éléments

2. ...     print val

3. un

4. (2, 3)

5. 4

- continue : interrompt l'itération courante, et reprend la boucle à l'itération suivante,

- break : interrompt complètement la boucle.

la logique des boucles Python est assez différente des langages C[++]/fortran, pour lesquels l'itération porte sur les indices plutôt que sur les éléments eux-mêmes.

- Boucle while : while condition : se répète tant que la condition est vraie, ou jusqu'à une sortie explicite avec break.

Attention aux boucles infinies, dont la condition d'exécution reste invariablement vraie (typiquement un critère de convergence qui n'est jamais atteint). On peut toujours s'en protéger en testant en outre sur un nombre maximal (raisonnable) d'itérations :

1. niter = 0

2. while (error > 1e-6) and (niter < 100):

3. error = ...   # A priori, error va décroître, et la boucle s'interrompre

4. niter += 1    # Mais on n'est jamais assez prudent !

Exercices :

- Intégration : méthode des rectangles ★

- Fizz Buzz ★

- PGCD : algorithme d'Euclide ★★

II-D - Les chaînes de caractères

II-D-1 - Indexation

Les chaînes de caractères sont des objets itérables - c.-à-d. constitués d'éléments (ici les caractères) sur lesquels il est possible de « boucler » (p.ex. avec for) - et immuables - c.-à-d. dont les éléments individuels ne peuvent pas être modifiés intrinsèquement.

Comme en C[++], l'indexation en Python commence à 0 : le 1er élément d'une liste est l'élément n° 0, le 2e est le n° 1, etc. Les n éléments d'une liste sont donc indexés de 0 à n-1.

1. >>> alpha = 'abcdefghijklmnopqrstuvwxyz'

2. >>> len(alpha)

3. 26

4. >>> alpha[0]    # 1er élément (l'indexation commence à 0)

5. 'a'

6. >>> alpha[-1]   # = alpha[26-1=25], dernier élément (-2: avant-dernier, etc.)

7. 'z'

II-D-2 - Sous-liste (slice)

Des portions d'une chaîne peuvent être extraites en utilisant des slice (« tranches »), de notation générique [start=0]:[stop=len][:step=1]. P.ex.

1. >>> alpha[3:7]  # De l'élément n°3 (inclus) au n°7 (exclu), soit 7-3=4 éléments

2. 'defg'

3. >>> alpha[:3]   # Du n°0 (défaut) au n°3 (exclu), soit 3 éléments

4. 'abc'

5. >>> alpha[-3:]  # Du n°26-3=23 (inclus) au dernier inclus (défaut)

6. 'xyz'

7. >>> alpha[8:23:3] # Du n°8 (inclus) au n°23 (exclu), tous les 3 éléments

8. 'iloru'

9. >>> alpha[::5]    # Du 1er au dernier élément (défauts), tous les 5 éléments

10. 'afkpuz'

II-D-3 - Méthodes

Les chaînes de caractères disposent de nombreuses fonctionnalités - appelées « méthodes » en POO - facilitant leur manipulation :

1. >>> enfant, peluche = "Calvin", 'Hobbes'      # Affectations mutiples

2. >>> titre = enfant + ' et ' + peluche; titre  # +: Concaténation de chaînes

3. 'Calvin et Hobbes'

4. >>> titre.replace('et', '&')  # Remplacement de sous-chaînes (→ nouvelle chaîne)

5. 'Calvin & Hobbes'

6. >>> titre                     # titre est immuable et reste inchangé

7. 'Calvin et Hobbes'

8. >>> ' & '.join(titre.split(' et ')) # Découpage (split) et jonction (join)

9. 'Calvin & Hobbes'

10. >>> 'Hobbes' in titre               # in: Test d'inclusion

11. True

12. >>> titre.find("Hobbes")            # str.find: Recherche de sous-chaîne

13. 10

14. >>> titre.center(30, '-')

15. '-------Calvin et Hobbes-------'

16. >>> dir(str)                        # Liste toutes les méthodes des chaînes

II-D-4 - Formatage

Le système de formatage permet un contrôle précis de la conversion de variables en chaînes de caractères. Il s'appuie essentiellement sur la méthode str.format() :

1. >>> "{0} a {1} ans".format('Calvin', 6)  # args

2. 'Calvin a 6 ans'

3. >>> "{} a {} ans".format('Calvin', 6)    # Raccourci

4. 'Calvin a 6 ans'

5. >>> "{nom} a {age} ans".format(nom='Calvin', age=6)  # kwargs

6. 'Calvin a 6 ans'

7. >>> pi = 3.1415926535897931

8. >>> "{x:f}  {x:.2f}  {y:f}  {y:g}".format(x=pi, y=pi*1e9)  # Options de formatage

9. '3.141593  3.14  3141592653.589793  3.14159e+09'

print() affiche à l'écran (plus spécifiquement la sortie standard) la conversion d'une variable en chaîne de caractères :

1. >>> print "Calvin and Hobbes\nScientific progress goes 'boink'"

2. Calvin and Hobbes

3. Scientific progress goes 'boink'

4. >>> print "{0:2d} fois {1:2d} font {2}".format(3, 4, 3*4) # Formatage et affichage

5. 3 fois  4 font 12

Exercice :

- Tables de multiplication ★

II-E - Objets itérables

Les chaînes de caractères, listes, tuples et dictionnaires sont les objets itérables de base en Python. Les listes et dictionnaires sont modifiables (« mutables ») - leurs éléments constitutifs peuvent être changés à la volée - tandis que chaînes de caractères et les tuples sont immuables.

- Accès indexé : conforme à celui des chaînes de caractères

1. >>> l = range(1, 10, 2); l  # De 1 (inclus) à 10 (exclu) par pas de 2

2. [1, 3, 5, 7, 9]

3. >>> len(l)          # Nb d'éléments dans la liste (i varie de 0 à 4)

4. 5

5. >>> l[0], l[-2]     # 1er et avant-dernier élément (l'indexation commence à 0)

6. (1, 7)

7. >>> l[5]            # Erreur: indice hors-bornes

8. IndexError: list index out of range

9. >>> d = dict(a=1, b=2) # Création du dictionnaire {'a':1, 'b':2}

10. >>> d['a']          # Accès à une entrée via sa clé

11. 1

12. >>> d['c']          # Erreur: clé inexistante !

13. KeyError: 'c'

14. >>> d['c'] = 3; d   # Ajout d'une clé et sa valeur

15. {'a': 1, 'c': 3, 'b': 2}

16. >>> # Noter qu'un dictionnaire N'est PAS ordonné !

- Sous-listes (slices) :

1. >>> l[1:-1]        # Du 2e ('1') *inclus* au dernier ('-1') *exclu*

2. [3, 5, 7]

3. >>> l[1:-1:2]      # Idem, tous les 2 éléments

4. [3, 7]

5. >>> l[::2]         # Tous les 2 éléments (*start=0* et *stop=len* par défaut)

6. [1, 5, 9]

- Modification d'éléments d'une liste (chaînes et tuples sont immuables) :

1. >>> l[0] = 'a'; l            # Remplacement du 1er élément

2. ['a', 3, 5, 7, 9]

3. >>> l[1::2] = ['x', 'y']; l  # Remplacement d'éléments par *slices*

4. ['a', 'x', 5, 'y', 9]

5. >>> l + [1, 2]; l            # Concaténation (l reste inchangé)

6. ['a', 'x', 5, 'y', 9, 1, 2]

7. ['a', 'x', 5, 'y', 9]

8. >>> l += [1, 2]; l           # Concaténation sur place (l est modifié)

9. ['a', 'x', 5, 'y', 9, 1, 2]

10. >>> l.append('z'); l         # Ajout d'un élément en fin de liste

11. ['a', 'x', 5, 'y', 9, 1, 2, 'z']

12. >>> l.extend([-1, -2]); l    # Extension par une liste

13. ['a', 'x', 5, 'y', 9, 1, 2, 'z', -1, -2]

14. >>> del l[-6:]; l            # Efface les 6 derniers éléments de la liste

15. ['a', 'x', 5, 'y']

Attention à la modification des objets mutables :

1. >>> l = range(3)  # l pointe vers la liste [0, 1, 2]

2. >>> m = l; m      # m est un *alias* de la liste l: c'est le même objet

3. [0, 1, 2]

4. >>> id(l); id(m); m is l

5. 171573452         # id({obj}) retourne le n° d'identification en mémoire

6. 171573452         # m et l ont le même id:

7. True              # ils correspondent donc bien au même objet en mémoire

8. >>> l[0] = 'a'; m # puisque l a été modifiée, il en est de même de m

9. ['a', 1, 2]

10. >>> m = l[:]      # copie de tous les éléments de l dans une *nouvelle* liste m (clonage)

11. >>> id(l); id(m); m is l

12. 171573452

13. 171161228         # m a un id différent de l: il s'agit de 2 objets distincts

14. False             # (contenant éventuellement la même chose !)

15. >>> del l[-1]; m  # les éléments de m n'ont pas été modifiés

16. ['a', 1, 2]

- Liste en compréhension : elle permet la construction d'une liste à la volée

1. >>> [ i**2 for i in range(5) ] # Carré de tous les éléments de [0, ..., 4]

2. [0, 1, 4, 9, 16]

3. >>> [ 2*i for i in range(10) if (i%3 != 0) ]  # Compréhension conditionnelle

4. [2, 4, 8, 10, 14, 16]

5. >>> [ 10*i+j for i in range(3) for j in range(4) ]     # Double compréhension

6. [0, 1, 2, 3, 10, 11, 12, 13, 20, 21, 22, 23]

7. >>> [ [ 10*i+j for i in range(3) ] for j in range(4) ] # Compréhensions imbriquées

8. [[0, 10, 20], [1, 11, 21], [2, 12, 22], [3, 13, 23]]

9. >>> { i: i**2 for i in range(1, 5) }  # Dictionnaire en compréhension

10. {1: 1, 2: 4, 3: 9, 4: 16}

- Utilitaires sur les itérables :

1. >>> humans = ['Calvin', 'Wallace', 'Boule']

2. >>> for i in range(len(humans)):  # Boucle sur les indices de humans

3. ...     print i, humans[i]        # Accès explicite, pas pythonique :-(

4. 0 Calvin

5. 1 Wallace

6. 2 Boule

7. >>> for i, name in enumerate(humans):  # Boucle sur (indice, valeur) de humans

8. ...     print i, name                  # Pythonique :-D

9. 0 Calvin

10. 1 Wallace

11. 2 Boule

12. >>> animals = ['Hobbes', 'Gromit', 'Bill']

13. >>> for boy, dog in zip(humans, animals):  # Boucle simultanée sur 2 listes (ou +)

14. ...     print boy, 'et', dog

15. Calvin et Hobbes

16. Wallace et Gromit

17. Boule et Bill

18. >>> sorted(zip(humans, animals))  # Tri, ici sur le 1er élément de chaque tuple de la liste

19. [('Boule', 'Bill'), ('Calvin', 'Hobbes'), ('Wallace', 'Gromit')]

Exercices :

- Crible d'Ératosthène ★

- Carré magique ★★

II-F - Fonctions

Une fonction est un regroupement d'instructions impératives - assignations, branchements, boucles, etc. - s'appliquant sur des arguments d'entrée. C'est le concept central de la programmation impérative.

def permet de définir une fonction :

def fonction(arg1, arg2, ..., option1=valeur1, option2=valeur2, ...):

Les « args » sont des arguments nécessaires (c.-à-d. obligatoires), tandis que les « kwargs » - arguments de type option=valeur - sont optionnels, puisqu'ils possèdent une valeur par défaut. Si la fonction doit retourner une valeur, celle-ci est spécifiée par le mot-clé return.

Exemples :

1. def temp_f2c(tf):

2. """

3. Convertit une température en d° Fahrenheit `tf` en d° Celsius.

5. Exemple:

6. >>> temp_f2c(104)

7. 40.0

8. """

10. tc = (tf - 32.)/1.8       # Fahrenheit ? Celsius

12. return tc

Dans la définition d'une fonction, la première chaîne de caractères (appelé docstring) servira de documentation pour la fonction, accessible de l'interpréteur via p.ex. help(temp_f2c), ou temp_f2c? sous ipython. Elle se doit d'être tout à la fois pertinente, concise et complète. Elle peut également inclure des exemples d'utilisation (doctests, voir Développement piloté par les tests).

1. def mean_power(alist, power=1):

2. """

3. Retourne la racine `power` de la moyenne des éléments de `alist` à

4. la puissance `power`:

6. .. math:: \mu = (\frac{1}{N}\sum_{i=0}^{N-1} x_i^p)^{1/p}

8. `power=1` correspond à la moyenne arithmétique, `power=2` au *Root

9. Mean Squared*, etc.

11. Exemples:

12. >>> mean_power([1, 2, 3])

13. 2.0

14. >>> mean_power([1, 2, 3], power=2)

15. 2.160246899469287

16. """

18. s = 0.                  # Initialisation de la variable *s* comme *float*

19. for val in alist:       # Boucle sur les éléments de *alist*

20. s += val ** power   # *s* est augmenté de *val* puissance *power*

21. # *mean* = (somme valeurs / nb valeurs)**(1/power)

22. mean = (s / len(alist)) ** (1 / power)  # ATTENTION aux divisions euclidiennes !

24. return mean

Il faut noter plusieurs choses importantes :

- Python est un langage à typage dynamique, p.ex., le type des arguments d'une fonction n'est pas fixé a priori. Dans l'exemple précédent, alist peut être une list, un tuple ou tout autre itérable contenant des éléments pour lesquels les opérations effectuées - somme, exponentiation, division par un entier - ont été préalablement définies (p.ex. des entiers, des complexes, des matrices, etc.) ;

- le typage est fort, c.-à-d. que le type d'une variable ne peut pas changer à la volée. Ainsi, "abra" + "cadabra" a un sens (concaténation de chaînes), mais pas 3 + "cochons" (entier + chaîne) ;

- la définition d'une fonction se fait dans un « espace parallèle » où les variables ont une portée (scope) locale(2). Ainsi, la variable s définie dans la fonction mean_power n'interfère pas avec le « monde extérieur » ; inversement, la définition de mean_power ne connaît a priori rien d'autre que les variables explicitement définies dans la liste des arguments ou localement.

Exercice :

- Suite de Syracuse (fonction) ★

II-G - Bibliothèques et scripts

II-G-1 - Bibliothèques externes

Une bibliothèque (ou module) est un code fournissant des fonctionnalités supplémentaires - p.ex. des fonctions prédéfinies - à Python. Ainsi, le module math définit les fonctions et constantes mathématiques usuelles (sqrt(), pi, etc.).

Une bibliothèque est « importée » avec la commande import module. Les fonctionnalités supplémentaires sont alors accessibles dans l'espace de noms module via module.fonction :

1. >>> sqrt(2)                   # sqrt n'est pas une fonction standard de python

2. NameError: name 'sqrt' is not defined

3. >>> import math               # Importe tout le module 'math'

4. >>> dir(math)                 # Liste les fonctionnalités de 'math'

5. ['__doc__', '__name__', '__package__', 'acos', 'acosh', 'asin',

6. 'asinh', 'atan', 'atan2', 'atanh', 'ceil', 'copysign', 'cos', 'cosh',

7. 'degrees', 'e', 'exp', 'fabs', 'factorial', 'floor', 'fmod', 'frexp',

8. 'fsum', 'hypot', 'isinf', 'isnan', 'ldexp', 'log', 'log10', 'log1p',

9. 'modf', 'pi', 'pow', 'radians', 'sin', 'sinh', 'sqrt', 'tan', 'tanh',

10. 'trunc']

11. >>> math.sqrt(math.pi)        # Les fonctionnalités sont disponibles sous 'math'

12. 1.7724538509055159

13. >>> import math as M          # Importe 'math' dans l'espace 'M'

14. >>> M.sqrt(M.pi)

15. 1.7724538509055159

16. >>> from math import sqrt, pi # Importe uniquement 'sqrt' et 'pi' dans l'espace courant

17. >>> sqrt(pi)

18. 1.7724538509055159

Il est possible d'importer toutes les fonctionnalités d'une bibliothèque dans l'espace de noms courant :

1. >>> from math import *    # Argh ! Pas pythonique :-(

2. >>> sqrt(pi)

3. 1.7724538509055159

Cette pratique est cependant fortement déconseillée du fait des confusions dans les espaces de noms qu'elle peut entraîner :

>>> from cmath import *
>>> sqrt(-1)   # Quel sqrt: le réel ou le complexe ?

Nous verrons par la suite quelques exemples de modules de la Bibliothèque standard, ainsi que des Bibliothèques numériques de base orientées analyse numérique.

Exercice :

Flocon de Koch (programmation récursive) ★★★

II-G-2 - Bibliothèques personnelles et scripts

Vous pouvez définir vos propres bibliothèques en regroupant les fonctionnalités au sein d'un même fichier monfichier.py.

- Si ce fichier est importé (p.ex. import monfichier), il agira comme une bibliothèque.

- Si ce fichier est exécuté - p.ex. python ./monfichier.py - il agira comme un script.

Toutes les instructions d'un module qui ne sont pas encapsulées dans le __main__ (voir plus bas) sont interprétées et exécutées lors de l'import du module. Elles doivent donc en général se limiter à la définition de variables, de fonctions et de classes (en particulier, éviter les affichages ou les calculs longs).

Un code Python peut donc être :

- un module - s'il n'inclut que des définitions mais pas d'instruction exécutable en dehors d'un éventuel __main__ ;

- un exécutable - s'il inclut un __main__ ou des instructions exécutables ;

- ou les deux à la fois.

Exemple

Le code mean_power.py peut être importé comme une bibliothèque (p.ex. import mean_power) dans un autre code Python, ou bien être exécuté depuis la ligne de commande (p.ex. python mean_power.py), auquel cas la partie __main__ sera exécutée.

- #! (Hash-bang) : la première ligne d'un script défini l'interpréteur à utiliser(3) :

#!/usr/bin/env python

- Un fichier incluant des caractères non ASCII (p.ex. caractères accentués, ou symboles UTF tels que ±) doit définir le système d'encodage, généralement utf-8 :

# -*- coding: utf-8 -*-

Notez que les noms de variables, fonctions, etc. doivent être purement ASCII (a-zA-Z0-9_). De manière générale, favorisez la langue anglaise (variables, commentaires, affichages).

- """doc""" : la chaîne de documentation de la bibliothèque (docstring, PEP 257), qui sera utilisée comme aide en ligne du module (help(mean_power)), doit être la 1re instruction du script.

- from __future__ import division : permet de ne pas considérer les divisions entre entiers comme euclidiennes par défaut.

- if __name__ == '__main__' : permet de séparer le __main__ (c.-à-d. le corps du programme, à exécuter lors d'une utilisation en script) des définitions de fonctions et classes, permettant une utilisation en module.

II-H - Exceptions

Lorsqu'il rencontre une erreur dans l'exécution d'une instruction, l'interpréteur Python génère (raise) une erreur (Exception), de nature différente selon la nature de l'erreur : KeyError, ValueError, AttributeError, NameError, TypeError, IOError, NotImplementedError, etc. La levée d'une erreur n'est cependant pas nécessairement fatale, puisque Python dispose d'un mécanisme de gestion des erreurs.

Il est d'usage en Python d'utiliser la philosophie EAFP[4] : plutôt que de tester explicitement toutes les conditions de validité d'une instruction, on « tente sa chance » d'abord, quitte à gérer les erreurs a posteriori. Cette gestion des Exception se fait par la construction try ... except.

1. def lireEntier():

2. while True:

3. chaine = raw_input('Entrez un entier: ') # Lecture du clavier

4. try:

5. # La conversion en type entier génère `ValueError` si nécessaire

6. return int(chaine)

7. except ValueError:             # Gestion de l'exception ValueError

8. print "'{}' n'est pas un entier".format(chaine)

1. >>> lireEntier()

2. Entrez un entier: toto

3. 'toto' n'est pas un entier

4. Entrez un entier: 3,4

5. '3,4' n'est pas un entier

6. Entrez un entier: 4

7. 4

Dans l'élaboration d'un programme, gérez explicitement les erreurs que vous auriez pu tester a priori et pour lesquels il existe une solution de repli, et laissez passer les autres (ce qui provoquera éventuellement l'interruption du programme).

Évitez à tout prix les except nus, c.-à-d. ne spécifiant pas la ou les exceptions à gérer, car ils intercepteraient alors toutes les exceptions, y compris celles que vous n'aviez pas prévues ! Trouvez l'erreur dans le code suivant :

1. y = 2

2. try:

3. x = z                 # Copie y dans x

4. print "Tout va bien"

5. except:

6. print "Rien ne va plus"

Vos procédures doivent également générer des exceptions (documentées) - avec l'instruction raise Exception - si elles ne peuvent conclure leur action, à charge pour la procédure appelante de les gérer si besoin :

1. def diff_sqr(x, y):

2. """

3. Return x**2 - y**2 for x >= y, raise ValueError otherwise.

5. Exemples:

6. >>> diff_sqr(5, 3)

7. 16

8. >>> diff_sqr(3, 5)

9. Traceback (most recent call last):

10. ...

11. ValueError: x=3 < y=5

12. """

14. if x < y:

15. raise ValueError("x={} < y={}".format(x, y))

17. return x**2 - y**2

Avant de se lancer dans un calcul long et complexe, on peut vouloir tester la validité de certaines hypothèses fondamentales, soit par une structure if ... raise, ou plus facilement à l'aide d'assert (qui, si l'hypothèse n'est pas vérifiée, génère une AssertionError) :

1. def diff_sqr(x, y):

2. """

3. Returns x**2 - y**2 for x >= y, AssertionError otherwise.

4. """

6. assert x >= y, "x={} < y={}".format(x, y) # Test et msg d'erreur

7. return x**2 - y**2

La règle générale à retenir concernant la gestion des erreurs :

Fail early, fail often, fail better !

Exercice :

Jeu du plus ou moins (exceptions) ★

II-I - Classes

Un objet est une entité de programmation, disposant de ses propres états et fonctionnalités. C'est le concept central de la Programmation Orientée Objet.

Au concept d'objet sont liées les notions de :

- Classe : il s'agit d'un modèle d'objet, dans lequel sont définis ses propriétés usuelles. P.ex. la classe Forme peut représenter une forme plane caractérisée par sa couleur, et disposant de fonctionnalités propres, p.ex. change_couleur() ;

- Instanciation : c'est le fait générer un objet concret (une instance) à partir d'un modèle (une classe). P.ex. rosie = Forme('rose') crée une instance rosie à partir de la classe Forme et d'une couleur (chaîne de caractères 'rose') ;

- Attributs : variables internes décrivant l'état de l'objet. P.ex., rosie.couleur donne la couleur de la Forme rosie ;

- Méthodes : fonctions internes, s'appliquant en premier lieu sur l'objet lui-même (self), décrivant les capacités de l'objet. P.ex. rosie.change_couleur('bleu') change la couleur de la Forme rosie.

Toutes les méthodes d'une classe doivent au moins prendre self - représentant l'objet lui-même - comme premier argument.

- Surcharge d'opérateurs : cela permet de redéfinir les opérateurs et fonctions usuels (+, abs(), str(), etc.), pour simplifier l'écriture d'opérations sur les objets. Ainsi, on peut redéfinir les opérateurs de comparaison (<, >=, etc.) dans la classe Forme pour que les opérations du genre forme1 < forme2 aient un sens (p.ex. en comparant les aires).
Liste des méthodes standard et des surcharges d'opérateur ;

- Héritage de classe : il s'agit de définir une classe à partir d'une (ou plusieurs) classe(s) parente(s). La nouvelle classe hérite des attributs et méthodes de sa (ses) parente(s), que l'on peut alors modifier ou compléter. P.ex. la classe Rectangle hérite de la classe Forme (elle partage la notion de couleur et d'aire), et lui ajoute des méthodes propres à la notion de rectangle (p.ex. formule explicite de l'aire, étirement).

Toutes les classes doivent au moins hériter de la classe principale object.

Exemple de définition de classe

1. class Forme(object):  # *object* est la classe dont dérivent toutes les autres

3. """Une forme plane, avec éventuellement une couleur."""

5. def __init__(self, couleur=None):

6. """Initialisation d'une Forme, sans couleur par défaut."""

8. if couleur is None:

9. self.couleur = 'indéfinie'

10. else:

11. self.couleur = couleur

13. def __str__(self):

14. """

15. Surcharge de la fonction `str()`: l'affichage *informel* de

16. l'objet dans l'interpréteur, p.ex. `print a` sera résolu comme

17. `a.__str__()`

19. Retourne une chaîne de caractères.

20. """

22. return "Forme encore indéfinie de couleur {}".format(self.couleur)

24. def change_couleur(self, newcolor):

25. """Change la couleur de la Forme."""

27. self.couleur = newcolor

29. def aire(self):

30. """

31. Renvoi l'aire de la Forme.

33. L'aire ne peut pas être calculée dans le cas où la forme n'est

34. pas encore spécifiée: c'est ce que l'on appelle une méthode

35. 'abstraite', qui pourra être précisée dans les classes filles.

36. """

38. raise NotImplementedError(

39. "Impossible de calculer l'aire d'une forme indéfinie.")

41. def __cmp__(self, other):

42. """

43. Comparaison de deux Formes sur la base de leur aire.

45. Surcharge des opérateurs de comparaison de type `{self} <

46. {other}`: la comparaison sera résolue comme

47. `self.__cmp__(other)` et le résultat sera correctement

48. interprété.

50. .. WARNING:: cette construction n'est plus supportée en Python3.

51. """

53. return cmp(self.aire(), other.aire())  # Opérateur de comparaison

Exemple d'héritage de classe

1. class Rectangle(Forme):

3. """

4. Un Rectangle est une Forme particulière.

6. La classe fille hérite des attributs et méthodes de la

7. classe mère, mais peut les surcharger (i.e. en changer la

8. définition), ou en ajouter de nouveaux :

10. - les méthodes `Rectangle.change_couleur()` et

11. `Rectangle.__cmp__()` dérivent directement de

12. `Forme.change_couleur()` et `Forme.__cmp__()`;

13. - `Rectangle.__str__()` surcharge `Forme.__str__()`;

14. - `Rectangle.aire()` définit la méthode jusqu'alors abstraite

15. `Forme.aire()`;

16. - `Rectangle.allonger()` est une nouvelle méthode propre à

17. `Rectangle`.

18. """

20. def __init__(self, longueur, largeur, couleur=None):

21. """

22. Initialisation d'un Rectangle longueur × largeur, sans couleur par

23. défaut.

24. """

26. # Initialisation de la classe parente (nécessaire pour assurer

27. # l'héritage)

28. Forme.__init__(self, couleur)

30. # Attributs propres à la classe Rectangle

31. self.longueur = longueur

32. self.largeur = largeur

34. def __str__(self):

35. """Surcharge de `Forme.__str__()`."""

37. return "Rectangle {}x{}, de couleur {}".format(

38. self.longueur, self.largeur, self.couleur)

40. def aire(self):

41. """

42. Renvoie l'aire du Rectangle.

44. Cette méthode définit la méthode abstraite `Forme.area()`,

45. pour les Rectangles uniquement.

46. """

48. return self.longueur * self.largeur

50. def allonger(self, facteur):

51. """Multiplie la *longueur* du Rectangle par un facteur"""

53. self.longueur *= facteur

Il est traditionnel de commencer les noms de classes avec une majuscule (Forme), et les noms d'instances de classes (les variables) avec une minuscule (rosie).

Exemples

- Formes (POO)

- Cercle circonscrit (POO, argparse)

Études de cas

- turtle.Vec2D

- fractions.Fraction

Exercices :

- Animaux (POO/TDD) ★

- Jeu de la vie (POO) ★★

II-J - Entrées-sorties

II-J-1 - Intéractif

Comme nous avons pu le voir précédemment, l'affichage à l'écran se fait par print, la lecture du clavier par raw_input.

II-J-2 - Fichiers

La gestion des fichiers (lecture et écriture) se fait à partir de la fonction open() retournant un objet de type file :

1. # ========== ÉCRITURE ==========

2. outfile = open("carres.dat", 'w') # Ouverture du fichier "carres.dat" en écriture

3. for i in range(1, 10):

4. outfile.write("{}  {}\n".format(i, i**2)) # Noter la présence du '\n' (non automatique)

5. outfile.close()                   # Fermeture du fichier (nécessaire)

7. # ========== LECTURE ==========

8. infile = open("carres.dat")  # Ouverture du fichier "carres.dat" en lecture

9. for line in infile:          # Boucle sur les lignes du fichier

10. if line.strip().startswith('#'): # Ne pas considérer les lignes "commentées"

11. continue

12. try:                     # Essayons de lire 2 entiers sur cette ligne

13. x, x2 = [ int(tok) for tok in line.split() ]

14. except ValueError:       # Gestion des erreurs

15. print "Cannot decipher line '{}'".format(line)

16. continue

17. print "{}**3 = {}".format(x, x**3)

II-K - Éléments passés sous silence

Cette (brève) introduction à Python se limite à des fonctionnalités relativement simples du langage. De nombreuses fonctionnalités du langage n'ont pas été abordées :

- v ;

- arguments anonymes : *args et **kwargs ;

- fonction anonyme : lambda x, y: x + y ;

- itérateurs et générateurs : yield ;

- gestion de contexte : with ... as (PEP 243) ;

- décorateurs : fonction sur une fonction ou une classe (@staticmethod, etc.) ;

- héritages multiples et méthodes de résolution ;

- etc.

Ces fonctionnalités peuvent évidemment être très utiles, mais ne sont généralement pas strictement indispensables pour une première utilisation de Python dans un contexte scientifique.

II-K-1 - Python 3.x

Pour des raisons historiques autant que pratiques(4), ce cours présente le langage Python dans sa version 2.x. Pourtant, le développement actuel de Python se fait uniquement sur la branche 3.x, qui constitue une remise à plat non rétrocompatible du langage, et la branche 2.x ne sera a priori plus supporté au-delà de 2020 (PEP 466). Il est donc préférable, si vous vous lancez dans un développement substantiel, de passer aussi rapidement que possible à Python 3 (voir Python 2 vs. python 3).

Python 3 apporte quelques changements fondamentaux, notamment :

- print n'est plus un mot-clé mais une fonction : print(...) ;

- l'opérateur / ne réalise plus la division euclidienne entre les entiers, mais toujours la division réelle ;

- un support complet (mais encore complexe) des chaînes Unicode ;

- un nouveau système de formatage des chaînes de caractères (f-string du PEP 498 à partir de Python 3.6).

III - Bibliothèque standard

Python dispose d'une très riche bibliothèque de modules étendant les capacités du langage dans de nombreux domaines : nouveaux types de données, interactions avec le système, gestion des fichiers et des processus, protocoles de communication (Internet, mail, FTP, etc.), multimédia, etc.

- The Python Standard Library (v2.7)

- Python Module of the Week

III-A - Gestion des arguments/options de la ligne de commande

Utilisation de sys.argv

Le module sys permet un accès direct aux arguments de la ligne de commande, via la liste sys.argv : sys.argv[0] contient le nom du script exécuté, sys.argv[1] le nom du 1er argument (s'il existe), etc. P.ex. :

1. # Gestion simplifiée d'un argument entier sur la ligne de commande

2. import sys

4. if sys.argv[1:]:  # Présence d'au moins un argument sur la ligne de commande

5. try:

6. n = int(sys.argv[1])  # Essayer de lire le 1er argument comme un entier

7. except ValueError:

8. raise ValueError("L'argument '{}' n'est pas un entier"

9. .format(sys.argv[1]))

10. else:                        # Pas d'argument sur la ligne de commande

11. n = 101                  # Valeur par défaut

Module argparse

Pour une gestion avancée des arguments et/ou options de la ligne de commande, il est préférable d'utiliser le module argparse. P.ex. :

1. import argparse

3. parser = argparse.ArgumentParser(

4. usage="%(prog)s [-p/--plot] [-i/--input coordfile | x1,y1 x2,y2 x3,y3]",

5. description=__doc__)

6. parser.add_argument('coords', nargs='*', type=str, metavar='x,y',

7. help="Coordinates of point")

8. parser.add_argument('-i', '--input', nargs='?', type=file,

9. help="Coordinate file (one 'x,y' per line)")

10. parser.add_argument('-p', '--plot', action="store_true", default=False,

11. help="Draw the circumscribed circle")

12. parser.add_argument('--version', action='version', version=__version__)

14. args = parser.parse_args()

Cette solution génère automatiquement une aide en ligne, p.ex. :

1. $ python circonscrit.py -h

2. usage: circonscrit.py [-p/--plot] [-i/--input coordfile | x1,y1 x2,y2 x3,y3]

4. Compute the circumscribed circle to 3 points in the plan.

6. positional arguments:

7. x,y                   Coordinates 'x,y' of point

9. optional arguments:

10. -h, --help            show this help message and exit

11. -i [INPUT], --input [INPUT]

12. Coordinate file (one 'x,y' per line)

13. -p, --plot            Draw the circumscribed circle

14. --version             show program's version number and exit

III-B - [c]Pickle : sérialisation des données

Les modules pickle/cPickle permettent la sauvegarde pérenne d'objets python (« sérialisation »).

1. >>> d = dict(a=1, b=2, c=3)

2. >>> l = range(10000)

3. >>> import cPickle as pickle         # 'cPickle' est + rapide que 'pickle'

4. >>> pkl = open('archive.pkl', 'w')   # Overture du fichier en écriture

5. >>> pickle.dump((d, l), pkl, protocol=-1)  # Sérialisation du tuple (d, l)

6. >>> pkl.close()                      # *IMPORTANT !* Fermeture du fichier

7. >>> d2, l2 = pickle.load(open('archive.pkl')) # Désérialisation (relecture)

8. >>> (d == d2) and (l == l2)

9. True

les pickles ne sont pas un format d'échange de données. Ils sont spécifiques à python, et peuvent dépendre de la machine utilisée.

III-C - Batteries included

Quelques modules de la bibliothèque standard qui peuvent être d'intérêt :

- math : accès aux fonctions mathématiques réelles

>>> math.asin(math.sqrt(2) / 2) / math.pi * 180
45.00000000000001

- cmath : accès aux fonctions mathématiques complexes

>>> cmath.exp(cmath.pi * 1j) + 1
1.2246467991473532e-16j

- a ;

- collections définit de nouveaux types spécialisés, p.ex. collections.OrderedDict, un dictionnaire ordonné ;

- itertools fournit des générateurs de boucle (itérateurs) et de combinatoire :

>>> [ ''.join(item) for item in itertools.combinations('ABCD', 2) ]
['AB', 'AC', 'AD', 'BC', 'BD', 'CD']

- interactions avec le système :

- sys, os : interface système,

- shutil : opérations sur les fichiers (copy, move, etc.),

- subprocess : exécution de commandes système,

- glob : métacaractères du shell (p.ex. toto?.*) ;

- e : re ;

- gestion du temps (time) et des dates (datetime, calendar) ;

- fichiers compressés et archives : gzip, bz2, zipfile, tarfile ;

- Lecture et sauvegarde des données (outre pickle/cPickle) :

- csv : lecture/sauvegarde de fichiers CSV,

- ConfigParser : fichiers de configuration,

- json : lightweight data interchange format ;

- lecture d'une URL (p.ex. page web) : urllib2.

III-D - Text/Graphical User Interfaces

- TUI : curses

- GUI : Tkinter,

bibliothèques externes :

- TUI : termcolor (texte coloré ANSII), https://pypi.python.org/pypi/blessed (mise en page)

- GUI : PyGI (GTK3), PyQt / pySide (Qt), wxPython (wxWidgets)

IV - Bibliothèques numériques de base

IV-A - Numpy

numpy est une bibliothèque numérique apportant le support efficace de larges tableaux multidimensionnels, et de routines mathématiques de haut niveau (fonctions spéciales, algèbre linéaire, statistiques, etc.).

La convention d'import utilisé dans les exemples est « import numpy as N ».

Liens :

- Numpy User Guide

- Numpy Reference

IV-A-1 - Tableaux

Un numpy.ndarray (généralement appelé array) est un tableau multidimensionnel homogène : tous les éléments doivent avoir le même type, en général numérique. Les différentes dimensions sont appelées des axes, tandis que le nombre de dimensions - 0 pour un scalaire, 1 pour un vecteur, 2 pour une matrice, etc. - est appelé le rang.

1. >>> import numpy as N       # Import de la bibliothèque numpy avec le surnom N

2. >>> a = N.array([1, 2, 3])  # Création d'un array 1D à partir d'une liste d'entiers

3. >>> a.ndim      # Rang du tableau

4. 1               # Vecteur (1D)

5. >>> a.shape     # Format du tableau: par définition, len(shape)=ndim

6. (3,)            # Vecteur 1D de longueur 3

7. >>> a.dtype     # Type des données du tableau

8. dtype('int32')  # Python 'int' = numpy 'int32' = C 'long'

9. >>> # Création d'un tableau 2D de float (de 0. à 12.) de shape 4×3

10. >>> b = N.arange(12, dtype=float).reshape(4, 3); b

11. array([[  0.,   1.,   2.],

12. [  3.,   4.,   5.],

13. [  6.,   7.,   8.],

14. [  9.,  10.,  11.]])

15. >>> b.shape              # Nb d'éléments le long de chacune des dimensions

16. (4, 3)                   # 4 lignes, 3 colonnes

17. >>> b.size               # Nb *total* d'éléments dans le tableau

18. 12                       # Par définition, size = prod(shape)

19. >>> b.dtype

20. dtype('float64')         # Python 'float' = numpy 'float64' = C 'double'

IV-A-1-a - Création de tableaux

- numpy.array() : convertit une liste d'éléments homogènes ou coercibles

1. >>> N.array([[1, 2],[3., 4.]]) # Liste de listes d'entiers et de réels

2. array([[ 1.,  2.],             # Tableau 2D de réels

3. [ 3.,  4.]])

- numpy.zeros() (resp. numpy.ones() et numpy.full()) : crée un tableau de format donné rempli de zéros (resp. de uns et d'une valeur fixe)

1. >>> N.zeros((2, 1))  # Shape (2, 1): 2 lignes, 1 colonne, float par défaut

2. array([[ 0.],

3. [ 0.]])

4. >>> N.ones((1, 2), dtype=bool) # Shape (1, 2): 1 ligne, 2 colonnes, type booléen

5. array([[True, True]], dtype=bool)

6. >>> N.full((2, 2), N.NaN)

7. array([[ nan,  nan],

8. [ nan,  nan]])

- numpy.arange() : crée une séquence de nombres, en spécifiant éventuellement le start, le end et le step (similaire à range() pour les listes)

1. >>> N.arange(10, 30, 5) # De 10 à 30 (exclu) par pas de 5, type entier par défaut

2. array([10, 15, 20, 25])

3. >>> N.arange(0.5, 2.1, 0.3) # Accepte des réels en argument, DANGER !

4. array([ 0.5,  0.8,  1.1,  1.4,  1.7,  2. ])

- numpy.linspace() : répartition uniforme d'un nombre fixe de points entre un start et un end (préférable à numpy.arange() sur des réels).

>>> N.linspace(0, 2*N.pi, 5) # 5 nb entre 0 et 2p *inclus*, type réel par défaut
array([ 0.,  1.57079633,  3.14159265,  4.71238898,  6.28318531])

- numpy.meshgrid() est similaire à numpy.linspace() en 2D ou plus :

1. >>> # 5 points entre 0 et 2 en "x", et 3 entre 0 et 1 en "y"

2. >>> x = N.linspace(0, 2, 5); x           # Tableau 1D des x, (5,)

3. array([ 0. ,  0.5,  1. ,  1.5,  2. ])

4. >>> y = N.linspace(0, 1, 3); y           # Tableau 1D des y, (3,)

5. array([ 0. ,  0.5,  1. ])

6. >>> xx, yy = N.meshgrid(x, y)            # Tableaux 2D des x et des y

7. >>> xx                                   # Tableau 2D des x, (3, 5)

8. array([[ 0. ,  0.5,  1. ,  1.5,  2. ],

9. [ 0. ,  0.5,  1. ,  1.5,  2. ],

10. [ 0. ,  0.5,  1. ,  1.5,  2. ]])

11. >>> yy                                   # Tableau 2D des y, (3, 5)

12. array([[ 0. ,  0. ,  0. ,  0. ,  0. ],

13. [ 0.5,  0.5,  0.5,  0.5,  0.5],

14. [ 1. ,  1. ,  1. ,  1. ,  1. ]])

- numpy.mgrid permet de générer des rampes d'indices (entiers) ou de coordonnées (réels) de rang arbitraire avec une notation évoluée faisant appel aux Index tricks. Équivalent à numpy.linspace() en 1D et similaire (mais différent) à numpy.meshgrid() en 2D.

1. >>> N.mgrid[0:4, 1:6:2]   # Grille 2D d'indices (entiers)

2. array([[[0, 0, 0],        # 0:4 = [0, 1, 2, 3] selon l'axe 0

3. [1, 1, 1],

4. [2, 2, 2],

5. [3, 3, 3]],

6. [[1, 3, 5],        # 1:6:2 = [1, 3, 5] selon l'axe 1

7. [1, 3, 5],

8. [1, 3, 5],

9. [1, 3, 5]]])

10. >>> N.mgrid[0:2*N.pi:5j]  # Rampe de coordonnées (réels): 5 nb de 0 à 2p (inclus)

11. array([ 0.,  1.57079633,  3.14159265,  4.71238898,  6.28318531])

12. >>> # 3 points entre 0 et 1 selon l'axe 0, et 5 entre 0 et 2 selon l'axe 1

13. >>> z = N.mgrid[0:1:3j, 0:2:5j]; z

14. array([[[ 0. ,  0. ,  0. ,  0. ,  0. ],  # Axe 0 variable, axe 1 constant

15. [ 0.5,  0.5,  0.5,  0.5,  0.5],

16. [ 1. ,  1. ,  1. ,  1. ,  1. ]],

17. [[ 0. ,  0.5,  1. ,  1.5,  2. ],  # Axe 0 constant, axe 1 variable

18. [ 0. ,  0.5,  1. ,  1.5,  2. ],

19. [ 0. ,  0.5,  1. ,  1.5,  2. ]]])

20. >>> z.shape

21. (2, 3, 5)    # 2 plans 2D (x, y) de 3 lignes (y) × 5 colonnes (x)

22. >>> N.mgrid[0:1:5j, 0:2:7j, 0:3:9j].shape

23. (3, 5, 7, 9) # 3 volumes 3D (x, y, z) de 5 plans (z) × 7 lignes (y) × 9 colonnes (x)

Attention à l'ordre de variation des indices dans les tableaux multidimensionnels, et aux différences entre numpy.meshgrid() et numpy.mgrid.

- numpy.random.rand() crée un tableau d'un format donné de réels aléatoires dans [0, 1[; numpy.random.randn() génère un tableau d'un format donné de réels tirés aléatoirement d'une distribution gaussienne (normale) standard .
Le sous-module numpy.random fournit des générateurs de nombres aléatoires pour de nombreuses distributions discrètes et continues.

IV-A-1-b - Manipulations sur les tableaux

Les array 1D sont indexables comme les listes standard. En dimension supérieure, chaque axe est indexable indépendamment.

1. >>> x = N.arange(10);   # Rampe 1D

2. >>> x[1::3] *= -1; x    # Modification sur place ("in place")

3. array([ 0, -1,  2,  3, -4,  5,  6, -7,  8,  9])

Slicing

Les sous-tableaux de rang < N d'un tableau de rang N sont appelées slices : le (ou les) axe(s) selon le(s)quel(s) la slice a été découpée, devenu(s) de longueur 1, est (sont) éliminé(s).

1. >>> y = N.arange(2*3*4).reshape(2, 3, 4); y  # 2 plans, 3 lignes, 4 colonnes

2. array([[[ 0,  1,  2,  3],

3. [ 4,  5,  6,  7],

4. [ 8,  9, 10, 11]],

5. [[12, 13, 14, 15],

6. [16, 17, 18, 19],

7. [20, 21, 22, 23]]])

8. >>> y[0, 1, 2]  # 1er plan (axe 0), 2e ligne (axe 1), 3e colonne (axe 2)

9. 6               # Scalaire, shape *()*, ndim 0

10. >>> y[0, 1]     # = y[0, 1, :] 1er plan (axe 0), 2e ligne (axe 1)

11. array([4, 5, 6, 7])        # Shape (4,)

12. >>> y[0]        # = y[0, :, :] 1er plan (axe 0)

13. array([[ 0,  1,  2,  3],

14. [ 4,  5,  6,  7],

15. [ 8,  9, 10, 11]])  # Shape (3, 4)

16. >>> y[0][1][2]  # = y[0, 1, 2] en ~4× plus lent (slices successives)

17. 6

18. >>> y[:, -1]    # = y[:, 2, :] Dernière slice selon le 2e axe

19. array([[ 8,  9, 10, 11],

20. [20, 21, 22, 23]])  # Shape (2, 4)

21. >>> y[..., 0]   # = y[:, :, 0] 1re slice selon le dernier axe

22. array([[ 0,  4,  8],

23. [12, 16, 20]])      # Shape (2, 3)

24. >>> # On peut vouloir garder explicitement la dimension "tranchée"

25. >>> y[..., 0:1] # 1re slice selon le dernier axe *en gardant le rang originel*

26. array([[[ 0],

27. [ 4],

28. [ 8]],

29. [[12],

30. [16],

31. [20]]])

32. >>> y[..., 0:1].shape

33. (2, 3, 1)  # Le dernier axe a été conservé, il ne contient pourtant qu'un seul élément

Modification de format

numpy.ndarray.reshape() modifie le format d'un tableau sans modifier le nombre total d'éléments :

1. >>> y = N.arange(6).reshape(2, 3); y # Shape (6,) ? (2, 3) (*size* inchangé)

2. array([[0, 1, 2],

3. [3, 4, 5]])

4. >>> y.reshape(2, 4)        # Format incompatible (*size* serait modifié)

5. ValueError: total size of new array must be unchanged

numpy.ndarray.ravel() « déroule » tous les axes et retourne un tableau de rang 1 :

1. >>> y.ravel()              # *1st axis slowest, last axis fastest*

2. array([ 0, 1, 2, 3, 4, 5]) # Shape (2, 3) ? (6,) (*size* inchangé)

3. >>> y.ravel('F')           # *1st axis fastest, last axis slowest* (ordre Fortran)

4. array([0, 3, 1, 4, 2, 5])

numpy.ndarray.transpose() transpose deux axes, par défaut les derniers (raccourci : numpy.ndarray.T) :

1. >>> y.T         # Transposition = y.transpose() (voir aussi *rollaxis*)

2. array([[0, 3],

3. [1, 4],

4. [2, 5]])

numpy.ndarray.squeeze() élimine tous les axes de dimension 1. numpy.expand_dims() ajoute un axe de dimension 1 en position arbitraire. Cela est également possible en utilisant notation slice avec numpy.newaxis.

1. >>> y[..., 0:1].squeeze()               # Élimine *tous* les axes de dimension 1

2. array([0, 3])

3. >>> N.expand_dims(y[..., 0], -1).shape  # Ajoute un axe de dim. 1 en dernière position

4. (2, 1)

5. >>> y[:, N.newaxis].shape               # Ajoute un axe de dim. 1 en 2de position

6. (2, 1, 3)

numpy.resize() modifie le format en modifiant le nombre total d'éléments :

1. >>> N.resize(N.arange(4), (2, 4)) # Complétion avec des copies du tableau

2. array([[0, 1, 2, 3],

3. [0, 1, 2, 3]])

4. >>> N.resize(N.arange(4), (4, 2))

5. array([[0, 1],

6. [2, 3],

7. [0, 1],

8. [2, 3]])

Exercice :

- Inversion de matrice ★

Stacking

1. >>> a = N.arange(5); a

2. array([0, 1, 2, 3, 4])

3. >>> N.hstack((a, a))  # Stack horizontal (le long des colonnes)

4. array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4])

5. >>> N.vstack((a, a))  # Stack vertical (le long des lignes)

6. array([[0, 1, 2, 3, 4],

7. [0, 1, 2, 3, 4]])

8. >>> N.dstack((a, a))  # Stack en profondeur (le long des plans)

9. array([[[0, 0],

10. [1, 1],

11. [2, 2],

12. [3, 3],

13. [4, 4]]])

Broadcasting

L'array broadcasting définit les règles selon lesquelles deux tableaux de formats différents peuvent éventuellement s'apparier :

1. Deux tableaux sont compatibles (broadcastable) si, pour chaque axe, soit les tailles sont égales, soit l'une d'elles est exactement égale à 1. P.ex. (5, 3) et (1, 3) sont des formats broadcastable, (5, 3) et (5, 1) également, mais (5, 3) et (3, 1) ne le sont pas ;

2. Si un tableau a un axe de taille 1, le tableau sera dupliqué à la volée autant de fois que nécessaire selon cet axe pour atteindre la taille de l'autre tableau le long de cet axe. P.ex. un tableau (2, 1, 3) pourra être transformé en tableau (2, 5, 3) en le dupliquant 5 fois le long du 2e axe (axis=1) ;

3. La taille selon chaque axe après broadcast est égale au maximum de toutes les tailles d'entrée le long de cet axe. P.ex. (5, 3, 1) × (1, 3, 4) → (5, 3, 4) ;

4. Si un des tableaux a un rang (ndim) inférieur à l'autre, alors son format (shape) est précédé d'autant de 1 que nécessaire pour atteindre le même rang. P.ex. (5, 3, 1) × (4,) = (5, 3, 1) × (1, 1, 4) → (5, 3, 4).

1. >>> a = N.arange(6).reshape(2, 3); a # Shape (2, 3)

2. array([[0, 1, 2],

3. [3, 4, 5]])

4. >>> b = N.array([10, 20, 30]); b     # Shape (3,)

5. array([10, 20, 30])

6. >>> a + b                     # Shape (3,) ~ (1, 3) ? (2, 3) = (1, 3) copié 2 fois

7. array([[10, 21, 32],

8. [13, 24, 35]])         # Shape (2, 3)

9. >>> c = N.array([10, 20]); c  # Shape (2,)

10. array([10, 20])

11. >>> a + c                     # Shape (2,) ~ (1, 2) incompatible avec (2, 3)

12. ValueError: shape mismatch: objects cannot be broadcast to a single shape

13. >>> c[:, N.newaxis]           # = c.reshape(-1, 1) Shape (2, 1)

14. array([[10],

15. [20]])

16. >>> a + c[:, N.newaxis]       # Shape (2, 1) ? (2, 3) = (2, 1) copié 3 fois

17. array([[10, 11, 12],

18. [23, 24, 25]])

Voir également cette présentation.

Indexation évoluée

1. >>> a = N.linspace(-1, 1, 5); a

2. array([-1. , -0.5,  0. ,  0.5,  1. ])

3. >>> a >= 0                # Test logique : tableau de booléens

4. array([False, False,  True,  True,  True], dtype=bool)

5. >>> (a >= 0).nonzero()    # Indices des éléments ne s'évaluant pas à False

6. (array([2, 3, 4]),)       # Indices des éléments >= 0

7. >>> a[(a >= 0).nonzero()] # Indexation par un tableau d'indices, pas pythonique :-(

8. array([ 0. ,  0.5,  1. ])

9. >>> a[a >= 0]             # Indexation par un tableau de booléens, pythonique :-D

10. array([ 0. ,  0.5,  1. ])

11. >>> a[a < 0] -= 10; a     # = N.where(a < 0, a - 10, a)

12. array([-11. , -10.5,   0. ,   0.5,   1. ])

IV-A-1-c - Opérations de base

1. >>> a = N.arange(3); a # Shape (3,), type *int*

2. array([0, 1, 2])

3. >>> b = 1.             # ~ Shape (), type *float*

4. >>> c = a + b; c       # *Broadcasting*: () ? (1,) ? (3,)

5. array([ 1.,  2.,  3.]) # *Upcasting*: int ? float

6. >>> a += 1; a          # Modification *in place* (plus efficace si possible)

7. array([ 1.,  2.,  3.])

8. >>> a.mean()           # *ndarray* dispose de nombreuses méthodes numériques de base

9. 2.0

Opérations sur les axes

1. >>> x = N.random.permutation(6).reshape(3, 2); x # 3 lignes, 2 colonnes

2. array([[3, 4],

3. [5, 1],

4. [0, 2]])

5. >>> x.min()          # Minimum global (comportement par défaut : `axis=None`)

6. 0

7. >>> x.min(axis=0)    # Minima le long de l'axe 0 (i.e. l'axe des lignes)

8. array([0, 1])        # ce sont les minima colonne par colonne : (*3*, 2) ? (2,)

9. >>> x.min(axis=1)    # Minima le long de l'axe 1 (i.e. l'axe des colonnes)

10. array([3, 1, 0])     # ce sont les minima ligne par ligne (3, *2*) ? (3,)

11. >>> x.min(axis=1, keepdims=True) # Idem mais en *conservant* le format originel

12. array([[3],

13. [1],

14. [0]])

15. >>> x.min(axis=(0, 1)) # Minima le long des axes 0 *et* 1 (c.-à-d. ici tous les axes)

16. 0

Opérations matricielles

Les opérations de base s'appliquent sur les éléments des tableaux, et n'ont pas une signification matricielle par défaut :

1. >>> m = N.arange(4).reshape(2, 2); m # Tableau de rang 2

2. array([[0, 1],

3. [2, 3]])

4. >>> i = N.identity(2, dtype=int); i  # Tableau "identité" de rang 2 (type entier)

5. array([[1, 0],

6. [0, 1]])

7. >>> m * i             # Attention ! opération * sur les éléments

8. array([[0, 0],

9. [0, 3]])

10. >>> N.dot(m, i)       # Multiplication *matricielle* des tableaux : M × I = M

11. array([[0, 1],

12. [2, 3]])

Il est possible d'utiliser systématiquement les opérations matricielles en manipulant des numpy.matrix plutôt que de numpy.ndarray :

1. >>> N.matrix(m) * N.matrix(i)  # Opération * entre matrices

2. matrix([[0, 1],

3. [2, 3]])

Le sous-module numpy.linalg fournit des outils spécifiques au calcul matriciel (inverse, déterminant, valeurs propres, etc.).

Ufuncs

numpy fournit de nombreuses fonctions mathématiques de base (numpy.exp(), numpy.atan2(), etc.) s'appliquant directement sur les éléments des tableaux d'entrée :

1. >>> x = N.linspace(0, 2*N.pi, 5)  # [0, p/2, p, 3p/2, 2p]

2. >>> y = N.sin(x); y               # sin(x) = [0, 1, 0, -1, 0]

3. array([  0.00000000e+00,   1.00000000e+00,   1.22460635e-16,

4. -1.00000000e+00,  -2.44921271e-16])

5. >>> y == [0, 1, 0, -1, 0]         # Test d'égalité stricte (élément par élément)

6. array([ True,  True, False,  True, False], dtype=bool) # Attention aux calculs en réels (float) !

7. >>> N.all(N.sin(x) == [0, 1, 0, -1, 0]) # Test d'égalité stricte de tous les éléments

8. False

9. >>> N.allclose(y, [0, 1, 0, -1, 0])     # Test d'égalité numérique de tous les éléments

10. True

Exercices :

- Median Absolute Deviation ★

- Distribution du pull ★★★

IV-A-2 - Tableaux évolués

Types composés

Par définition, tous les éléments d'un tableau homogène doivent être du même type. Cependant, outre les types scalaires élémentaires - bool, int, float, complex, str, etc. - numpy supporte les types composés, c.-à-d. incluant plusieurs sous-éléments de types différents :

1. >>> dt = N.dtype([('nom', 'S10'),    # 1er  élément : une chaîne de 10 caractères

2. ...               ('age', 'i'),      # 2e élément : un entier

3. ...               ('taille', 'd')])  # 3e élément : un réel (double)

4. >>> arr = N.array([('Calvin', 6, 1.20), ('Hobbes', 5, 1.80)], dtype=dt); arr

5. array([('Calvin', 6, 1.2), ('Hobbes', 6, 1.8)],

6. dtype=[('nom', '|S10'), ('age', '<i4'), ('taille', '<f8')])

7. >>> arr[0]                           # Accès par élément

8. ('Calvin', 6, 1.2)

9. >>> arr['nom']                       # Accès par sous-type

10. array(['Calvin', 'Hobbes'], dtype='|S10')

11. >>> rec = arr.view(N.recarray); arr  # Vue de type 'recarray'

12. rec.array([('Calvin', 6, 1.2), ('Hobbes', 5, 1.8)],

13. dtype=[('nom', '|S10'), ('age', '<i4'), ('taille', '<f8')])

14. >>> rec.nom                          # Accès direct par attribut

15. chararray(['Calvin', 'Hobbes'], dtype='|S10')

Les tableaux structurés sont très puissants pour manipuler des données (semi-)hétérogènes, p.ex. les entrées du catalogue CSV des objets de Messier Messier.csv :

1. # Messier, NGC, Magnitude, Size [arcmin], Distance [pc], RA [h], Dec [deg], Constellation, Season, Name

2. # Attention : les données n'ont pas vocation à être très précises !

3. # D'après http://astropixels.com/messier/messiercat.html

4. M,NGC,Type,Mag,Size,Distance,RA,Dec,Con,Season,Name

5. M1,1952,Sn,8.4,5.0,1930.0,5.575,22.017,Tau,winter,Crab Nebula

6. M2,7089,Gc,6.5,12.9,11600.0,21.558,0.817,Aqr,autumn,

7. M3,5272,Gc,6.2,16.2,10400.0,13.703,28.383,CVn,spring,

8. M4,6121,Gc,5.6,26.3,2210.0,16.393,-26.533,Sco,summer,

1. >>> dt = N.dtype([('M', 'S3'),        # N° catalogue Messier

2. ...               ('NGC', 'i'),       # N° New General Catalogue

3. ...               ('Type', 'S2'),     # Code type

4. ...               ('Mag', 'f'),       # Magnitude

5. ...               ('Size', 'f'),      # Taille [arcmin]

6. ...               ('Distance', 'f'),  # Distance [pc]

7. ...               ('RA', 'f'),        # Ascension droite [h]

8. ...               ('Dec', 'f'),       # Déclinaison [deg]

9. ...               ('Con', 'S3'),      # Code constellation

10. ...               ('Season', 'S6'),   # Saison

11. ...               ('Name', 'S30')])   # Nom alternatif

12. >>> messier = N.genfromtxt("Messier.csv", dtype=dt, delimiter=',', comments='#')

13. >>> messier[1]

14. ('M1', 1952, 'Sn',  8.39999962,  5.,  1930.,  5.57499981,  22.0170002, 'Tau', 'winter', 'Crab Nebula')

15. >>> N.nanmean(messier['Mag'])

16. 7.4927273

Tableaux masqués

Le sous-module numpy.ma ajoute le support des tableaux masqués (Masked Arrays). Imaginons un tableau (4, 5) de réels (positifs ou négatifs), sur lequel nous voulons calculer pour chaque colonne la moyenne des éléments positifs uniquement :

1. >>> x = N.random.randn(4, 5); x

2. array([[-0.55867715,  1.58863893, -1.4449145 ,  1.93265481, -0.17127422],

3. [-0.86041806,  1.98317832, -0.32617721,  1.1358607 , -1.66150602],

4. [-0.88966893,  1.36185799, -1.54673735, -0.09606195,  2.23438981],

5. [ 0.35943269, -0.36134448, -0.82266202,  1.38143768, -1.3175115 ]])

6. >>> x[x >= 0]        # Donne les éléments >0 du tableau, sans leur indice

7. array([ 1.58863893,  1.93265481,  1.98317832,  1.1358607 ,  1.36185799,

8. 2.23438981,  0.35943269,  1.38143768])

9. >>> (x >= 0).nonzero() # Donne les indices ([i], [j]) des éléments positifs

10. (array([0, 0, 1, 1, 2, 2, 3, 3]), array([1, 3, 1, 3, 1, 4, 0, 3]))

11. >>> y = N.ma.masked_less(x, 0); y  # Tableau où les éléments <0 sont masqués

12. masked_array(data =

13. [[-- 1.58863892701 -- 1.93265481164 --]   # Données

14. [-- 1.98317832359 -- 1.13586070417 --]

15. [-- 1.36185798574 -- -- 2.23438980788]

16. [0.359432688656 -- -- 1.38143767743 --]],

17. mask =

18. [[ True False  True False  True]          # Bit de masquage

19. [ True False  True False  True]

20. [ True False  True  True False]

21. [False  True  True False  True]],

22. fill_value = 1e+20)

23. >>> m0 = y.mean(axis=0); m0                # Moyenne sur les lignes (axe 0)

24. masked_array(data = [0.359432688656 1.64455841211 -- 1.48331773108 2.23438980788],

25. mask = [False False  True False False],

26. fill_value = 1e+20)                 # Le résultat est un *Masked Array*

27. >>> m0.filled(-1)                          # Conversion en tableau normal

28. array([ 0.35943269,  1.64455841, -1.        ,  1.48331773,  2.23438981])

Les tableaux évolués de numpy sont parfois suffisants, mais pour une utilisation avancée, il peut être plus pertinent d'invoquer les bibliothèques dédiées Pandas & xarray.

IV-A-3 - Entrées/sorties

numpy peut lire - numpy.loadtxt() - ou sauvegarder - numpy.savetxt() - des tableaux dans un simple fichier ASCII :

1. >>> x = N.linspace(-1, 1, 100)

2. >>> N.savetxt('archive_x.dat', x)   # Sauvegarde dans le fichier 'archive_x.dat'

3. >>> y = N.loadtxt("archive_x.dat")  # Relecture à partir du fichier 'archive_x.dat'

4. >>> (x == y).all()                  # Test d'égalité stricte

5. True

numpy.loadtxt() supporte les types composés, mais ne supporte pas les données manquantes, utiliser alors la fonction numpy.genfromtxt(), plus robuste mais plus lente.

Le format texte n'est pas optimal pour de gros tableaux : il peut alors être avantageux d'utiliser le format binaire .npy, beaucoup plus compact (mais non human readable) :

1. >>> x = N.linspace(-1, 1, 1e6)

2. >>> N.save('archive_x.npy', x)   # Sauvegarde dans le fichier 'archive_x.npy'

3. >>> y = N.load("archive_x.npy")  # Relecture à partir du fichier 'archive_x.npy'

4. >>> (x == y).all()

5. True

Il est enfin possible de sérialiser les tableaux à l'aide de la bibliothèque standard [c]Pickle.

IV-A-4 - Sous-modules

numpy fournit en outre quelques fonctionnalités supplémentaires, parmi lesquelles les sous-modules suivants :

- numpy.fft : Discrete Fourier Transform ;

- numpy.random : valeurs aléatoires ;

- numpy.polynomial : manipulation des polynômes (racines, polynômes orthogonaux, etc.).

IV-A-5 - Performances

Premature optimization is the root of all evil - Donald Knuth

Même si numpy apporte un gain significatif en performance par rapport à du Python standard, il peut être possible d'améliorer la vitesse d'exécution par l'utilisation de bibliothèques externes dédiées, p.ex. :

- NumExpr est un évaluateur optimisé d'expressions numériques :

1. >>> a = N.arange(1e6)

2. >>> b = N.arange(1e6)

3. >>> %timeit a*b - 4.1*a > 2.5*b

4. 100 loops, best of 3: 11.4 ms per loop

5. >>> %timeit numexpr.evaluate("a*b - 4.1*a > 2.5*b")

6. 100 loops, best of 3: 1.97 ms per loop

7. >>> %timeit N.exp(-a)

8. 10 loops, best of 3: 60.1 ms per loop

9. >>> timeit numexpr.evaluate("exp(-a)") # Multi-threaded

10. 10 loops, best of 3: 19.3 ms per loop

- Bottleneck est une collection de fonctions accélérées, notamment pour des tableaux contenant des NaN ou pour des statistiques glissantes ;

- theano, pour optimiser l'évaluation des expressions mathématiques sur les tableaux numpy, notamment par l'utilisation des GPU et de code C généré à la volée.

Voir également Profilage et optimisation.

IV-B - Scipy

scipy est une bibliothèque numérique(5) d'algorithmes et de fonctions mathématiques, basée sur les tableaux numpy.ndarray, complétant ou améliorant (en termes de performances) les fonctionnalités de numpy.

N'oubliez pas de citer scipy & co. dans vos publications et présentations utilisant ces outils.

IV-B-1 - Tour d'horizon

- scipy.special : fonctions spéciales (fonctions de Bessel, erf, gamma, etc.).

- scipy.integrate : intégration numérique (intégration numérique ou d'équations différentielles).

- scipy.optimize : méthodes d'optimisation (minimisation, moindres carrés, zéros d'une fonction, etc.).

- scipy.interpolate : interpolation (interpolation, splines).

- scipy.fftpack : transformées de Fourier.

- scipy.signal : traitement du signal (convolution, corrélation, filtrage, ondelettes, etc.).

- scipy.linalg : algèbre linéaire.

- scipy.stats : statistiques (fonctions et distributions statistiques).

- scipy.ndimage : traitement d'images multidimensionnelles.

- scipy.io : entrées/sorties.

Liens :

- Scipy Reference

- Scipy Cookbook

Voir également :

- Scikits : modules plus spécifiques étroitement liés à scipy, parmi lesquels :

- scikit-learn : machine learning ;

- scikit-image : image processing ;

- statsmodel : modèles statistiques (tutorial).

- Scipy Topical softwares

Exercices :

Quadrature, Schéma de Romberg ★★, Méthode de Runge-Kutta ★★

IV-B-2 - Quelques exemples complets

- Interpolation

- Integration (intégrales numériques, équations différentielles)

- odeint notebook

- Zombie Apocalypse

- Optimisation (moindres carrés, ajustements, recherche de zéros)

- Traitement du signal (splines, convolution, filtrage)

- Algèbre linéaire (systèmes linéaires, moindres carrés, décompositions)

- Statistiques (variables aléatoires, distributions, tests)

IV-C - Matplotlib

Matplotlib est une bibliothèque graphique de visualisation 2D (avec un support pour la 3D, l'animation et l'interactivité), permettant des sorties de haute qualité « prêtes à publier ». C'est à la fois une bibliothèque de haut niveau, fournissant des fonctions de visualisation « clés en main » (échelle logarithmique, histogramme, courbes de niveau, etc., voir la galerie), et de bas niveau, permettant de modifier tous les éléments graphiques de la figure (titre, axes, couleurs et styles des lignes, etc., voir Anatomie d'une figure).

IV-C-1 - pylab vs. pyplot

Il existe deux interfaces pour deux types d'utilisation :

- pylab : interface procédurale, originellement très similaire à MATLAB™ et généralement réservée à l'analyse interactive :

1. >>> from pylab import *         # DÉCONSEILLÉ DANS UN SCRIPT !

2. >>> x = linspace(-pi, pi, 100)  # pylab importe numpy dans l'espace courant

3. >>> plot(x, sin(x), 'b-', label="Sinus")     # Trace la courbe y = sin(x)

4. >>> plot(x, cos(x), 'r:', label="Cosinus")   # Trace la courbe y = cos(x)

5. >>> xlabel("x [rad]")           # Ajoute le nom de l'axe des x

6. >>> ylabel("y")                 # Ajoute le nom de l'axe des y

7. >>> title("Sinus et Cosinus")   # Ajoute le titre de la figure

8. >>> legend()                    # Ajoute une légende

9. >>> savefig("simple.png")       # Enregistre la figure en PNG

- matplotlib.pyplot : interface orientée objet, préférable pour les scripts :

1. import numpy as N

2. import matplotlib.pyplot as P

4. x = N.linspace(-N.pi, N.pi, 100)

6. fig, ax = P.subplots()  # Création d'une figure contenant un seul système d'axes

7. ax.plot(x, N.sin(x), c='b', ls='-', label="Sinus")    # Courbe y = sin(x)

8. ax.plot(x, N.cos(x), c='r', ls=':', label="Cosinus")  # Courbe y = cos(x)

9. ax.set_xlabel("x [rad]")          # Nom de l'axe des x

10. ax.set_ylabel("y")                # Nom de l'axe des y

11. ax.set_title("Sinus et Cosinus")  # Titre de la figure

12. ax.legend()                       # Légende

13. fig.savefig("simple.png")         # Sauvegarde en PNG

Dans les deux cas, le résultat est le même :

Par la suite, nous nous concentrerons sur l'interface OO matplotlib.pyplot, plus puissante et flexible.

IV-C-2 - Figure et axes

L'élément de base est le système d'axes matplotlib.axes.Axes, qui définit et réalise la plupart des éléments graphiques (tracé de courbes, définition des axes, annotations, etc.). Un ou plusieurs de ces systèmes d'axes sont regroupés au sein d'une matplotlib.figure.Figure.

Ainsi, pour générer une figure contenant 2 (vertical) × 3 (horizontal) = 6 systèmes d'axes (numérotés de 1 à 6) :

1. fig = P.figure()

2. for i in range(1, 4):

3. ax = fig.add_subplot(2, 3, i, xticks=[], yticks=[])

4. ax.text(0.5, 0.5, "subplot(2, 3, {})".format(i),

5. ha='center', va='center', size='large')

6. for i in range(3, 5):

7. ax = fig.add_subplot(2, 2, i, xticks=[], yticks=[])

8. ax.text(0.5, 0.5, "subplot(2, 2, {})".format(i),

9. ha='center', va='center', size='large')

Pour des mises en page plus complexes, il est possible d'utiliser le kit gridspec, p.ex. :

1. from matplotlib.gridspec import GridSpec

3. fig = P.figure()

4. fig.suptitle("grid = GridSpec(2, 3)", fontsize='x-large')

6. grid = GridSpec(2, 3)

7. ax1 = fig.add_subplot(grid[0, :-1], xticks=[], yticks=[])

8. ax1.text(0.5, 0.5, "grid[0, :-1]", ha='center', va='center', size='large')

9. ax2 = fig.add_subplot(grid[:, -1], xticks=[], yticks=[])

10. ax3.text(0.5, 0.5, "grid[:, -1]", ha='center', va='center', size='large')

11. ax3 = fig.add_subplot(grid[1, 0], xticks=[], yticks=[])

12. ax3.text(0.5, 0.5, "grid[1, 0]", ha='center', va='center', size='large')

13. ax4 = fig.add_subplot(grid[1, 1], xticks=[], yticks=[])

14. ax4.text(0.5, 0.5, "grid[1, 1]", ha='center', va='center', size='large')

Enfin, il est toujours possible (mais peu pratique) de créer soi-même le système d'axes dans les coordonnées relatives à la figure :

1. fig = P.figure()

2. ax0 = fig.add_axes([0, 0, 1, 1], frameon=False,

3. xticks=N.linspace(0, 1, 11), yticks=N.linspace(0, 1, 11))

4. ax0.grid(True, ls='--')

5. ax1 = fig.add_axes([0.1, 0.2, 0.8, 0.6], xticks=[], yticks=[], fc='0.9')

6. ax1.text(0.5, 0.5, "[0.1, 0.2, 0.8, 0.6]", ha='center', va='center', size='large')

7. ax2 = fig.add_axes([0.2, 0.3, 0.4, 0.1], xticks=[], yticks=[], fc='0.8')

8. ax2.text(0.5, 0.5, "[0.2, 0.3, 0.4, 0.1]", ha='center', va='center', size='large')

IV-C-3 - Sauvegarde et affichage interactif

La méthode matplotlib.figure.Figure.savefig() permet de sauvegarder la figure dans un fichier dont le format est automatiquement défini par son extension, png (raster), [e]ps, pdf, svg (vector), etc., via différents backends.

Il est également possible d'afficher la figure dans une fenêtre interactive avec la commande matplotlib.pyplot.show() :

Utiliser ipython --pylab pour l'utilisation interactive des figures dans une session ipython.

IV-C-4 - Anatomie d'une figure

L'interface OO matplotlib.pyplot donne accès à tous les éléments d'une figure (titre, axes, légende, etc.), qui peuvent alors être ajustés (couleur, police, taille, etc.).

Anatomie d'une figure.

N'oubliez pas de citer matplotlib (notamment [Matplotlib07]) dans vos publications et présentations utilisant cet outil.

Liens :

- User's Guide

- Gallery

- Tutorial matplotlib

- Tutoriel matplotlib

Voir également :

- MPLD3, un backend matplotlib interactif basé sur la bibliothèque web 3D.js ;

- basemap et cartopy, bibliothèques de cartographie sphérique ;

- Seaborn, une surcouche de visualisation statistique à matplotlib et Pandas & xarray ;

- HoloViews, une surcouche de visualisation et d'analyse à matplotlib ;

- ggplot, une surcouche orientée Grammar of Graphics à matplotlib ;

- Bokeh, une bibliothèque graphique alternative à matplotlib plus orientée web/temps réel.

Exemples :

figure.py, filtres2ndOrdre.py

Exercices :

Quartet d'Anscombe ★, Ensemble de Julia ★★, Diagramme de bifurcation : la suite logistique ★★

IV-C-5 - Visualisation 3D

Matplotlib fournit d'emblée une interface mplot3d pour des figures 3D assez simples :

Exemple de figure matplotlib 3D.

Pour des visualisations plus complexes, mayavi.mlab est une bibliothèque graphique de visualisation 3D s'appuyant sur Mayavi.

Imagerie par résonance magnétique.

N'oubliez pas de citer mayavi dans vos publications et présentations utilisant cet outil.

Voir également :

- VPython : 3D Programming for Ordinary Mortals ;

- Glowscript : VPython dans le navigateur.

V - Bibliothèques scientifiques avancées

V-A - Pandas et xarray

P est une bibliothèque pour la structuration et l'analyse avancée de données hétérogènes (PANel DAta). Elle fournit notamment :

- des structures de données relationnelles (« labellisées »), avec une indexation simple ou hiérarchique (c.-à-d. à plusieurs niveaux) ;

- des méthodes d'alignement et d'agrégation des données, avec gestion des données manquantes ;

- un support performant des labels temporels (p.ex. dates, de par son origine dans le domaine de l'économétrie), et des statistiques « glissantes » ;

- de nombreuses fonctions d'entrée/sortie, d'analyse statistiques et de visualisation.

Les fonctionnalités de pandas sont très riches et couvrent de nombreux aspects (données manquantes, dates, analyse statistiques, etc.) : il n'est pas question de toutes les aborder ici. Avant de vous lancer dans une manipulation qui vous semble complexe, bien inspecter la documentation, très complète (p.ex. les recettes du cookbook), pour vérifier qu'elle n'est pas déjà implémentée ou documentée, ou pour identifier l'approche la plus efficace.

La convention utilisée ici est « import pandas as PD ».

Les bibliothèques pandas et xarray sont encore en phase de développement assez intense, et vont probablement être amenées à évoluer significativement, et pas nécessairement de manière rétrocompatible. Nous travaillons ici sur les versions :

- pandas 0.20.x

- xarray 0.9.x

V-A-1 - Structures

Références : Introduction to Data Structures

Pandas dispose de deux grandes structures de données(6) :

Nom de la structure

Rang

Description

pandas.Series

1

Vecteur de données homogènes labellisées

pandas.DataFrame

2

Tableau structuré de colonnes homogènes

1. >>> PD.Series(range(3))  # Série d'entiers sans indexation

2. 0    0

3. 1    1

4. 2    2

5. dtype: int64

6. >>> PD.Series(N.random.randn(3), index=list('abc'))  # Série de réels indexés

7. a   -0.553480

8. b    0.081297

9. c   -1.845835

10. dtype: float64

11. >>> PD.DataFrame(N.random.randn(3, 4))

12. 0         1         2         3

13. 0  1.570977 -0.677845  0.094364 -0.362031

14. 1 -0.136712  0.762300  0.068611  1.265816

15. 2 -0.697760  0.791288  0.449645 -1.105062

16. >>> PD.DataFrame([(1, N.exp(1), 'un'), (2, N.exp(2), 'deux'), (3, N.exp(3), 'trois')],

17. ...              index=list('abc'), columns='i val nom'.split())

18. i        val    nom

19. a  1   2.718282     un

20. b  2   7.389056   deux

21. c  3  20.085537  trois

Pour mettre en évidence la puissance de Pandas, nous utiliserons le catalogue des objets Messier vu précédemment. Le fichier peut être importé à l'aide de la fonction pandas.read_csv(), et le dataframe résultant est labellisé à la volée par la colonne M (pandas.DataFrame.set_index()) :

1. >>> messier = PD.read_csv("Messier.csv", comment='#')  # Lecture du fichier CSV

2. >>> messier.set_index('M', inplace=True)  # Indexation sur la colonne "M"

3. >>> messier.info()                        # Informations générales

4. <class 'pandas.core.frame.DataFrame'>

5. Index: 110 entries, M1 to M110

6. Data columns (total 10 columns):

7. NGC         108 non-null object

8. Type        110 non-null object

9. Mag         110 non-null float64

10. Size        110 non-null float64

11. Distance    110 non-null float64

12. RA          110 non-null float64

13. Dec         110 non-null float64

14. Con         110 non-null object

15. Season      110 non-null object

16. Name        31 non-null object

17. dtypes: float64(5), object(5)

18. memory usage: 9.5+ KB

19. >>> messier.head(3)  # Par défaut les 5 premières lignes

20. NGC Type  Mag  Size  Distance      RA     Dec  Con  Season         Name

21. M

22. M1  1952   Sn  8.4   5.0    1930.0   5.575  22.017  Tau  winter  Crab Nebula

23. M2  7089   Gc  6.5  12.9   11600.0  21.558   0.817  Aqr  autumn          NaN

24. M3  5272   Gc  6.2  16.2   10400.0  13.703  28.383  CVn  spring          NaN

Un dataframe est caractérisé par son indexation pandas.DataFrame.index et ses colonnes pandas.DataFrame.columns (de type pandas.Index ou pandas.MultiIndex), et les valeurs des données pandas.DataFrame.values :

1. >>> messier.index    # Retourne un Index

2. Index([u'M1', u'M2', u'M3', ..., u'M108', u'M109', u'M110'],

3. dtype='object', name=u'M', length=110)

4. >>> messier.columns  # Retourne un Index

5. Index([u'NGC', u'Type', u'Mag', ..., u'Con', u'Season', u'Name'],

6. dtype='object')

7. >>> messier.dtypes   # Retourne une Series indexée sur le nom des colonnes

8. NGC          object

9. Type         object

10. Mag         float64

11. Size        float64

12. Distance    float64

13. RA          float64

14. Dec         float64

15. Con          object

16. Season       object

17. Name         object

18. dtype: object

19. >>> messier.values

20. array([['1952', 'Sn', 8.4, ..., 'Tau', 'winter', 'Crab Nebula'],

21. ['7089', 'Gc', 6.5, ..., 'Aqr', 'autumn', nan],

22. ...,

23. ['3992', 'Ba', 9.8, ..., 'UMa', 'spring', nan],

24. ['205', 'El', 8.5, ..., 'And', 'autumn', nan]], dtype=object)

25. >>> messier.shape

26. (110, 10)

Une description statistique sommaire des colonnes numériques est obtenue par pandas.DataFrame.describe() :

1. >>> messier.drop(['RA', 'Dec'], axis=1).describe()

2. Mag        Size      Distance

3. count  110.000000  110.000000  1.100000e+02

4. mean     7.492727   17.719091  4.574883e+06

5. std      1.755657   22.055100  7.141036e+06

6. min      1.600000    0.800000  1.170000e+02

7. 25%      6.300000    6.425000  1.312500e+03

8. 50%      7.650000    9.900000  8.390000e+03

9. 75%      8.900000   17.300000  1.070000e+07

10. max     10.200000  120.000000  1.840000e+07

V-A-2 - Accès aux données

Référence : Indexing and Selecting Data

L'accès par colonne retourne une pandas.Series (avec la même indexation) pour une colonne unique, ou un nouveau pandas.DataFrame pour plusieurs colonnes :

1. >>> messier.NGC  # Équivalent à messier['NGC']

2. M

3. M1      1952

4. M2      7089

5. ...

6. M109    3992

7. M110     205

8. Name: NGC, Length: 110, dtype: object

9. >>> messier[['RA', 'Dec']]  # = messier.filter(items=('RA', 'Dec'))

10. RA     Dec

11. M

12. M1     5.575  22.017

13. M2    21.558   0.817

14. ...      ...     ...

15. M109  11.960  53.383

16. M110   0.673  41.683

17. [110 rows x 2 columns]

L'accès par slice retourne un nouveau dataframe :

1. >>> messier[:6:2]       # Lignes 0 (inclus) à 6 (exclu) par pas de 2

2. NGC Type  Mag  Size  Distance      RA     Dec  Con  Season         Name

3. M

4. M1  1952   Sn  8.4   5.0    1930.0   5.575  22.017  Tau  winter  Crab Nebula

5. M3  5272   Gc  6.2  16.2   10400.0  13.703  28.383  CVn  spring          NaN

6. M5  5904   Gc  5.6  17.4    7520.0  15.310   2.083  Ser  summer          NaN

L'accès peut également se faire par labels via pandas.DataFrame.loc :

1. >>> messier.loc['M31']  # Retourne une Series indexée par les noms de colonne

2. NGC                    224

3. Type                    Sp

4. ...

5. Season              autumn

6. Name      Andromeda Galaxy

7. Name: M31, Length: 10, dtype: object

8. >>> messier.loc['M31', ['Type', 'Name']]           # Retourne une Series

9. Type                  Sp

10. Name    Andromeda Galaxy

11. Name: M31, dtype: object

12. >>> messier.loc[['M31', 'M51'], ['Type', 'Name']]  # Retourne un DataFrame

13. Type              Name

14. M

15. M31   Sp  Andromeda Galaxy

16. M51   Sp  Whirlpool Galaxy

17. >>> messier.loc['M31':'M33', ['Type', 'Name']]     # De M31 à M33 inclu

18. Type               Name

19. M

20. M31   Sp   Andromeda Galaxy

21. M32   El                NaN

22. M33   Sp  Triangulum Galaxy

De façon symétrique, l'accès peut se faire par position (n° de ligne/colonne) via pandas.DataFrame.iloc, p.ex. :

1. >>> messier.iloc[30:33, [1, -1]]  # Ici, l'indice 33 n'est PAS inclus !

2. Type               Name

3. M

4. M31   Sp   Andromeda Galaxy

5. M32   El                NaN

6. M33   Sp  Triangulum Galaxy

7. >>> messier.iloc[30:33, messier.columns.get_loc('Name')]  # Mix des 2 approches

8. M

9. M31     Andromeda Galaxy

10. M32                  NaN

11. M33    Triangulum Galaxy

12. Name: Name, dtype: object

Les fonctions pandas.DataFrame.at et pandas.DataFrame.iat permettent d'accéder rapidement aux données individuelles :

1. >>> messier.at['M31', 'NGC']  # 20× plus rapide que messier.loc['M31']['NGC']

2. '224'

3. >>> messier.iat[30, 0]        # 20× plus rapide que messier.iloc[0][0]

4. '224'

Noter qu'il existe une façon de filtrer les données sur les colonnes ou les labels :

1. >>> messier.filter(regex='M.7', axis='index').filter('RA Dec'.split())

2. RA     Dec

3. M

4. M17  18.347 -16.183

5. M27  19.993  22.717

6. M37   5.873  32.550

7. M47   7.610 -14.500

8. M57  18.893  33.033

9. M67   8.840  11.817

10. M77   2.712   0.033

11. M87  12.513  12.400

12. M97  11.247  55.017

Comme pour numpy, il est possible d'opérer une sélection booléenne :

1. >>> messier.loc[messier['Con'] == 'UMa', ['NGC', 'Name']]

2. NGC             Name

3. M

4. M40   Win4       Winnecke 4

5. M81   3031    Bode's Galaxy

6. M82   3034     Cigar Galaxy

7. M97   3587       Owl Nebula

8. M101  5457  Pinwheel Galaxy

9. M108  3556              NaN

10. M109  3992              NaN

11. >>> messier[messier['Season'].isin('winter spring'.split())].head(3)

12. NGC Type  Mag  Size  Distance      RA     Dec  Con  Season         Name

13. M

14. M1   1952   Sn  8.4   5.0    1930.0   5.575  22.017  Tau  winter  Crab Nebula

15. M3   5272   Gc  6.2  16.2   10400.0  13.703  28.383  CVn  spring          NaN

16. M35  2168   Oc  5.3  28.0     859.0   6.148  24.333  Gem  winter          NaN

17. >>> messier.loc[lambda df: N.radians(df.Size / 60) * df.Distance < 1].Name

18. M

19. M27           Dumbbell Nebula

20. M40                Winnecke 4

21. M57               Ring Nebula

22. M73                       NaN

23. M76    Little Dumbbell Nebula

24. M78                       NaN

25. M97                Owl Nebula

26. Name: Name, dtype: object

27. >>> messier.query("(Mag < 5) & (Size > 60)").Name

28. M

29. M7          Ptolemy's Cluster

30. M24    Sagittarius Star Cloud

31. M31          Andromeda Galaxy

32. M42     Great Nebula in Orion

33. M44           Beehive Cluster

34. M45                  Pleiades

35. Name: Name, dtype: object

Sélection

Syntaxe

Résultat

Colonne unique

df.col or df[col]

pandas.Series

Liste de colonnes

df[[c1, ...]]

pandas.DataFrame

Lignes par tranche

df[slice]

pandas.DataFrame

Label unique

df.loc[label]

pandas.Series

Liste de labels

df.loc[[lab1, ...]]

pandas.DataFrame

Labels par tranche

df.loc[lab1:lab2]

pandas.DataFrame

Ligne entière par n°

df.iloc[i]

pandas.Series

Ligne partielle par n°

df.iloc[i, [j,...]]

pandas.Series

Valeur par label

df.at[lab, col]

Scalaire

Valeur par n°

df.iat[i, j]

Scalaire

Ligne par sél. booléenne

df.loc[sel] or df[sel] or df.query("sel")

pandas.DataFrame

pandas.DataFrame.drop() permet d'éliminer une ou plusieurs colonnes d'un dataframe :

1. >>> messier.drop(['RA', 'Dec'], axis=1).head(3)  # Élimination de colonnes

2. NGC Type  Mag  Size  Distance  Con  Season  Name

3. M

4. M1  1952   Sn  8.4   5.0    1930.0  Tau  winter  Crab Nebula

5. M2  7089   Gc  6.5  12.9   11600.0  Aqr  autumn          NaN

6. M3  5272   Gc  6.2  16.2   10400.0  CVn  spring          NaN

pandas.DataFrame.dropna() et pandas.DataFrame.fillna() permettent de gérer les données manquantes (NaN) :

1. >>> messier.dropna(axis=0, how='any', subset=['NGC', 'Name']).head(3)

2. NGC Type  Mag  Size  Distance      RA     Dec  Con  Season               Name

3. M

4. M1  1952   Sn  8.4   5.0    1930.0   5.575  22.017  Tau  winter        Crab Nebula

5. M6  6405   Oc  4.2  25.0     491.0  17.668 -32.217  Sco  summer  Butterfly Cluster

6. M7  6475   Oc  3.3  80.0     245.0  17.898 -34.817  Sco  summer  Ptolemy's Cluster

7. >>> messier.fillna('', inplace=True)  # Remplace les NaN à la volée

8. >>> messier.head(3)

9. NGC Type  Mag  Size  Distance      RA     Dec  Con  Season         Name

10. M

11. M1  1952   Sn  8.4   5.0    1930.0   5.575  22.017  Tau  winter  Crab Nebula

12. M2  7089   Gc  6.5  12.9   11600.0  21.558   0.817  Aqr  autumn

13. M3  5272   Gc  6.2  16.2   10400.0  13.703  28.383  CVn  spring

Référence : Working with missing data

Attention par défaut, beaucoup d'opérations retournent une copie de la structure, sauf si l'opération se fait « sur place » (inplace=True). D'autres opérations d'accès retournent seulement une nouvelle vue des mêmes données.

1. >>> df = PD.DataFrame(N.arange(12).reshape(3, 4),

2. ...                   index=list('abc'), columns=list('ABCD')); df

3. A  B   C   D

4. a  0  1   2   3

5. b  4  5   6   7

6. c  8  9  10  11

7. >>> df.drop('a', axis=0)

8. A  B   C   D

9. b  4  5   6   7

10. c  8  9  10  11

11. >>> colA = df['A']  # Nouvelle vue de la colonne 'A'

12. >>> colA += 1       # Opération sur place

13. >>> df              # la ligne 'a' est tjs là, la colonne 'A' a été modifiée

14. A  B   C   D

15. a  1  1   2   3

16. b  5  5   6   7

17. c  9  9  10  11

Lien : Returning a view versus a copy

Indexation hiérarchique

Références : Multiindex / Advanced Indexing

pandas.MultiIndex offre une indexation hiérarchique, permettant de stocker et manipuler des données avec un nombre arbitraire de dimensions dans des structures plus simples.

1. >>> saisons = messier.reset_index()                     # Élimine l'indexation actuelle

2. >>> saisons.set_index(['Season', 'Type'], inplace=True) # MultiIndex

3. >>> saisons.head(3)

4. M   NGC  Mag  Size  Distance      RA     Dec  Con         Name

5. Season Type

6. winter Sn    M1  1952  8.4   5.0    1930.0   5.575  22.017  Tau  Crab Nebula

7. autumn Gc    M2  7089  6.5  12.9   11600.0  21.558   0.817  Aqr

8. spring Gc    M3  5272  6.2  16.2   10400.0  13.703  28.383  CVn

Les informations contenues sont toujours les mêmes, mais structurées différemment :

1. >>> saisons.loc['spring'].head(3)  # Utilisation du 1er label

2. M   NGC  Mag  Size    Distance      RA     Dec  Con        Name

3. Type

4. Gc     M3  5272  6.2  16.2     10400.0  13.703  28.383  CVn

5. Ds    M40  Win4  8.4   0.8       156.0  12.373  58.083  UMa  Winnecke 4

6. El    M49  4472  8.4   8.2  18400000.0  12.497   8.000  Vir

7. >>> saisons.loc['spring', 'El'].head(3)  # Utilisation des 2 labels

8. M   NGC  Mag  Size    Distance      RA    Dec  Con Name

9. Season Type

10. spring El    M49  4472  8.4   8.2  18400000.0  12.497   8.00  Vir

11. El    M59  4621  9.6   4.2  18400000.0  12.700  11.65  Vir

12. El    M60  4649  8.8   6.5  18400000.0  12.728  11.55  Vir

La fonction pandas.DataFrame.xs() permet des sélections sur les différents niveaux d'indexation :

1. >>> saisons.xs('spring', level='Season').head(3) # = saisons.loc['spring']

2. M   NGC  Mag  Size    Distance      RA     Dec  Con        Name

3. Type

4. Gc     M3  5272  6.2  16.2     10400.0  13.703  28.383  CVn

5. Ds    M40  Win4  8.4   0.8       156.0  12.373  58.083  UMa  Winnecke 4

6. El    M49  4472  8.4   8.2  18400000.0  12.497   8.000  Vir

7. >>> saisons.xs('El', level='Type').head(3) # Sélection sur le 2e niveau

8. M   NGC  Mag  Size    Distance      RA     Dec  Con Name

9. Season

10. autumn  M32   221  8.1   7.0    920000.0   0.713  40.867  And

11. spring  M49  4472  8.4   8.2  18400000.0  12.497   8.000  Vir

12. spring  M59  4621  9.6   4.2  18400000.0  12.700  11.650  Vir

Le (multi)index n'est pas nécessairement trié à sa création, pandas.sort_index() permet d'y remédier :

1. >>> saisons[['M', 'NGC', 'Name']].head()

2. M   NGC         Name

3. Season Type

4. winter Sn    M1  1952  Crab Nebula

5. autumn Gc    M2  7089

6. spring Gc    M3  5272

7. summer Gc    M4  6121

8. Gc    M5  5904

9. >>> saisons[['M', 'NGC', 'Name']].sort_index().head()

10. M   NGC                    Name

11. Season Type

12. autumn El     M32   221

13. El    M110   205

14. Gc      M2  7089

15. Gc     M15  7078  Great Pegasus Globular

16. Gc     M30  7099

V-A-3 - Manipulation des données

Références : Essential Basic Functionality

Comme dans numpy, il est possible de modifier les valeurs, ajouter/retirer des colonnes ou des lignes, tout en gérant les données manquantes.

l'interopérabilité entre pandas et numpy est totale, toutes les fonctions Numpy peuvent prendre une structure Pandas en entrée, et s'appliquer aux colonnes appropriées :

1. >>> N.mean(messier, axis=0)  # Moyenne sur les lignes ? Series indexée sur les colonnes

2. Mag         7.492727e+00

3. Size        1.771909e+01

4. Distance    4.574883e+06

5. RA          1.303774e+01

6. Dec         9.281782e+00

7. dtype: float64

1. >>> N.random.seed(0)

2. >>> df = PD.DataFrame(

3. {'one': PD.Series(N.random.randn(3), index=['a', 'b', 'c']),

4. 'two': PD.Series(N.random.randn(4), index=['a', 'b', 'c', 'd']),

5. 'three': PD.Series(N.random.randn(3), index=['b', 'c', 'd'])})

6. >>> df

7. one     three       two

8. a  1.764052       NaN  2.240893

9. b  0.400157 -0.151357  1.867558

10. c  0.978738 -0.103219 -0.977278

11. d       NaN  0.410599  0.950088

12. >>> df['four'] = df['one'] + df['two']; df  # Création d'une nouvelle colonne

13. one     three       two      four

14. a  1.764052       NaN  2.240893  4.004946

15. b  0.400157 -0.151357  1.867558  2.267715

16. c  0.978738 -0.103219 -0.977278  0.001460

17. d       NaN  0.410599  0.950088       NaN

18. >>> df.sub(df.loc['b'], axis='columns') # Soustraction d'une ligne à toutes les colonnes (axis=1)

19. one     three       two      four

20. a  1.363895       NaN  0.373335  1.737230

21. b  0.000000  0.000000  0.000000  0.000000

22. c  0.578581  0.048138 -2.844836 -2.266255

23. d       NaN  0.561956 -0.917470       NaN

24. >>> df.sub(df['one'], axis='index') # Soustraction d'une colonne à toutes les lignes (axis=0 ou 'rows')

25. one     three       two      four

26. a  0.0       NaN  0.476841  2.240893

27. b  0.0 -0.551514  1.467401  1.867558

28. c  0.0 -1.081957 -1.956016 -0.977278

29. d  NaN       NaN       NaN       NaN

1. >>> df.sort_values(by='a', axis=1) # Tri des colonnes selon les valeurs de la ligne 'a'

2. one       two     three

3. a  1.764052  2.240893       NaN

4. b  0.400157  1.867558 -0.151357

5. c  0.978738 -0.977278 -0.103219

6. d       NaN  0.950088  0.410599

7. >>> df.min(axis=1) # Valeur minimale le long des colonnes

8. a    1.764052

9. b   -0.151357

10. c   -0.977278

11. d    0.410599

12. dtype: float64

13. >>> df.idxmin(axis=1) # Colonne des valeurs minimales le long des colonnes

14. a      one

15. b    three

16. c      two

17. d    three

18. dtype: object

1. >>> df.mean(axis=0)  # Moyenne sur toutes les lignes (gestion des données manquantes)

2. one      1.047649

3. three    0.052007

4. two      1.020315

5. dtype: float64

Si les bibliothèques d'optimisation de performances Bottleneck et NumExpr sont installées, pandas en bénéficiera de façon transparente.

V-A-4 - Regroupement et agrégation de données

Histogramme et discrétisation

Compter les objets Messier par constellation avec pandas.value_counts() :

1. >>> PD.value_counts(messier['Con']).head(3)

2. Sgr    15

3. Vir    11

4. Com     8

5. Name: Con, dtype: int64

Partitionner les objets en trois groupes de magnitude (par valeurs : pandas.cut(), par quantiles : pandas.qcut()), et les compter :

1. >>> PD.value_counts(PD.cut(messier['Mag'], 3)).sort_index() # Par valeurs

2. (1.591, 4.467]     6

3. (4.467, 7.333]    40

4. (7.333, 10.2]     64

5. Name: Mag, dtype: int64

6. >>> PD.value_counts(PD.qcut(messier['Mag'], [0, .3, .6, 1])).sort_index()  # Par quantiles

7. (1.599, 6.697]    36

8. (6.697, 8.4]      38

9. (8.4, 10.2]       36

10. Name: Mag, dtype: int64

Group-by

Référence : Group By: split-apply-combine

Pandas offre la possibilité de regrouper les données selon divers critères (pandas.DataFrame.groupby()), de les agréger au sein de ces groupes et de stocker les résultats dans une structure avec indexation hiérarchique (pandas.DataFrame.agg()).

1. >>> seasonGr = messier.groupby('Season')  # Retourne un DataFrameGroupBy

2. >>> seasonGr.groups

3. {'autumn': Index([u'M2', u'M15', ..., u'M103', u'M110'],

4. dtype='object', name=u'M'),

5. 'spring': Index([u'M3', u'M40', ..., u'M108', u'M109'],

6. dtype='object', name=u'M'),

7. 'summer': Index([u'M4', u'M5', ..., u'M102', u'M107'],

8. dtype='object', name=u'M'),

9. 'winter': Index([u'M1', u'M35', ..., u'M79', u'M93'],

10. dtype='object', name=u'M')}

11. >>> seasonGr.size()

12. Season

13. autumn    14

14. spring    38

15. summer    40

16. winter    18

17. dtype: int64

18. >>> seasonGr.get_group('winter').head(3)

19. Con     Dec  Distance  Mag   NGC         Name     RA  Size Type

20. M

21. M1   Tau  22.017    1930.0  8.4  1952  Crab Nebula  5.575   5.0   Sn

22. M35  Gem  24.333     859.0  5.3  2168               6.148  28.0   Oc

23. M36  Aur  34.133    1260.0  6.3  1960               5.602  12.0   Oc

24. >>> seasonGr['Size'].agg([N.mean, N.std])  # Taille moyenne et stddev par groupe

25. mean        std

26. Season

27. autumn  24.307143  31.472588

28. spring   7.197368   4.183848

29. summer  17.965000  19.322400

30. winter  34.261111  29.894779

31. >>> seasonGr.agg({'Size': N.max, 'Mag': N.min})

32. Mag   Size

33. Season

34. autumn  3.4  120.0

35. spring  6.2   22.0

36. summer  3.3   90.0

37. winter  1.6  110.0

1. >>> magGr = messier.groupby(

2. ...     [PD.qcut(messier['Mag'], [0, .3, .6, 1],

3. ...              labels='Bright Medium Faint'.split()),

4. ...     "Season"])

5. >>> magGr['Mag', 'Size'].agg({'Mag': ['count', 'mean'],

6. ...                           'Size': [N.mean, N.std]})

7. Mag                 Size

8. count      mean       mean        std

9. Mag    Season

10. Bright autumn     6  5.316667  45.200000  40.470878

11. spring     1  6.200000  16.200000        NaN

12. summer    15  5.673333  30.840000  26.225228

13. winter    13  5.138462  42.923077  30.944740

14. Faint  autumn     4  9.225000   8.025000   4.768910

15. spring    30  9.236667   5.756667   2.272578

16. summer     7  8.971429   7.814286   9.135540

17. winter     3  8.566667   9.666667   6.429101

18. Medium autumn     4  7.500000   9.250000   3.304038

19. spring     7  7.714286  12.085714   5.587316

20. summer    18  7.366667  11.183333   4.825453

21. winter     2  7.550000  14.850000   8.697413

Tableau croisé (Pivot table)

Référence : Reshaping and Pivot Tables

Calculer la magnitude et la taille moyennes des objets Messier selon leur type avec pandas.DataFrame.pivot_table() :

1. >>> messier['Mag Size Type'.split()].pivot_table(columns='Type')

2. Type   As    Ba         Di  ...     Pl   Sn         Sp

3. Mag   9.0  9.85   7.216667  ...  9.050  8.4   8.495652

4. Size  2.8  4.80  33.333333  ...  3.425  5.0  15.160870

V-A-5 - Visualisations

Exemple :

Démonstration Pandas/Seaborn (notebook : pokemon.ipynb) sur le jeu de données : Pokemon.csv.

Références :

- Visualization

- Seaborn: statistical data visualization

Autres exemples de visualisation de jeux de données complexes (utilisation de pandas et seaborn)

- Iris Dataset

- Histoire des sujets télévisés

Liens

- Pandas Tutorial

- Pandas Cookbook

- Pandas Lessons for New Users

- Practical Data Analysis

- Modern Pandas

- Various tutorials

Exercices :

- Exercices for New Users

V-A-6 - Xarray

xarray est une bibliothèque pour la structuration de données homogènes de dimension arbitraire. Suivant la philosophie de la bibliothèque Pandas dont elle est issue (et dont elle dépend), elle permet notamment de nommer les différentes dimensions (axes) des tableaux (p.ex. x.sum(axis='time')), d'indexer les données (p.ex. x.loc['M31']), de naturellement gérer les opérations de broadcasting, des opérations de regroupement et d'agrégation (p.ex. x.groupby(time.dayofyear).mean()), une gestion plus facile des données manquantes et d'alignement des tableaux indexés (p.ex. align(x, y, join='outer')).

pandas excels at working with tabular data. That suffices for many statistical analyses, but physical scientists rely on N-dimensional arrays - which is where xarray comes in.

xarray fournit deux structures principales, héritées du format netCDF :

- xarray.DataArray, un tableau N-D indexé généralisant le pandas.Series ;

- xarray.DataSet, un dictionnaire regroupant plusieurs DataArray alignés selon une ou plusieurs dimensions, et similaire au pandas.DataFrame.

La convention utilisée ici est « import xarray as X ».

1. >>> N.random.seed(0)

2. >>> data = X.DataArray(N.arange(3*4, dtype=float).reshape((3, 4)), # Tableau de données

3. ...                    dims=('x', 'y'),            # Nom des dimensions

4. ...                    coords={'x': list('abc')},  # Indexation des coordonnées en 'x'

5. ...                    name='mesures',             # Nom du tableau

6. ...                    attrs=dict(author='Y. Copin'))  # Métadonnées

7. >>> data

8. <xarray.DataArray 'mesures' (x: 3, y: 4)>

9. array([[ 0.,  1.,  2.,  3.],

10. [ 4.,  5.,  6.,  7.],

11. [ 8.,  9., 10., 11.]])

12. Coordinates:

13. * x        (x) |S1 'a' 'b' 'c'

14. Dimensions without coordinates: y

15. Attributes:

16. author:     Y. Copin

17. >>> data.to_pandas()     # Conversion en DataFrame à indexation simple

18. y    0    1     2     3

19. x

20. a  0.0  1.0   2.0   3.0

21. b  4.0  5.0   6.0   7.0

22. c  8.0  9.0  10.0  11.0

23. >>> data.to_dataframe()  # Conversion en DataFrame multi-indexé (hiérarchique)

24. mesures

25. x y

26. a 0      0.0

27. 1      1.0

28. 2      2.0

29. 3      3.0

30. b 0      4.0

31. 1      5.0

32. 2      6.0

33. 3      7.0

34. c 0      8.0

35. 1      9.0

36. 2     10.0

37. 3     11.0

38. >>> data.dims

39. ('x', 'y')

40. >>> data.coords

41. Coordinates:

42. * x        (x) |S1 'a' 'b' 'c'

43. >>> data.values

44. array([[  0.,   1.,   2.,   3.],

45. [  4.,   5.,   6.,   7.],

46. [  8.,   9.,  10.,  11.]])

47. >>> data.attrs

48. OrderedDict([('author', 'Y. Copin')])

1. >>> data[:, 1]             # Accès par indices

2. <xarray.DataArray 'mesures' (x: 3)>

3. array([ 1.,  5.,  9.])

4. Coordinates:

5. * x        (x) |S1 'a' 'b' 'c'

6. >>> data.loc['a':'b', -1]  # Accès par labels

7. <xarray.DataArray 'mesures' (x: 2)>

8. array([ 3.,  7.])

9. Coordinates:

10. * x        (x) |S1 'a' 'b'

11. >>> data.sel(x=['a', 'c'], y=2)

12. <xarray.DataArray 'mesures' (x: 2)>

13. array([  2.,  10.])

14. Coordinates:

15. * x        (x) |S1 'a' 'c'

1. >>> data.mean(dim='x')  # Moyenne le long de la dimension 'x' = data.mean(axis=0)

2. <xarray.DataArray 'mesures' (y: 4)>

3. array([ 4.,  5.,  6.,  7.])

4. Dimensions without coordinates: y

1. >>> data2 = X.DataArray(N.arange(6).reshape(2, 3) * 10,

2. ...                     dims=('z', 'x'), coords={'x': list('abd')})

3. >>> data2.to_pandas()

4. x   a   b   d

5. z

6. 0   0  10  20

7. 1  30  40  50

8. >>> data.to_pandas()  # REMINDER

9. y    0    1     2     3

10. x

11. a  0.0  1.0   2.0   3.0

12. b  4.0  5.0   6.0   7.0

13. c  8.0  9.0  10.0  11.0

14. >>> data2.values + data.values  # Opération sur des tableaux numpy incompatibles

15. ValueError: operands could not be broadcast together with shapes (2,3) (3,4)

16. >>> data2 + data                # Alignement automatique sur les coord. communes !

17. <xarray.DataArray (z: 2, x: 2, y: 4)>

18. array([[[  0.,   1.,   2.,   3.],

19. [ 14.,  15.,  16.,  17.]],

20. [[ 30.,  31.,  32.,  33.],

21. [ 44.,  45.,  46.,  47.]]])

22. Coordinates:

23. * x        (x) object 'a' 'b'

24. Dimensions without coordinates: z, y

1. >>> data['isSmall'] = data.sum(dim='y') < 10; data  # Booléen "Somme sur y < 10"

2. <xarray.DataArray 'mesures' (x: 3, y: 4)>

3. array([[  0.,   1.,   2.,   3.],

4. [  4.,   5.,   6.,   7.],

5. [  8.,   9.,  10.,  11.]])

6. Coordinates:

7. * x        (x) |S1 'a' 'b' 'c'

8. isSmall  (x) bool True False False

9. Dimensions without coordinates: y

10. >>> data.groupby('isSmall').mean(dim='x')  # Regroupement et agrégation

11. <xarray.DataArray 'mesures' (isSmall: 2, y: 4)>

12. array([[ 6.,  7.,  8.,  9.],

13. [ 0.,  1.,  2.,  3.]])

14. Coordinates:

15. * isSmall  (isSmall) object False True

16. Dimensions without coordinates: y

Exemples plus complexes :

- Examples Data

N'oubliez pas de citer xarray dans vos publications et présentations.

V-B - Astropy

Astropy est une bibliothèque astronomique maintenue par la communauté et visant à fédérer les efforts jusque-là disparates. Elle offre en outre une interface unifiée à des bibliothèques affiliées plus spécifiques.

V-B-1 - Tour d'horizon

- Structures de base :

- astropy.constants : constantes fondamentales (voir également scipy.constants) ;

- astropy.units : unités et quantités dimensionnées ;

- astropy.nddata : extension des numpy.ndarray (incluant métadonnées, masque, unité, incertitude, etc.) ;

- astropy.table : tableaux hétérogènes ;

- astropy.time : manipulation du temps et des dates ;

- astropy.coordinates : systèmes de coordonnées ;

- astropy.wcs : World Coordinate System ;

- astropy.modeling : modèles et ajustements ;

- astropy.analytic_functions : fonctions analytiques.

- Entrées/sorties :

- astropy.io.fits : fichiers FITS ;

- astropy.io.ascii : tables ASCII ;

- astropy.io.votable : XML VO-tables ;

- astropy.io.misc : divers ;

- astropy.vo : accès au Virtual Observatory.

- Calculs astronomiques :

- astropy.cosmology : calculs cosmologiques ;

- astropy.convolution : convolution et filtrage ;

- astropy.visualization : visualisation de données ;

- astropy.stats : méthodes astrostatistiques.

V-B-2 - Démonstration

Démonstration Astropy (astropy.ipynb)

Nous présentons ici quelques possibilités de la bibliothèque Astropy.

Référence : cette démonstration est très largement inspirée de la partie Astropy du cours Python Euclid 2016.

In [1] :

1. # Mimic Python-3.x

2. from __future__ import print_function, division

4. import numpy as N

5. import matplotlib.pyplot as P

6. try:

7. import seaborn

8. seaborn.set_color_codes()  # Override default matplotlib colors 'b', 'r', 'g', etc.

9. except ImportError:

10. pass

12. # Interactive figures

13. # %matplotlib notebook

14. # Static figures

15. %matplotlib inline

V-B-2-a - Fichiers FITS

Le format FITS (Flexible Image Transport System) constitue le format de données historique (et encore très utilisé) de la communauté astronomique. Il permet le stockage simultané de données - sous forme de tableaux numériques multidimensionnels (spectre 1D, image 2D, cube 3D, etc.) ou de tables de données structurées (texte ou binaires) - et des métadonnées associées - sous la forme d'un entête ASCII nommé header. Il autorise en outre de combiner au sein d'un même fichier différents segments de données (extensions, p.ex. le signal et la variance associée) sous la forme de HDU (Header-Data Units).

Le fichier FITS de test est disponible ici : image.fits (données Herschel Space Observatory).

V-B-2-a-i - Lire un fichier FITS

In [2] :

1. from astropy.io import fits as F

3. filename = "image.fits"

4. hdulist = F.open(filename)

Hdulist est un objet HDUList de type liste regroupant les différents HDU du fichier :

In [3] :

hdulist.info()
Filename: image.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU     151   ()
1    image       ImageHDU        52   (273, 296)   float64
2    error       ImageHDU        20   (273, 296)   float64
3    coverage    ImageHDU        20   (273, 296)   float64
4    History     ImageHDU        23   ()
5    HistoryScript  BinTableHDU     39   105R x 1C    [300A]
6    HistoryTasks  BinTableHDU     46   77R x 4C     [1K, 27A, 1K, 9A]
7    HistoryParameters  BinTableHDU     74   614R x 10C   [1K, 20A, 7A, 46A, 1L, 1K, 1L, 74A, 11A, 41A]

Chaque HDU contient:

- un attribut data pour la partie données sous la forme d'un numpy.array ou d'une structure équivalente à un tableau à type structuré ;

- un attribut header pour la partie métadonnées sous la forme “ KEY = value / comment ”.

In [4] :

imhdu = hdulist['image']
print(type(imhdu.data), type(imhdu.header))

Out [4] :

<type 'numpy.ndarray'> <class 'astropy.io.fits.header.Header'>

Il est également possible de lire directement les données et les métadonnées de l'extension image :

In [5] :

ima, hdr = F.getdata(filename, 'image', header=True)
print(type(ima), type(hdr))

Out [5] :

<type 'numpy.ndarray'> <class 'astropy.io.fits.header.Header'>

data contient donc les données numériques, ici un tableau 2D :

In [6] : N.info(ima)

Out [6] :

class:  ndarray
shape:  (296, 273)
strides:  (2184, 8)
itemsize:  8
aligned:  True
contiguous:  True
fortran:  False
data pointer: 0x7f43bce3bec0
byteorder:  big
byteswap:  True
type: >f8

L'entête hdr est un objet de type Header similaire à un OrderedDict (dictionnaire ordonné).

In [7] : hdr[:5] # Les 5 premières clés de l'entête

Out [7] :

XTENSION= 'IMAGE   '           / Java FITS: Wed Aug 14 11:37:21 CEST 2013
BITPIX  =                  -64
NAXIS   =                    2 / Dimensionality
NAXIS1  =                  273
NAXIS2  =                  296

Les axes des tableaux FITS et NumPy arrays sont inversés !

In [8] :

print("FITS: ", (hdr['naxis1'], hdr['naxis2']))  # format de l'image FITS
print("Numpy:", ima.shape)                       # format du tableau numpy

Out [8] :

FITS:  (273, 296)
Numpy: (296, 273)

V-B-2-a-ii - World Coordinate System

L'entête d'un fichier FITS peut notamment inclure une description détaillée du système de coordonnées lié aux données, le World Coordinate System.

In [9] :

1. from astropy import wcs as WCS

3. wcs = WCS.WCS(hdr)  # Décrypte le WCS à partir de l'entête

4. print(wcs)

Out [9] :

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'
CRVAL : 30.07379502155236  -24.903630299920962
CRPIX : 134.0  153.0
NAXIS    : 273 296

In [10] :

1. ra, dec = wcs.wcs_pix2world(0, 0, 0)  # Coordonnées réelles du px (0, 0)

2. print("World:", ra, dec)

3. x, y = wcs.wcs_world2pix(ra, dec, 0)  # Coordonnées dans l'image de la position (ra, dec)

4. print("Image:", x, y)

Out [10] :

World: 30.318687003 -25.1567606072
Image: -3.21165316564e-12 -1.05160324892e-12

V-B-2-a-iii - Visualisation dans matplotlib

Les tableaux 2D (image) se visualisent à l'aide de la commande imshow :

- cmap : table des couleurs ;

- vmin, vmax : valeurs minimale et maximale des données visualisées ;

- origin : position du pixel (0, 0) (‘lower' = en bas à gauche).

In [11] :

1. fig, ax = P.subplots()

2. ax.imshow(ima, cmap='viridis', origin='lower', interpolation='None', vmin=-2e-2, vmax=5e-2)

3. ax.set_xlabel("x [px]")

4. ax.set_ylabel("y [px]")

5. ax.set_title(filename);

Il est possible d'ajouter d'autres systèmes de coordonnées via WCS.

Requiert WCSaxes, qui n'est pas encore intégré d'office dans astropy.

In [12] :

1. import astropy.visualization as VIZ

2. from astropy.visualization.mpl_normalize import ImageNormalize

4. fig = P.figure(figsize=(8, 8))

5. ax = fig.add_subplot(1, 1, 1, projection=wcs)

7. # 10th and 99th percentiles

8. vmin, vmax = VIZ.AsymmetricPercentileInterval(10, 99).get_limits(ima)

9. # Linear normalization

10. norm = ImageNormalize(vmin=vmin, vmax=vmax, stretch=VIZ.LinearStretch())

12. ax.imshow(ima, cmap='gray', origin='lower', interpolation='None', norm=norm)

14. # Coordonnées équatoriales en rouge

15. ax.coords['ra'].set_axislabel(u'α [J2000]')

16. ax.coords['dec'].set_axislabel(u'δ [J2000]')

18. ax.coords['ra'].set_ticks(color='r')

19. ax.coords['dec'].set_ticks(color='r')

21. ax.coords.grid(color='r', ls='-', lw=2)

23. # Coordonnées galactiques en vert

24. overlay = ax.get_coords_overlay('galactic')

26. overlay['l'].set_axislabel('Galactic Longitude')

27. overlay['b'].set_axislabel('Galactic Latitude')

29. overlay['l'].set_ticks(color='g')

30. overlay['b'].set_ticks(color='g')

32. overlay.grid(color='g', ls='--', lw=2)

V-B-2-a-iv - Tables

Outre les tables FITS, astropy permet lire et écrire des Tables dans de nombreux formats ASCII usuels en astronomie (LaTeX, HTML, CDS, SExtractor, etc.).

Le fichier ASCII de test est disponible ici : sources.dat.

In [13] :

1. from astropy.io import ascii

3. catalog = ascii.read('sources.dat')

4. catalog.info()

Out [13] :

<Table length=167>
name      dtype
------------ -------
ra float64
dec float64
x float64
y float64
raPlusErr float64
decPlusErr float64
raMinusErr float64
decMinusErr float64
xPlusErr float64
yPlusErr float64
xMinusErr float64
yMinusErr float64
flux float64
fluxPlusErr float64
fluxMinusErr float64
background float64
bgPlusErr float64
bgMinusErr float64
quality float64

/usr/local/lib/python2.7/dist-packages/astropy/table/column.py:263: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison return self.data.__eq__(other)

/usr/local/lib/python2.7/dist-packages/astropy/table/info.py:93: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal if np.all(info[name] == ''):

In [14] :

1. catalog.sort('flux')  # Ordonne les sources par 'flux' croissant

2. catalog.reverse()     # Ordonne les sources par 'flux' décroissant

4. #catalog.show_in_notebook(display_length=5)

5. catalog[:5]           # Les 5 sources les plus brillantes du catalogue

Out [14] : <Table length=5>

ra float64

dec
float64

x float64

y float64

raPlusErr float64

30.0736543481

-24.8389847181

133.076596062

190.787365523

0.000137260656107

30.0997563127

-25.030193106

118.886083699

76.0608399428

0.000200560356177

30.2726211379

-25.0771584874

24.9485847511

47.803308157

0.000495906932571

29.8763509609

-24.7518860739

240.583543382

242.969136003

0.000363715086376

29.8668948822

-24.8539846811

245.642862323

181.701895847

0.000391871197269

decPlusErr
float64

RaMinusErr
float64

DecMinusE
rr float64

XplusErr
float64

decPlusErr
float64

0.000124562849667

0.000137260656103

0.000124562850274

0.0747378051241

0.000124562849667

0.00018172429472

0.000200560356813

0.000181724292183

0.109035119289

0.00018172429472

0.000449161858974

0.00049590696235

0.000449161837629

0.269501203678

0.000449161858974

0.000330300804947

0.000363715070463

0.000330300815023

0.198183052816

0.000330300804947

0.000355578581107

0.000391871177911

0.000355578584919

0.213348701017

0.000355578581107

YplusErr
float64

XminusErr
float64

YminusErr
float64

Flux
float64

Flux
PlusErr float64

0.0747378051241

0.0747378051241

0.0747378051241

157.862010775

6.55458899877

0.109035119289

0.109035119289

0.109035119289

118.448614593

7.17503267878

0.269501203678

0.269501203678

0.269501203678

115.452316575

17.2858546025

0.198183052816

0.198183052816

0.198183052816

95.2024034098

10.4819460796

0.213348701017

0.213348701017

0.213348701017

73.0101287468

8.65367562734

Flux
MinusErr float64

Background
float64

BgPlusErr
float64

BgMinusErr
float64

Quality
float64

6.55458899877

-0.00492962840045

0.00280563968109

0.00280563968109

24.0841967062

7.17503267878

0.00226922676457

0.00310958937187

0.00310958937187

16.5084425251

17.2858546025

0.0108994941326

0.00603800958334

0.00603800958334

6.67900541976

10.4819460796

-0.0013197765479

0.00466051306854

0.00466051306854

9.08251222505

8.65367562734

0.0147626682141

0.00330155226713

0.00330155226713

8.43689223988

Histogramme des flux, avec choix automatique du nombre de bins :

In [15] :

1. import astropy.visualization as VIZ

3. fig, ax = P.subplots()

4. VIZ.hist(catalog['flux'], bins='freedman', ax=ax, histtype='stepfilled')

5. ax.set(xlabel="Flux", title="%d sources" % len(catalog));

Après conversion des coordonnées RA-Dec de la table en coordonnées, on peut superposer la position des cinq sources les plus brillantes du catalogue à l'image précédente :

In [16] :

1. fig, ax = P.subplots()

2. ax.imshow(ima, cmap='gray', origin='lower', interpolation='None', vmin=-2e-2, vmax=5e-2)

3. ax.set(xlabel="x [px]", ylabel="y [px]", title=filename)

5. x, y = wcs.wcs_world2pix(catalog['ra'][:5], catalog['dec'][:5], 0)

6. ax.scatter(x, y, color='r', label="Brightest sources")

7. ax.legend();

V-B-2-a-v - Quantités et unités

Astropy permet de définir des Quantités dimensionnées et de gérer les conversions d'unités.

In [17] :

1. from astropy import units as U

2. from astropy import constants as K

4. print("Vitesse de la lumière: {} = {}".format(K.c, K.c.to("Mpc/Ga")))

Out [17] :

Vitesse de la lumière: 299792458.0 m / s = 306.601393788 Mpc / Ga

In [18] :

1. H0 = 70 * U.km / U.s / U.Mpc

2. print ("H0 = {} = {} = {} = {}".format(H0, H0.to("nm/(a*km)"), H0.to("Mpc/(Ga*Gpc)"), H0.decompose()))

Out [18] :

H0 = 70.0 km / (Mpc s) = 71.5898515538 nm / (a km) = 71.5898515538 Mpc / (Ga Gpc) = 2.26854550263e-18 1 / s

In [19] :

1. l = 100 * U.micron

2. print("{} = {}".format(l, l.to(U.GHz, equivalencies=U.spectral())))

4. s = 1 * U.mJy

5. print ("{} = {} à {}".format(s, s.to('erg/(cm**2 * s * angstrom)',

6. equivalencies=U.spectral_density(1 * U.micron)), 1 * U.micron))

Out [19] :

100.0 micron = 2997.92458 GHz
1.0 mJy = 2.99792458e-16 erg / (Angstrom cm2 s) à 1.0 micron

V-B-2-a-vi - Calculs cosmologiques

Astropy permet des calculs de base de cosmologie.

In [20] :

from astropy.cosmology import Planck15 as cosmo
print(cosmo)

Out [20] :

FlatLambdaCDM(name="Planck15", H0=67.7 km / (Mpc s), Om0=0.307, Tcmb0=2.725 K, Neff=3.05, m_nu=[ 0.    0.    0.06] eV, Ob0=0.0486)

In [21] :

1. z = N.logspace(-2, 1, 100)

3. fig = P.figure(figsize=(14, 6))

5. # Ax1: lookback time

6. ax1 = fig.add_subplot(1, 2, 1,

7. xlabel="Redshift", xscale='log')

8. ax1.plot(z, cosmo.lookback_time(z), 'b-')

9. ax1.set_ylabel("Lookback time [Ga]", color='b')

10. ax1.set_yscale('log')

12. ax1.xaxis.set_minor_locator(P.matplotlib.ticker.LogLocator(subs=range(2,10)))

13. ax1.yaxis.set_minor_locator(P.matplotlib.ticker.LogLocator(subs=range(2,10)))

14. ax1.grid(which='minor', color='w', linewidth=0.5)

16. # En parallèle : facteur d'échelle

17. ax1b = ax1.twinx()

18. ax1b.plot(z, cosmo.scale_factor(z), 'r-')

19. ax1b.set_ylabel("Scale factor", color='r')

20. ax1b.set_yscale('log')

21. ax1b.grid(False)

23. ht = (1/cosmo.H0).to('Ga')  # Hubble time

24. ax1.axhline(ht.value, c='b', ls='--')

25. ax1.annotate("Hubble time = {:.2f}".format(ht), (z[1], ht.value), (3,3),

26. textcoords='offset points', size='small');

28. # Ax2: distances de luminosité et angulaire

29. ax2 = fig.add_subplot(1, 2, 2,

30. xlabel="Redshift", xscale='log')

31. ax2.plot(z, cosmo.luminosity_distance(z), 'k-', label="Luminosity distance")

32. ax2.plot(z, cosmo.angular_diameter_distance(z), 'k--', label="Angular distance")

33. ax2.set_ylabel("Distance [Mpc]")

34. ax2.set_yscale('log')

35. ax2.legend()

37. ax2.xaxis.set_minor_locator(P.matplotlib.ticker.LogLocator(subs=range(2,10)))

38. ax2.yaxis.set_minor_locator(P.matplotlib.ticker.LogLocator(subs=range(2,10)))

39. ax2.grid(which='minor', color='w', linewidth=0.5)

41. # En parallèle, module de distance

42. ax2b = ax2.twinx()

43. ax2b.plot(z, cosmo.distmod(z), 'k-', visible=False)  # Just to get the scale

44. ax2b.set_ylabel("Distance modulus")

46. fig.subplots_adjust(wspace=0.3)

47. fig.suptitle(cosmo.name, fontsize='x-large');

Voir également :

- AstroBetter tutorials

N'oubliez pas de citer [Astropy13] ou de mentionner l'utilisation d'astropy dans vos publications et présentations.

V-C - Autres bibliothèques scientifiques

Python est maintenant très largement utilisé par la communauté scientifique, et dispose d'innombrables bibliothèques dédiées aux différents domaines de la physique, chimie, etc. :

- Astronomie : Kapteyn, AstroML, HyperSpy ;

- Mécanique quantique : QuTiP ;

- Électromagnétisme : EMpy ;

- PDE solver : FiPy, SfePy ;

- Analyse statistique bayésienne : PyStan ;

- Markov Chain Monte-Carlo : emcee, PyMC3 ;

- Machine Learning : mlpy, milk, MDP, et autres modules d'intelligence artificielle ;

- Calcul symbolique : sympy (voir également ce tutoriel sympy) et sage ;

- PyROOT ;

- High Performance Computing in Python ;

- Etc.

VI - Développer en Python

VI-A - Le zen du Python

Le zen du Python (PEP 20) est une série de 20 aphorismes(7) donnant les grands principes du Python :

>>> import this

1. Beautiful is better than ugly.

2. Explicit is better than implicit.

3. Simple is better than complex.

4. Complex is better than complicated.

5. Flat is better than nested.

6. Sparse is better than dense.

7. Readability counts.

8. Special cases aren't special enough to break the rules.

9. Although practicality beats purity.

10. Errors should never pass silently.

11. Unless explicitly silenced.

12. In the face of ambiguity, refuse the temptation to guess.

13. There should be one- and preferably only one -obvious way to do it.

14. Although that way may not be obvious at first unless you're Dutch.

15. Now is better than never.

16. Although never is often better than right now.

17. If the implementation is hard to explain, it's a bad idea.

18. If the implementation is easy to explain, it may be a good idea.

19. Namespaces are one honking great idea - let's do more of those!

Une traduction libre en français :

1. Préférer le beau au laid.

2. L'explicite à l'implicite.

3. Le simple au complexe.

4. Le complexe au compliqué.

5. Le déroulé à l'imbriqué.

6. L'aéré au compact.

7. La lisibilité compte.

8. Les cas particuliers ne le sont jamais assez pour violer les règles.

9. Même s'il faut privilégier l'aspect pratique à la pureté.

10. Ne jamais passer les erreurs sous silence.

11. Ou les faire taire explicitement.

12. Face à l'ambiguïté, ne pas se laisser tenter à deviner.

13. Il doit y avoir une - et si possible une seule - façon évidente de procéder.

14. Même si cette façon n'est pas évidente à première vue, à moins d'être Hollandais.

15. Mieux vaut maintenant que jamais.

16. Même si jamais est souvent mieux qu'immédiatement.

17. Si l'implémentation s'explique difficilement, c'est une mauvaise idée.

18. Si l'implémentation s'explique facilement, c'est peut-être une bonne idée.

19. Les espaces de nommage sont une sacrée bonne idée, utilisons-les plus souvent !

VI-A-1 - Us et coutumes

- Fail early, fail often, fail better! (raise)

- Easier to Ask for Forgiveness than Permission (try ... except)

- le Style Guide for Python Code (PEP 8)

- Idioms and Anti-Idioms in Python

- Code Like a Pythonista: Idiomatic Python

- Google Python Style Guide

- The Best of the Best Practices (BOBP) Guide for Python

Quelques conseils supplémentaires

- « Ne réinventez pas la roue, sauf si vous souhaitez en savoir plus sur les roues » (Jeff Atwood(8)) : cherchez si ce que vous voulez faire n'a pas déjà été fait (éventuellement en mieux…) pour vous concentrer sur votre valeur ajoutée, réutilisez le code (en citant évidemment vos sources), améliorez-le, et contribuez en retour si possible !

- Écrivez des programmes pour les humains, pas pour les ordinateurs : codez proprement, structurez vos algorithmes, commentez votre code, utilisez des noms de variable qui ont un sens, soignez le style et le formatage, etc.

- Codez proprement dès le début : ne croyez pas que vous ne relirez jamais votre code (ou même que personne n'aura jamais à le lire), ou que vous aurez le temps de le refaire mieux plus tard…

- « L'optimisation prématurée est la source de tous les maux » (Donald Knuth(9)) : mieux vaut un code lent mais juste et maintenable qu'un code rapide et faux ou incompréhensible. Dans l'ordre absolu des priorités :

1. Make it work ;

2. Make it right ;

3. Make it fast.

- Respectez le zen du python, il vous le rendra.

VI-A-2 - Principes de conception logicielle

La bonne conception d'un programme va permettre de gérer efficacement la complexité des algorithmes, de faciliter la maintenance (p.ex. correction des erreurs) et d'accroître les possibilités d'extension.

Modularité

Le code est structuré en répertoires, fichiers, classes, méthodes et fonctions. Les blocs ne font pas plus de quelques dizaines de lignes, les fonctions ne prennent que quelques arguments, la structure logique n'est pas trop complexe, etc.

En particulier, le code doit respecter le principe de responsabilité unique : chaque entité élémentaire (classe, méthode, fonction) ne doit avoir qu'une unique raison d'exister, et ne pas tenter d'effectuer plusieurs tâches sans rapport direct (p.ex. lecture d'un fichier de données et analyse des données).

Flexibilité

Une modification du comportement du code (p.ex. l'ajout d'une nouvelle fonctionnalité) ne nécessite de changer le code qu'en un nombre restreint de points.

Un code rigide devient rapidement difficile à faire évoluer, puisque chaque changement requiert un grand nombre de modifications.

Robustesse

La modification du code en un point ne change pas de façon inopinée le comportement dans une autre partie a priori non reliée.

Un code fragile est facile à modifier, mais chaque modification peut avoir des conséquences inattendues et le code tend à devenir instable.

Réutilisabilité

La réutilisation d'une portion de code ne demande pas de changement majeur, n'introduit pas trop de dépendances, et ne conduit pas à une duplication du code.

L'application de ces principes de développement dépend évidemment de l'objectif final du code :

- une bibliothèque centrale (utilisée par de nombreux programmes) favorisera la robustesse et la réutilisabilité aux dépens de la flexibilité : elle devra être particulièrement bien pensée, et ne pourra être modifiée qu'avec parcimonie ;

- inversement, un script d'analyse de haut niveau, d'utilisation restreinte, pourra être plus flexible, mais plus fragile et peu réutilisable.

VI-B - Développement piloté par les tests

Le Test Driven Development (TDD, ou en français « développement piloté par les tests ») est une méthode de programmation qui permet d'éviter des bogues a priori plutôt que de les résoudre a posteriori. Ce n'est pas une méthode propre à Python, elle est utilisée très largement par les programmeurs professionnels.

Le cycle préconisé par TDD comporte cinq étapes :

1. Écrire un premier test ;

2. Vérifier qu'il échoue (puisque le code qu'il teste n'existe pas encore), afin de s'assurer que le test est valide et exécuté ;

3. Écrire un code minimal pour passer le test ;

4. Vérifier que le test passe correctement ;

5. Éventuellement, « réusiner » le code (refactoring), c'est-à-dire l'améliorer (rapidité, lisibilité) tout en gardant les mêmes fonctionnalités.

Diviser pour mieux régner : chaque fonction, classe ou méthode est testée indépendamment. Ainsi, lorsqu'un nouveau morceau de code ne passe pas les tests qui y sont associés, il est certain que l'erreur provient de cette nouvelle partie et non des fonctions ou objets que ce morceau de code utilise. On distingue ainsi hiérarchiquement :

1. Les tests unitaires vérifient individuellement chacune des fonctions, méthodes, etc. ;

2. Les tests d'intégration évaluent les interactions entre différentes unités du programme ;

3. Les tests système assurent le bon fonctionnement du programme dans sa globalité.

Il est très utile de transformer toutes les vérifications réalisées au cours du développement et du débogage sous forme de tests, ce qui permet de les réutiliser lorsque l'on veut compléter ou améliorer une partie du code. Si le nouveau code passe toujours les anciens tests, on est alors sûr de ne pas avoir cassé les fonctionnalités précédentes (régressions).

Nous avons déjà vu aux TD précédents plusieurs façons de rédiger des tests unitaires :

- Les doctest sont des exemples (assez simples) d'exécution de code inclus dans les docstring des classes ou fonctions :

1. def mean_power(alist, power=1):

2. """

3. Retourne la racine `power` de la moyenne des éléments de `alist` à

4. la puissance `power`:

6. .. math:: \mu = (\frac{1}{N}\sum_{i=0}^{N-1} x_i^p)^{1/p}

8. `power=1` correspond à la moyenne arithmétique, `power=2` au *Root

9. Mean Squared*, etc.

11. Exemples:

12. >>> mean_power([1, 2, 3])

13. 2.0

14. >>> mean_power([1, 2, 3], power=2)

15. 2.160246899469287

16. """

18. s = 0.                  # Initialisation de la variable *s* comme *float*

19. for val in alist:       # Boucle sur les éléments de *alist*

20. s += val ** power   # *s* est augmenté de *val* puissance *power*

21. # *mean* = (somme valeurs / nb valeurs)**(1/power)

22. mean = (s / len(alist)) ** (1 / power)  # ATTENTION aux divisions euclidiennes !

24. return mean

Les doctests peuvent être exécutés de différentes façons (voir ci-dessous) :

- avec le module standard doctest : python -m doctest -v mean_power.py ;

- avec pytest : py.test --doctest-modules -v mean_power.py ;

- avec nose : nosetests --with-doctest -v mean_power.py ;

- les fonctions dont le nom commence par test_ et contenant des assert sont automatiquement détectées par pytest(10). Cette méthode permet d'effectuer des tests plus poussés que les doctests, éventuellement dans un fichier séparé du code à tester. P.ex. :

1. def test_empty_init():

2. with pytest.raises(TypeError):

3. youki = Animal()

6. def test_wrong_init():

7. with pytest.raises(ValueError):

8. youki = Animal('Youki', 'lalala')

11. def test_init():

12. youki = Animal('Youki', 600)

13. assert youki.masse == 600

14. assert youki.vivant

15. assert youki.estVivant()

16. assert not youki.empoisonne

Les tests sont exécutés via py.test programme.py ;

- le module unittest de la bibliothèque standard permet à peu près la même chose que pytest, mais avec une syntaxe souvent plus lourde. unittest est étendu par le module non standard nose.

VI-C - Outils de développement

Je fournis ici essentiellement des liens vers des outils pouvant être utiles pour développer en Python.

VI-C-1 - Integrated Development Environment

- idle, l'IDE intégré à Python.

- emacs + python-mode pour l'édition, et ipython pour l'exécution de code (voir Python Programming In Emacs).

- spyder.

- PythonToolkit.

- pyCharm (la version community est gratuite).

- Etc.

VI-C-2 - Vérification du code

Il s'agit d'outils permettant de vérifier a priori la validité stylistique et syntaxique du code, de mettre en évidence des constructions dangereuses, les variables non définies, etc. Ces outils ne testent pas nécessairement la validité des algorithmes et de leur mise en œuvre…

- pycodestyle (ex-pep8) et autopep8.

- pyflakes.

- pychecker.

- pylint.

VI-C-3 - Débogage

Les débogueurs permettent de se « plonger » dans un code en cours d'exécution ou juste après une erreur (analyse postmortem).

- Module de la bibliothèque standard : pdb
Pour déboguer un script, il est possible de l'exécuter sous le contrôle du débogueur pdb en s'interrompant dès la 1re instruction :

python -m pdb script.py
(Pdb)

Commandes (très similaires à gdb) :

- h[elp] [command] : aide en ligne ;

- q[uit] : quitter ;

- r[un] [args] : exécuter le programme avec les arguments ;

- d[own]/u[p] : monter/descendre dans le stack (empilement des appels de fonction) ;

- p expression : afficher le résultat de l'expression (pp : pretty-print) ;

- l[ist] [first[,last]] : afficher le code source autour de l'instruction courante (ll : long list) ;

- n[ext]/s[tep] : exécuter l'instruction suivante (sans y entrer/en y entrant) ;

- unt[il] : continuer l'exécution jusqu'à la ligne suivante (utile pour les boucles) ;

- c[ont[inue]] : continuer l'exécution (jusqu'à la prochaine interruption ou la fin du programme) ;

- r[eturn] : continuer l'exécution jusqu'à la sortie de la fonction ;

- b[reak] [[filename:]lineno | function[, condition]] : mettre en place un point d'arrêt (tbreak pour un point d'arrêt temporaire). Sans argument, afficher les points d'arrêts déjà définis ;

- disable/enable [bpnumber] : désactiver/réactiver tous ou un point d'arrêt ;

- cl[ear] [bpnumber] : éliminer tous ou un point d'arrêt ;

- ignore bpnumber [count] : ignorer un point d'arrêt une ou plusieurs fois ;

- condition bpnumber : ajouter une condition à un point d'arrêt ;

- commands [bpnumber] : ajouter des instructions à un point d'arrêt ;

- commandes ipython : %run monScript.py, %debug, %pdb
Si un script exécuté sous ipython (commande %run) génère une exception, il est possible d'inspecter l'état de la mémoire au moment de l'erreur avec la commande %debug, qui lance une session pdb au point d'arrêt. %pdb on lance systématiquement le débogueur à chaque exception.

L'activité de débogage s'intègre naturellement à la nécessité d'écrire des tests unitaires :

1. Trouver un bogue ;

2. Écrire un test qui aurait dû être validé en l'absence du bogue ;

3. Corriger le code jusqu'à validation du test.

Vous aurez alors finalement corrigé le bogue, et écrit un test s'assurant que ce bogue ne réapparaîtra pas inopinément.

VI-C-4 - Profilage et optimisation

Premature optimization is the root of all evil - Donald Knuth

Avant toute optimisation, s'assurer extensivement que le code fonctionne et produit les bons résultats dans tous les cas. S'il reste trop lent ou gourmand en mémoire pour vos besoins, il peut être nécessaire de l'optimiser.

Le profilage permet de déterminer le temps passé dans chacune des sous-fonctions d'un code (ou ligne par ligne : line profiler, ou selon l'utilisation de la mémoire : memory profiler), afin d'y identifier les parties qui gagneront à être optimisées.

- python -O, __debug__, assert
Il existe un mode « optimisé » de python (option -O), qui pour l'instant ne fait pas grand-chose (et n'est donc guère utilisé…) :

- la variable interne __debug__ passe de True à False ;

- les instructions assert ne sont pas évaluées.

- timeit et %timeit statement sous ipython :

1. In [1]: def t1(n):

2. ...:     l = []

3. ...:     for i in range(n):

4. ...:         l.append(i**2)

5. ...:     return l

6. ...:

7. ...: def t2(n):

8. ...:     return [ i**2 for i in xrange(n) ]

9. ...:

10. ...: def t3(n):

11. ...:     return N.arange(n)**2

12. ...:

13. In [2]: %timeit t1(10000)

14. 1000 loops, best of 3: 950 µs per loop

15. In [3]: %timeit t2(10000)

16. 1000 loops, best of 3: 599 µs per loop

17. In [4]: %timeit t3(10000)

18. 10000 loops, best of 3: 18.1 µs per loop

- cProfile et pstats, et %prun statement sous ipython :

$ python -m cProfile -o output.pstats monScript.py
$ python -m pstats output.pstats

- Tutoriel de profilage

Une fois identifiée la partie du code à optimiser, quelques conseils généraux :

- en cas de doute, favoriser la lisibilité aux performances ;

- utiliser des opérations sur les tableaux, plutôt que sur des éléments individuels (vectorization) : listes en compréhension, tableaux numpy (qui ont eux-mêmes été optimisés) ;

- cython est un langage de programmation compilé très similaire à Python. Il permet d'écrire des extensions en C avec la facilité de Python (voir notamment Working with Numpy) ;

- numba permet automagiquement de compiler à la volée (JIT) du pur code Python via le compilateur LLVM, avec une optimisation selon le CPU (éventuellement le GPU) utilisé, p.ex. :

1. from numba import guvectorize

3. @guvectorize(['void(float64[:], intp[:], float64[:])'], '(n),()->(n)')

4. def move_mean(a, window_arr, out):

5. ...

- à l'avenir, l'interpréteur CPython actuel sera éventuellement remplacé par pypy, basé sur une compilation JIT.

Lien : Performance tips

VI-C-5 - Documentation

- Outils de documentation, ou comment transformer automagiquement un code source bien documenté en une documentation fonctionnelle :

- Sphinx ;

- reStructuredText for Sphinx ;

- Awesome Sphinx ;

- apidoc (documentation automatique).

- Conventions de documentation :

- Docstring convention : PEP 257 ;

- Documenting Your Project Using Sphinx ;

- A Guide to NumPy/SciPy Documentation ;

- Sample doc (matplotlib).

Lien :

Documentation Tools

VI-C-6 - Python packages

Comment installer/créer des modules externes :

- pip ;

- Hitchhiker's Guide to Packaging ;

- Packaging a python library ;

- Cookiecutter est un générateur de squelettes de projet via des templates (pas uniquement pour Python) ;

- cx-freeze <http://cx-freeze.readthedocs.io/>, pour générer un exécutable à partir d'un script.

VI-C-7 - Système de gestion de versions

La gestion des versions du code permet de suivre avec précision l'historique des modifications du code (ou de tout autre projet), de retrouver les changements critiques, de développer des branches alternatives, de faciliter le travail collaboratif, etc.

Git est un VCS particulièrement performant (p.ex. utilisé pour le développement du noyau Linux(11)). Il est souvent couplé à un dépôt en ligne faisant office de dépôt de référence et de solution de sauvegarde, et offrant généralement des solutions d'intégration continue, p.ex. :

- les très célèbres GitHub et GitLab, gratuits pour les projets libres ;

- pour des projets liés à votre travail, je conseille plutôt des dépôts directement gérés par votre institution, p.ex. GitLab-IN2P3.

Git mérite un cours en soi, et devrait être utilisé très largement pour l'ensemble de vos projets (p.ex. rédaction d'articles, de thèse de cours, fichiers de configuration, tests numériques, etc.).

Quelques liens d'introduction :

- Pro-Git book, le livre « officiel » ;

- Git Immersion ;

- Git Magic.

VI-C-8 - Intégration continue

L'intégration continue est un ensemble de pratiques de développement logiciel visant à s'assurer de façon systématique que chaque modification du code n'induit aucune régression, et passe l'ensemble des tests. Cela passe généralement par la mise en place d'un système de gestion des sources, auquel est accolé un mécanisme automatique de compilation (build), de déploiement sur les différentes infrastructures, d'exécution des tests (unitaires, intégration, fonctionnels, etc.) et de mise à disposition des résultats, de mise en ligne de la documentation, etc.

La plupart des développements des logiciels open source majeurs se fait maintenant sous intégration continue en utilisant des services en ligne directement connectés au dépôt source. Exemple sur Astropy :

- Travis CI intégration continue ;

- Coveralls taux de couverture des tests unitaires ;

- Readthedocs documentation en ligne ;

- Depsy mise en valeur du développement logiciel dans le monde académique (measure the value of software that powers science).

VI-D - Python 2 vs. Python 3

Il existe de nombreux outils permettant de faciliter la transition 2.x → 3.x.

- L'interpréteur Python 2.7 dispose d'une option -3 mettant en évidence dans un code les parties qui devront être modifiées pour un passage à Python 3.

- Le script 2to3 permet également d'automatiser la conversion du code 2.x en 3.x.

- La bibliothèque standard __future__ permet d'utiliser nativement des constructions 3.x dans un code 2.x, p.ex. :

1. from __future__ import print_function  # Fonction print()

2. from __future__ import division        # Division non euclidienne

4. print(1/2)                             # Affichera '0.5'

- La bibliothèque non standard six fournit une couche de compatibilité 2.x-3.x, permettant de produire de façon transparente un code compatible simultanément avec les deux versions.

Liens

- Python 2 or python 3 ?

- Porting Python 2 Code to Python 3

- Python 2/3 compatibility

VII - Références supplémentaires

Voici une liste très partielle de documents Python disponibles en ligne. La majorité des liens sont en anglais, quelques-uns sont en français.

VII-A - Documentation générale

- Python

- Python Wiki

- Python Frequently Asked Questions

- Python 2.7 Quick Reference

- The Python Package Index

VII-B - Listes de liens

- Python facile (2005)

- Python : quelques références, trucs et astuces : (2014)

- Improving your programming style in Python (2014)

- Starter Kit (py4science) (2010)

- Learning Python For Data Science (2016)

- Awesome Python

- A curated list of courses on Python

ipython

- IPython tutorial

- IPython cookbook

- IPython en ligne

- IPython quick refsheet

Expressions rationnelles

- regex tester

Python 3.x

- 10 awesome features of Python that you can't use because you refuse to upgrade to Python 3

VII-C - Livres libres

- How to Think Like a Computer Scientist

- Wikibook

- Interactive edition

- Dive into Python

- 10 Free Python Programming Books

- A Python Book

- Start Programming with Python

- Learn Python the Hard Way

- Python for Informatics: Exploring Information

- Intermediate Python

- Apprendre à programmer avec Python

- Programmation Python

VII-D - Cours en ligne

VII-D-1 - Python

- Python Tutorial (v2.7)

- Tutoriel Python (v2.4)

- Apprenez à programmer en Python

- Débuter avec Python au lycée

- Présentation de Python

- Introduction à Python pour la programmation scientifique

- Begginer's Guide

- DIY python workshop

- Google's Python Class

- CheckIO, pour apprendre la programmation Python en s'amusant !

- Python Testing Tools

- Python Programming (Code Academy)

VII-D-2 - Scientifique

- Python Scientific Lecture Notes

- Handbook of the Physics Computing Course (2002)

- Practical Scientific Computing in Python

- Computational Physics with Python (avec exercices)

- SciPy tutorials (numpy, scipy, matplotlib, ipython) : 2011, 2012, 2013,

- Scipy Central (code snippets)

- Advance Scientific Programming in Python

- Lectures on Computational Economics (avec exercices)

- Intro to Python for Data Science (DataCamp avec vidéos & exercices)

- Python for Data Science

- Learning Python For Data Science

- Computational Statistics in Python

- Python Data Science Handbook

En français

- Formation à Python scientifique

- NumPy et SciPy

- L'informatique scientifique avec Python

- La programmation scientifique avec Python

Astronomie

- Practical Python for Astronomers

- Astropy tutorials

- Python for Euclid 2016

- Advanced software programming for astrophysics and astroparticle physics (1st ASTERICS-OBELICS International School, voir Timetable/Vue détaillée)

VII-D-3 - Snippets

- Python cheatsheets

- Scipy Central (code snippets)

VIII - Exemples

VIII-A - Mean power (fonction, argparse)

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. """

5. Exemple de script (shebang, docstring, etc.) permettant une

6. utilisation en module (`import mean_power`) et en exécutable (`python

7. mean_power.py -h`);

8. """

10. from __future__ import division  # Les divisions entre entiers ne sont pas euclidiennes

12. def mean_power(alist, power=1):

13. """

14. Retourne la racine `power` de la moyenne des éléments de `alist` à

15. la puissance `power`:

17. .. math:: \mu = (\frac{1}{N}\sum_{i=0}^{N-1} x_i^p)^{1/p}

19. `power=1` correspond à la moyenne arithmétique, `power=2` au *Root

20. Mean Squared*, etc.

22. Exemples:

23. >>> mean_power([1, 2, 3])

24. 2.0

25. >>> mean_power([1, 2, 3], power=2)

26. 2.160246899469287

27. """

29. s = 0.                  # Initialisation de la variable *s* comme *float*

30. for val in alist:       # Boucle sur les éléments de *alist*

31. s += val ** power   # *s* est augmenté de *val* puissance *power*

32. # *mean* = (somme valeurs / nb valeurs)**(1/power)

33. mean = (s / len(alist)) ** (1 / power)  # ATTENTION aux divisions euclidiennes !

35. return mean

38. if __name__ == '__main__':

40. # start-argparse

41. import argparse

43. parser = argparse.ArgumentParser()

44. parser.add_argument('list', nargs='*', type=float, metavar='nombres',

45. help="Liste de nombres à moyenner")

46. parser.add_argument('-i', '--input', nargs='?', type=file,

47. help="Fichier contenant les nombres à moyenner")

48. parser.add_argument('-p', '--power', type=float, default=1.,

49. help="'Puissance' de la moyenne (%default)")

51. args = parser.parse_args()

52. # end-argparse

54. if args.input:       # Lecture des coordonnées du fichier d'entrée

55. # Le fichier a déjà été ouvert en lecture par argparse (type=file)

56. try:

57. args.list = [float(x) for x in args.input

58. if not x.strip().startswith('#')]

59. except ValueError:

60. parser.error(

61. "Impossible de déchiffrer la ligne '{}' du fichier '{}'".format(

62. x, args.input))

64. # Vérifie qu'il y a au moins un nombre dans la liste

65. if not args.list:

66. parser.error("La liste doit contenir au moins un nombre")

68. # Calcul

69. moyenne = mean_power(alist, args.power)

71. # Affichage du résultat

72. print "La moyenne des {} nombres à la puissance {} est {}".format(

73. len(alist), args.power, moyenne)

Source : mean_power.py

VIII-B - Formes (POO)

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. """

5. Exemple de Programmation Orientée Objet.

6. """

8. __author__ = "Mathieu Leocmach <mathieu.leocmach@ens-lyon.fr>"

9. __version__ = "Time-stamp: <2014-10-03 10:54 mathieu.leocmach@ens-lyon.fr>"

12. # Définition d'une classe ==============================

14. class Forme(object):  # *object* est la classe dont dérivent toutes les autres

16. """Une forme plane, avec éventuellement une couleur."""

18. def __init__(self, couleur=None):

19. """Initialisation d'une Forme, sans couleur par défaut."""

21. if couleur is None:

22. self.couleur = 'indéfinie'

23. else:

24. self.couleur = couleur

26. def __str__(self):

27. """

28. Surcharge de la fonction `str()`: l'affichage *informel* de

29. l'objet dans l'interpréteur, p.ex. `print a` sera résolu comme

30. `a.__str__()`

32. Retourne une chaîne de caractères.

33. """

35. return "Forme encore indéfinie de couleur {}".format(self.couleur)

37. def change_couleur(self, newcolor):

38. """Change la couleur de la Forme."""

40. self.couleur = newcolor

42. def aire(self):

43. """

44. Renvoi l'aire de la Forme.

46. L'aire ne peut pas être calculée dans le cas où la forme n'est

47. pas encore spécifiée : c'est ce que l'on appelle une méthode

48. 'abstraite', qui pourra être précisée dans les classes filles.

49. """

51. raise NotImplementedError(

52. "Impossible de calculer l'aire d'une forme indéfinie.")

54. def __cmp__(self, other):

55. """

56. Comparaison de deux Formes sur la base de leur aire.

58. Surcharge des opérateurs de comparaison de type `{self} <

59. {other}`: la comparaison sera résolue comme

60. `self.__cmp__(other)` et le résultat sera correctement

61. interprété.

63. .. WARNING:: cette construction n'est plus supportée en Python3.

64. """

66. return cmp(self.aire(), other.aire())  # Opérateur de comparaison

69. class Rectangle(Forme):

71. """

72. Un Rectangle est une Forme particulière.

74. La classe fille hérite des attributs et méthodes de la

75. classe -mère, mais peut les surcharger (i.e. en changer la

76. définition), ou en ajouter de nouveaux:

78. - les méthodes `Rectangle.change_couleur()` et

79. `Rectangle.__cmp__()` dérivent directement de

80. `Forme.change_couleur()` et `Forme.__cmp__()`;

81. - `Rectangle.__str__()` surcharge `Forme.__str__()`;

82. - `Rectangle.aire()` définit la méthode jusqu'alors abstraite

83. `Forme.aire()`;

84. - `Rectangle.allonger()` est une nouvelle méthode propre à

85. `Rectangle`.

86. """

88. def __init__(self, longueur, largeur, couleur=None):

89. """

90. Initialisation d'un Rectangle longueur × largeur, sans couleur par

91. défaut.

92. """

94. # Initialisation de la classe parente (nécessaire pour assurer

95. # l'héritage)

96. Forme.__init__(self, couleur)

98. # Attributs propres à la classe Rectangle

99. self.longueur = longueur

100. self.largeur = largeur

102. def __str__(self):

103. """Surcharge de `Forme.__str__()`."""

105. return "Rectangle {}x{}, de couleur {}".format(

106. self.longueur, self.largeur, self.couleur)

108. def aire(self):

109. """

110. Renvoie l'aire du Rectangle.

112. Cette méthode définit la méthode abstraite `Forme.area()`,

113. pour les Rectangles uniquement.

114. """

116. return self.longueur * self.largeur

118. def allonger(self, facteur):

119. """Multiplie la *longueur* du Rectangle par un facteur"""

121. self.longueur *= facteur

124. if __name__ == '__main__':

126. s = Forme()                       # Forme indéfinie et sans couleur

127. print "s:", str(s)                # Interprété comme `s.__str__()`

128. s.change_couleur('rouge')         # On change la couleur

129. print "s après change_couleur:", str(s)

130. try:

131. print "Aire de s:", s.aire()  # La méthode abstraite lève une exception

132. except NotImplementedError as err:

133. print err

135. q = Rectangle(1, 4, 'vert')       # Rectangle 1×4 vert

136. print "q:", str(q)

137. print "Aire de q:", q.aire()

139. r = Rectangle(2, 1, 'bleu')       # Rectangle 2×1 bleu

140. print "r:", str(r)

141. print "Aire de r:", r.aire()

143. print "r >= q:", (r >= q)         # Interprété comme r.__cmp__(q)

145. r.allonger(2)                     # r devient un rectangle 4×1

146. print "Aire de r après l'avoir allongé d'un facteur 2:", r.aire()

147. print "r >= q:", (r >= q)

Source : formes.py

VIII-C - Cercle circonscrit (POO, argparse)

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. """

5. Calcule le cercle circonscrit à 3 points du plan.

7. Ce script sert d'illustration à plusieurs concepts indépendants:

9. - un exemple de script (shebang, docstring, etc.) permettant une

10. utilisation en module (`import circonscrit`) et en exécutable

11. (`python circonscrit.py -h`);

12. - des exemples de Programmation Orientée Objet : classe `Point` et la

13. classe héritière `Vector`;

14. - un exemple d'utilisation du module `argparse` de la bibliothèque

15. standard, permettant la gestion des arguments de la ligne de

16. commande ;

17. - l'utilisation de tests unitaires sous la forme de `doctest` (tests

18. inclus dans les *docstrings* des éléments à tester).

20. Pour exécuter les tests unitaires du module :

22. - avec doctest: `python -m doctest -v circonscrit.py`

23. - avec pytests: `py.test --doctest-modules -v circonscrit.py`

24. - avec nose:    `nosetests --with-doctest -v circonscrit.py`

25. """

27. __author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"

28. __version__ = "Time-stamp: <2014-01-12 22:19 ycopin@lyonovae03.in2p3.fr>"

30. # Définition d'une classe ==============================

33. class Point(object):  # *object* est la classe dont dérivent toutes les autres

35. """

36. Classe définissant un `Point` du plan, caractérisé par ses

37. coordonnées `x`,`y`.

38. """

40. def __init__(self, x, y):

41. """

42. Méthode d'instanciation à partir de deux coordonnées réelles.

44. >>> Point(0,1)          # doctest: +ELLIPSIS

45. <circonscrit.Point object at 0x...>

46. >>> Point(1+3j)

47. Traceback (most recent call last):

48. ...

49. TypeError: __init__() takes exactly 3 arguments (2 given)

50. """

52. try:  # Convertit les coords en `float`

53. self.x = float(x)

54. self.y = float(y)

55. except (ValueError, TypeError):

56. raise TypeError("Invalid input coordinates ({},{})".format(x, y))

58. def __str__(self):

59. """

60. Surcharge de la fonction `str()`: l'affichage *informel* de l'objet

61. dans l'interpréteur, p.ex. `print self` sera résolu comme

62. `self.__str__()`

64. Retourne une chaîne de caractères.

66. >>> print Point(1,2)

67. Point (x=1.0, y=2.0)

68. """

70. return "Point (x={p.x}, y={p.y})".format(p=self)

72. def isOrigin(self):

73. """

74. Teste si le point est à l'origine en testant la nullité des deux

75. coordonnées.

77. Attention aux éventuelles erreurs d'arrondis: il faut tester

78. la nullité à la précision numérique près.

80. >>> Point(1,2).isOrigin()

81. False

82. >>> Point(0,0).isOrigin()

83. True

84. """

86. import sys

88. eps = sys.float_info.epsilon  # Le plus petit float non nul

90. return ((abs(self.x) <= eps) and (abs(self.y) <= eps))

92. def distance(self, other):

93. """

94. Méthode de calcul de la distance du point (`self`) à un autre point

95. (`other`).

97. >>> A = Point(1,0); B = Point(1,1); A.distance(B)

98. 1.0

99. """

101. from math import hypot

103. return hypot(self.x - other.x, self.y - other.y)  # sqrt(dx**2 + dy**2)

106. # Définition du point origine O

107. O = Point(0, 0)

110. # Héritage de classe ==============================

113. class Vector(Point):

115. """

116. Un `Vector` hérite de `Point` avec des méthodes additionnelles

117. (p.ex. la négation d'un vecteur, l'addition de deux vecteurs, ou

118. la rotation d'un vecteur).

119. """

121. def __init__(self, A, B):

122. """

123. Définit le vecteur `AB` à partir des deux points `A` et `B`.

125. >>> Vector(Point(1,0), Point(1,1)) # doctest: +ELLIPSIS

126. <circonscrit.Vector object at 0x...>

127. >>> Vector(0, 1)

128. Traceback (most recent call last):

129. ...

130. AttributeError: 'int' object has no attribute 'x'

131. """

133. # Initialisation de la classe parente

134. Point.__init__(self, B.x - A.x, B.y - A.y)

136. # Attribut propre à la classe dérivée

137. self.sqnorm = self.x ** 2 + self.y ** 2  # Norme du vecteur au carré

139. def __str__(self):

140. """

141. Surcharge de la fonction `str()`: `print self` sera résolu comme

142. `Vector.__str__(self)` (et non pas comme

143. `Point.__str__(self)`)

145. >>> A = Point(1,0); B = Point(1,1); print Vector(A,B)

146. Vector (x=0.0, y=1.0)

147. """

149. return "Vector (x={v.x}, y={v.y})".format(v=self)

151. def __add__(self, other):

152. """

153. Surcharge de l'opérateur binaire `{self} + {other}`: l'instruction

154. sera résolue comme `self.__add__(other)`.

156. On construit une nouvelle instance de `Vector` à partir des

157. coordonnées propres à l'objet `self` et à l'autre opérande

158. `other`.

160. >>> A = Point(1,0); B = Point(1,1)

161. >>> print Vector(A,B) + Vector(B,O) # = Vector(A,O)

162. Vector (x=-1.0, y=0.0)

163. """

165. return Vector(O, Point(self.x + other.x, self.y + other.y))

167. def __sub__(self, other):

168. """

169. Surcharge de l'opérateur binaire `{self} - {other}`: l'instruction

170. sera résolue comme `self.__sub__(other)`.

172. Attention: ne surcharge pas l'opérateur unaire `-{self}`, géré

173. par `__neg__`.

175. >>> A = Point(1,0); B = Point(1,1)

176. >>> print Vector(A,B) - Vector(A,B) # Différence

177. Vector (x=0.0, y=0.0)

178. >>> -Vector(A,B)                    # Négation

179. Traceback (most recent call last):

180. ...

181. TypeError: bad operand type for unary -: 'Vector'

182. """

184. return Vector(O, Point(self.x - other.x, self.y - other.y))

186. def __eq__(self, other):

187. """

188. Surcharge du test d'égalité `{self}=={other}`: l'instruction sera

189. résolue comme `self.__eq__(other)`.

191. >>> Vector(O,Point(0,1)) == Vector(Point(1,0),Point(1,1))

192. True

193. """

195. # On teste ici la nullité de la différence des 2

196. # vecteurs. D'autres tests auraient été possibles -- égalité

197. # des coordonnées, nullité de la norme de la différence,

198. # etc. -- mais on tire profit de la méthode héritée

199. # `Point.isOrigin()` testant la nullité des coordonnées (à la

200. # précision numérique près).

201. return (self - other).isOrigin()

203. def __abs__(self):

204. """

205. Surcharge la fonction `abs()` pour retourner la norme du vecteur.

207. >>> abs(Vector(Point(1,0), Point(1,1)))

208. 1.0

209. """

211. # On pourrait utiliser sqrt(self.sqnorm), mais c'est pour

212. # illustrer l'utilisation de la méthode héritée

213. # `Point.distance`...

214. return Point.distance(self, O)

216. def rotate(self, angle, deg=False):

217. """

218. Rotation (dans le sens trigonométrique) du vecteur par un `angle`,

219. exprimé en radians ou en degrés.

221. >>> Vector(Point(1,0),Point(1,1)).rotate(90,deg=True) == Vector(O,Point(-1,0))

222. True

223. """

225. from cmath import rect  # Bibliothèque de fonctions complexes

227. # On calcule la rotation en passant dans le plan complexe

228. z = complex(self.x, self.y)

229. phase = angle if not deg else angle / 57.29577951308232  # [rad]

230. u = rect(1., phase)  # exp(i*phase)

231. zu = z * u  # Rotation complexe

233. return Vector(O, Point(zu.real, zu.imag))

236. def circumscribedCircle(M, N, P):

237. """

238. Calcule le centre et le rayon du cercle circonscrit aux points

239. M,N,P.

241. Retourne: (centre [Point],rayon [float])

243. Lève une exception `ValueError` si le rayon ou le centre du cercle

244. circonscrit n'est pas défini.

246. >>> M = Point(-1,0); N = Point(1,0); P = Point(0,1)

247. >>> C,r = circumscribedCircle(M,N,P) # Centre O, rayon 1

248. >>> print C.distance(O), r

249. 0.0 1.0

250. >>> circumscribedCircle(M,O,N)       # Indéfini

251. Traceback (most recent call last):

252. ...

253. ValueError: Undefined circumscribed circle radius.

254. """

256. from math import sqrt

258. MN = Vector(M, N)

259. NP = Vector(N, P)

260. PM = Vector(P, M)

262. # Rayon du cercle circonscrit

263. m = abs(NP)  # |NP|

264. n = abs(PM)  # |PM|

265. p = abs(MN)  # |MN|

267. d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)

268. if d > 0:

269. rad = m * n * p / sqrt(d)

270. else:

271. raise ValueError("Undefined circumscribed circle radius.")

273. # Centre du cercle circonscrit

274. d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)

275. if d == 0:

276. raise ValueError("Undefined circumscribed circle center.")

278. om2 = Vector(O, M).sqnorm  # |OM|**2

279. on2 = Vector(O, N).sqnorm  # |ON|**2

280. op2 = Vector(O, P).sqnorm  # |OP|**2

282. x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d

283. y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d

285. return (Point(x0, y0), rad)  # (centre [Point], R [float])

288. if __name__ == '__main__':

290. # start-argparse

291. import argparse

293. parser = argparse.ArgumentParser(

294. usage="%(prog)s [-p/--plot] [-i/--input coordfile | x1,y1 x2,y2 x3,y3]",

295. description=__doc__)

296. parser.add_argument('coords', nargs='*', type=str, metavar='x,y',

297. help="Coordinates of point")

298. parser.add_argument('-i', '--input', nargs='?', type=file,

299. help="Coordinate file (one 'x,y' per line)")

300. parser.add_argument('-p', '--plot', action="store_true", default=False,

301. help="Draw the circumscribed circle")

302. parser.add_argument('--version', action='version', version=__version__)

304. args = parser.parse_args()

305. # end-argparse

307. if args.input:  # Lecture des coordonnées du fichier d'entrée

308. # Le fichier a déjà été ouvert en lecture par argparse (type=file)

309. args.coords = [coords for coords in args.input

310. if not coords.strip().startswith('#')]

312. if len(args.coords) != 3:  # Vérifie le nb de points

313. parser.error("Specify 3 points by their coordinates 'x,y' (got {})"

314. .format(len(args.coords)))

316. points = []  # Liste des points

317. for i, arg in enumerate(args.coords, start=1):

318. try:  # Déchiffrage de l'argument 'x,y'

319. x, y = (float(t) for t in arg.split(','))

320. except ValueError:

321. parser.error(

322. "Cannot decipher coordinates #{}: '{}'".format(i, arg))

324. points.append(Point(x, y))  # Création du point et ajout à la liste

325. print "#{:d}: {}".format(i, points[-1])  # Affichage du dernier point

327. # Calcul du cercle cisconscrit (lève une ValueError en cas de problème)

328. center, radius = circumscribedCircle(*points)  # Délistage

329. print "Circumscribed circle: {}, radius: {}".format(center, radius)

331. if args.plot:  # Figure

332. import matplotlib.pyplot as P

334. fig = P.figure()

335. ax = fig.add_subplot(1, 1, 1, aspect='equal')

336. # Points

337. ax.plot([p.x for p in points], [p.y for p in points], 'ko')

338. for i, p in enumerate(points, start=1):

339. ax.annotate("#{}".format(i), (p.x, p.y),

340. xytext=(5, 5), textcoords='offset points')

341. # Cercle circonscrit

342. c = P.matplotlib.patches.Circle((center.x, center.y), radius=radius,

343. fc='none', ec='k')

344. ax.add_patch(c)

345. ax.plot(center.x, center.y, 'r+')

347. P.show()

1. $ ./circonscrit.py -h

2. usage: circonscrit.py [-p/--plot] [-i/--input coordfile | x1,y1 x2,y2 x3,y3]

4. positional arguments:

5. x,y                   Coordinates of point

7. optional arguments:

8. -h, --help            show this help message and exit

9. -i [INPUT], --input [INPUT]

10. Coordinate file (one 'x,y' per line)

11. -p, --plot            Draw the circumscribed circle

12. --version             show program's version number and exit

14. $ ./circonscrit.py -p 0,0 1,0 1,1

Source : circonscrit.py

VIII-D - Matplotlib

VIII-D-1 - Figure (relativement) simple

Exemple de figure (charte graphique : seaborn)

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2014-12-15 16:06:44 ycopin>

5. """

6. Exemple un peu plus complexe de figure, incluant 2 axes, légendes, axes, etc.

7. """

9. import numpy as N

10. import matplotlib.pyplot as P

11. try:

12. import seaborn              # Amélioration de la charte graphique

13. except ImportError:

14. print u"Seaborn n'est pas accessible, charte graphique par défaut."

16. x = N.linspace(-N.pi, 3*N.pi, 2*360)

18. # Signal carré

19. y = N.sign(N.sin(x))            # = ± 1

21. # 3 premiers termes de la décomposition en série de Fourier

22. y1 = 4/N.pi * N.sin(x)          # Fondamentale

23. y2 = 4/N.pi * N.sin(3*x) / 3    # 1re harmonique

24. y3 = 4/N.pi * N.sin(5*x) / 5    # 2de harmonique

25. # Somme des 3 premières composantes

26. ytot = y1 + y2 + y3

28. # Figure

29. fig = P.figure()                # Création de la Figure

31. # 1er axe: composantes

32. ax1 = fig.add_subplot(2, 1, 1,  # 1er axe d'une série de 2 × 1

33. ylabel="Composantes",

34. title=u"Décomposition en série de Fourier")

35. ax1.plot(x, y1, label="Fondamental")

36. ax1.plot(x, y2, label=u"1re harmonique")

37. ax1.plot(x, y3, label=u"2de harmonique")

38. ax1.legend(loc="upper left", fontsize="x-small")

40. # 2nd axe: décomposition

41. ax2 = fig.add_subplot(2, 1, 2,  # 2d axe d'une série de 2 × 1

42. ylabel=u"Décomposition",

43. xlabel="x [rad]")

44. ax2.plot(x, y, lw=2, color='k', label=u"Signal carré")

45. ax2.plot(x, ytot, lw=2, ls=':', color='k', label=u"Somme des composantes")

46. ax2.legend(loc="upper left", fontsize="x-small")

48. # Sauvegarde de la figure (pas d'affichage intéractif)

49. fig.savefig("figure.png")

Source : figure.py

VIII-D-2 - Filtres du 2d ordre

Filtre passe-haut du 2d ordre.

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. import numpy as N

5. import matplotlib.pyplot as P

8. def passeBas(x, Q=1):

9. """

10. Filtre passe-bas en pulsation réduite *x* = omega/omega0, facteur de

11. qualité *Q*.

12. """

14. return 1 / (1 - x ** 2 + x / Q * 1j)

17. def passeHaut(x, Q=1):

19. return -x ** 2 / (1 - x ** 2 + x / Q * 1j)

22. def passeBande(x, Q=1):

24. return 1 / (1 + Q * (x - 1 / x) * 1j)

27. def coupeBande(x, Q=1):

29. return (1 - x ** 2) / (1 - x ** 2 + x / Q * 1j)

32. def gainNphase(f, dB=True):

33. """

34. Retourne le gain (éventuellement en dB) et la phase [rad] d'un

35. filtre de fonction de transfert complexe *f*.

36. """

38. g = N.abs(f)                        # Gain

39. if dB:                              # [dB]

40. g = 20 * N.log10(g)

41. p = N.angle(f)                      # [rad]

43. return g, p

46. def asympGain(x, pentes=(0, -40)):

48. lx = N.log10(x)

49. return N.where(lx < 0, pentes[0] * lx, pentes[1] * lx)

52. def asympPhase(x, phases=(0, -N.pi)):

54. return N.where(x < 1, phases[0], phases[1])

57. def diagBode(x, filtres, labels,

58. title='', plim=None, gAsymp=None, pAsymp=None):

59. """

60. Trace le diagramme de Bode -- gain [dB] et phase [rad] -- des filtres

61. de fonction de transfert complexe *filtres* en fonction de la pulsation

62. réduite *x*.

63. """

65. fig = P.figure()

66. axg = fig.add_subplot(2, 1, 1,        # Axe des gains

67. xscale='log',

68. ylabel='Gain [dB]')

69. axp = fig.add_subplot(2, 1, 2,        # Axe des phases

70. sharex=axg,

71. xlabel=r'x = $\omega$/$\omega_0$', xscale='log',

72. ylabel='Phase [rad]')

74. lstyles = ['--', '-', '-.', ':']

75. for f, label, ls in zip(filtres, labels, lstyles):  # Tracé des courbes

76. g, p = gainNphase(f, dB=True)       # Calcul du gain et de la phase

77. axg.plot(x, g, lw=2, ls=ls, label="Q=" + str(label))  # Gain

78. axp.plot(x, p, lw=2, ls=ls)                           # Phase

80. # Asymptotes

81. if gAsymp is not None:              # Gain

82. axg.plot(x, asympGain(x, gAsymp), 'k:', lw=2, label='_')

83. if pAsymp is not None:              # Phase

84. #axp.plot(x, asympPhase(x,pAsymp), 'k:')

85. pass

87. axg.legend(loc='best', prop=dict(size='small'))

89. # Labels des phases

90. axp.set_yticks(N.arange(-2, 2.1) * N.pi / 2)

91. axp.set_yticks(N.arange(-4, 4.1) * N.pi / 4, minor=True)

92. axp.set_yticklabels([r'$-\pi$', r'$-\pi/2$', r'$0$', r'$\pi/2$', r'$\pi$'])

93. # Domaine des phases

94. if plim is not None:

95. axp.set_ylim(plim)

97. # Ajouter les grilles

98. for ax in (axg, axp):

99. ax.grid()                           # x et y, majors

100. ax.grid(which='minor')              # x et y, minors

102. # Ajustements fins

103. gmin, gmax = axg.get_ylim()

104. axg.set_ylim(gmin, max(gmax, 3))

106. fig.subplots_adjust(hspace=0.1)

107. axg.xaxis.set_major_formatter(P.matplotlib.ticker.ScalarFormatter())

108. P.setp(axg.get_xticklabels(), visible=False)

110. if title:

111. axg.set_title(title)

113. return fig

115. if __name__ == '__main__':

117. #P.rc('mathtext', fontset='stixsans')

119. x = N.logspace(-1, 1, 1000)              # de 0.1 à 10 en 1000 pas

121. # Facteurs de qualité

122. qs = [0.25, 1 / N.sqrt(2), 5]            # Valeurs numériques

123. labels = [0.25, r'$1/\sqrt{2}$', 5]      # Labels

125. # Calcul des fonctions de transfert complexes

126. pbs = [ passeBas(x, Q=q) for q in qs ]

127. phs = [ passeHaut(x, Q=q) for q in qs ]

128. pcs = [ passeBande(x, Q=q) for q in qs ]

129. cbs = [ coupeBande(x, Q=q) for q in qs ]

131. # Création des 4 diagrammes de Bode

132. figPB = diagBode(x, pbs, labels, title='Filtre passe-bas',

133. plim=(-N.pi, 0),

134. gAsymp=(0, -40), pAsymp=(0, -N.pi))

135. figPH = diagBode(x, phs, labels, title='Filtre passe-haut',

136. plim=(0, N.pi),

137. gAsymp=(40, 0), pAsymp=(N.pi, 0))

138. figPC = diagBode(x, pcs, labels, title='Filtre passe-bande',

139. plim=(-N.pi / 2, N.pi / 2),

140. gAsymp=(20, -20), pAsymp=(N.pi / 2, -N.pi / 2))

141. figCB = diagBode(x, cbs, labels, title='Filtre coupe-bande',

142. plim=(-N.pi / 2, N.pi / 2),

143. gAsymp=(0, 0), pAsymp=(0, 0))

145. P.show()

Source : filtres2ndOrdre.py

IX - Exercices

Les exercices sont de difficulté variable, de ★ (simple) à ★★★ (complexe).

IX-A - Introduction

IX-A-1 - Intégration : méthode des rectangles ★

La méthode des rectangles permet d'approximer numériquement l'intégrale d'une fonction f :

On définit la fonction sq renvoyant le carré d'un nombre par (cf. Fonctions) :

def sq(x) :
return x*x

Écrire un programme calculant l'intégrale de cette fonction entre a=0 et b=1, en utilisant une subdivision en n=100 pas dans un premier temps. Quelle est la précision de la méthode, et comment dépend-elle du nombre de pas ?

IX-A-2 - Fizz Buzz ★

Écrire un programme jouant au Fizz Buzz jusqu'à 99 :

1 2 Fizz! 4 Buzz! Fizz! 7 8 Fizz! Buzz! 11 Fizz! 13 14 Fizz Buzz! 16...

IX-A-3 - PGCD : algorithme d'Euclide ★★

Écrire un programme calculant le PGCD de deux nombres (p.ex. 306 et 756) par l'algorithme d'Euclide.

IX-A-4 - Tables de multiplication ★

Écrire un programme affichant les tables de multiplication :

1. 1 x 1 = 1

2. 1 x 2 = 2

3. ...

4. 9 x 9 = 81

IX-B - Manipulation de listes

IX-B-1 - Crible d'Ératosthène ★

Implémenter le crible d'Ératosthène pour afficher les nombres premiers compris entre 1 et un entier fixe, p.ex. :

Liste des entiers premiers <= 41
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]

IX-B-2 - Carré magique ★★

Un carré magique d'ordre n est un tableau carré n × n dans lequel on écrit une et une seule fois les nombres entiers de 1 à n², de sorte que la somme des n nombres de chaque ligne, colonne ou diagonale principale soit constante. P.ex. le carré magique d'ordre 5, où toutes les sommes sont égales à 65 :

11

18

25

2

9

10

12

19

21

3

4

6

13

20

22

23

5

7

14

16

17

24

1

8

15

Pour les carrés magiques d'ordre impair, on dispose de l'algorithme suivant - (i,j) désignant la case de la ligne i, colonne j du carré ; on se place en outre dans une indexation « naturelle » commençant à 1 :

1. La case (n,(n+1)/2) contient 1 ;

2. Si la case (i,j) contient la valeur k, alors on place la valeur k+1 dans la case (i+1,j+1) si cette case est vide, ou dans la case (i-1,j) sinon. On respecte la règle selon laquelle un indice supérieur à n est ramené à 1.

Programmer cet algorithme pour pouvoir construire un carré magique d'ordre impair quelconque.

IX-C - Programmation

IX-C-1 - Suite de Syracuse (fonction) ★

Écrire une fonction suite_syracuse(n) retournant la (partie non triviale de la) suite de Syracuse pour un entier n. Écrire une fonction temps_syracuse(n, altitude=False) retournant le temps de vol (éventuellement en altitude) correspondant à l'entier n. Tester ces fonctions sur n=15 :

1. >>> suite_syracuse(15)

2. [15, 46, 23, 70, 35, 106, 53, 160, 80, 40, 20, 10, 5, 16, 8, 4, 2, 1]

3. >>> temps_syracuse(15)

4. 17

5. >>> temps_syracuse(15, altitude=True)

6. 10

IX-C-2 - Flocon de Koch (programmation récursive) ★★★

En utilisant les commandes left, right et forward de la bibliothèque graphique standard turtle dans une fonction récursive, générer à l'écran un flocon de Koch d'ordre arbitraire.

Flocon de Koch d'ordre 3.

IX-C-3 - Jeu du plus ou moins (exceptions) ★

Écrire un jeu de « plus ou moins »:

1. Vous devez deviner un nombre entre 1 et 100.

2. Votre proposition: 27

3. C'est plus.

4. [...]

5. Vous avez trouvé en 6 coups !

La solution sera générée aléatoirement par la fonction random.randint(). Le programme devra être robuste aux entrées invalides (« toto », 120, etc.), et aux lâches abandons par interruption (KeyboardInterrupt).

IX-C-4 - Animaux (POO/TDD) ★

Téléchargez animaux.py et complétez les classes Animal et Chien pour qu'elles passent avec succès tous les tests (voir Développement piloté par les tests). On appellera les tests via la ligne de commande :

py.test animaux.py

IX-C-5 - Jeu de la vie (POO) ★★

On se propose de programmer l'automate cellulaire le plus célèbre, le Jeu de la vie.

Pour cela, vous créerez une classe Life qui contiendra la grille du jeu ainsi que les méthodes qui permettront son évolution. Vous initialiserez la grille aléatoirement à l'aide de la fonction random.choice(), et vous afficherez l'évolution de l'automate dans la sortie standard du terminal, p.ex. :

1. ...#..#.....##.......

2. .....###.............

3. #........#...........

4. .....#...#...........

5. ................##...

6. .....#.#......##..#..

7. ..............##.##..

8. ..............##.##..

9. ................#....

Pour que l'affichage soit agréable à l'œil, vous marquerez des pauses entre l'affichage de chaque itération grâce à la fonction time.sleep().

IX-D - Manipulation de tableaux (arrays)

IX-D-1 - Inversion de matrice ★

Créer un tableau carré réel rr\mathsf{r} aléatoire (numpy.random.randn()), calculer la matrice hermitienne (numpy.dot()), l'inverser (numpy.linalg.inv()), et vérifier que (numpy.eye()) à la précision numérique près (numpy.allclose()).

IX-D-2 - Median Absolute Deviation ★

En statistique, le Median Absolute Deviation (MAD) est un estimateur robuste de la dispersion d'un échantillon 1D : MAD = median(| x - median(x) |).

À l'aide des fonctions numpy.median() et numpy.abs(), écrire une fonction mad(x, axis=None) calculant le MAD d'un tableau, éventuellement le long d'un ou plusieurs de ses axes.

IX-D-3 - Distribution du pull ★★★

Le pull est une quantité statistique permettant d'évaluer la conformité des erreurs par rapport à une distribution de valeurs (typiquement les résidus d'un ajustement). Pour un échantillon et les erreurs associées , le pull est défini par :

- moyenne optimale (pondérée par la variance) :  ;

- erreur sur la moyenne pondérée :  ;

- définition du pull : , où et sont calculées sans le point i.

Si les erreurs sont correctes, la distribution du pull est centrée sur 0 avec une déviation standard de 1.

Écrire une fonction pull(x, dx) calculant le pull de tableaux 1D.

IX-E - Méthodes numériques

IX-E-1 - Quadrature et zéro d'une fonction ★

À l'aide des algorithmes disponibles dans scipy :

- calculer numériquement l'intégrale  ;

- résoudre numériquement l'équation .

IX-E-2 - Schéma de Romberg ★★

Écrire une fonction integ_romberg(f, a, b, epsilon=1e-6) permettant de calculer l'intégrale numérique de la fonction f entre les bornes a et b avec une précision epsilon selon la méthode de Romberg.

Tester sur des solutions analytiques et en comparant à scipy.integrate.romberg().

IX-E-3 - Méthode de Runge-Kutta ★★

Développer un algorithme permettant d'intégrer numériquement une équation différentielle du 1er ordre en utilisant la méthode de Runge-Kutta d'ordre quatre.

Tester sur des solutions analytiques et en comparant à scipy.integrate.odeint().

IX-F - Visualisation (matplotlib)

IX-F-1 - Quartet d'Anscombe ★

Après chargement des données, calculer et afficher les propriétés statistiques des quatre jeux de données du Quartet d'Anscombe :

- moyenne et variance des x et des y (numpy.mean() et numpy.var()) ;

- corrélation entre les x et les y (scipy.stats.pearsonr()) ;

- équation de la droite de régression linéaire y = ax + b (scipy.stats.linregress())

Quartet d'Anscombe

I

II

III

IV

x

y

x

y

x

y

x

y

10.0

8.04

10.0

9.14

10.0

7.46

8.0

6.58

8.0

6.95

8.0

8.14

8.0

6.77

8.0

5.76

13.0

7.58

13.0

8.74

13.0

12.74

8.0

7.71

9.0

8.81

9.0

8.77

9.0

7.11

8.0

8.84

11.0

8.33

11.0

9.26

11.0

7.81

8.0

8.47

14.0

9.96

14.0

8.10

14.0

8.84

8.0

7.04

6.0

7.24

6.0

6.13

6.0

6.08

8.0

5.25

4.0

4.26

4.0

3.10

4.0

5.39

19.0

12.50

12.0

10.84

12.0

9.13

12.0

8.15

8.0

5.56

7.0

4.82

7.0

7.26

7.0

6.42

8.0

7.91

5.0

5.68

5.0

4.74

5.0

5.73

8.0

6.89

Pour chacun des jeux de données, tracer y en fonction de x, ainsi que la droite de régression linéaire.

IX-F-2 - Diagramme de bifurcation : la suite logistique ★★

Écrivez une fonction qui calcule la valeur d'équilibre de la suite logistique pour un (nécessairement compris entre 0 et 1) et un paramètre (parfois noté ) donné.

Générez l'ensemble de ces points d'équilibre pour des valeurs de comprises entre 0 et 4 :

Diagramme de bifurcation.

N.B. Vous utiliserez la bibliothèque Matplotlib pour tracer vos résultats.

IX-F-3 - Ensemble de Julia ★★

Représentez l'ensemble de Julia pour la constante complexe  :

Ensemble de Julia pour c = 0.284 + 0.0122j.

On utilisera la fonction numpy.meshgrid() pour construire le plan complexe, et l'on affichera le résultat grâce à la fonction matplotlib.pyplot.imshow().

Voir également : Superposition d'ensembles de Julia

IX-G - Mise en œuvre de l'ensemble des connaissances acquises

IX-G-1 - Équation différentielle ★

À l'aide de la fonction scipy.integrate.odeint(), intégrer les équations du mouvement d'un boulet de canon soumis à des forces de frottement « turbulentes » (en ) :

Utiliser les valeurs numériques pour un boulet de canon de 36 livres :

1. g = 9.81       # Pesanteur [m/s2]

2. cx = 0.45      # Coefficient de frottement d'une sphère

3. rhoAir = 1.2   # Masse volumique de l'air [kg/m3]

4. rad = 0.1748/2 # Rayon du boulet [m]

5. rho = 6.23e3   # Masse volumique du boulet [kg/m3]

6. mass = 4./3.*N.pi*rad**3 * rho           # Masse du boulet [kg]

7. alpha = 0.5*cx*rhoAir*N.pi*rad**2 / mass # Coeff. de frottement / masse

8. v0 = 450.      # Vitesse initiale [m/s]

9. alt = 45.      # Inclinaison du canon [deg]

IX-H - Exercices en vrac

- Exercices de base

- Entraînez-vous !

- Learn Python The Hard Way

- Google Code Jam

- CheckIO

X - Solutions aux exercices

X-A - Méthode des rectangles

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2014-09-17 19:30:18 ycopin>

5. # Division réelle de type Python 3 - ligne admise

6. from __future__ import division

8. """

9. Définition de la fonction carré et calcul de son intégrale entre 0

10. et 1 dans le main par la méthode des rectangles (subdivision en 100

11. pas)

12. """

14. __author__ = "Adrien Licari <adrien.licari@ens-lyon.fr>"

16. # Définition de la fonction sq, admise au stade du TD 1

19. def sq(x):

20. return x * x

22. # Début du programme principal (main)

24. # On définit les bornes d'intégration a et b, le nombre de pas n

25. a = 0

26. b = 1

27. n = 100

29. # Largeur des rectangles dx

30. h = (b - a) / n      # Division réelle!

32. integ = 0         # Cette variable accumulera les aires des rectangles

34. # On sait déjà que l'on va calculer n aires de rectangles, donc une

35. # boucle for est appropriée

36. for i in range(n):            # Boucle de 0 à n-1

37. x = a + (i + 0.5) * h         # Abscisse du rectangle

38. integ += sq(x) * h        # On ajoute au total l'aire du rectangle

40. print "Intégrale de x**2 entre a =", a, "et b =", b, "avec n =", n

41. # On affiche le résultat numérique

42. print "Résultat numérique: ", integ

43. theorie = (b ** 3 - a ** 3) / 3                     # Résultat analytique

44. # On affiche le résultat analytique

45. print "Résultat analytique:", theorie

46. print "Erreur relative:", (integ / theorie - 1)  # On affiche l'erreur

$ python integ.py

Intégrale de x**2 entre a = 0 et b = 1 avec n = 100
Résultat numérique:  0.333325
Résultat analytique: 0.333333333333
Erreur relative: -2.49999999998e-05

Source : integ.py

X-B - Fizz Buzz

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2017-05-28 21:42 ycopin@lyonovae03.in2p3.fr>

5. """

6. Jeu du Fizz Buzz

7. """

9. __author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"

12. for i in range(1, 100):                    # Entiers de 1 à 99

13. if ((i % 3) == 0) and ((i % 5) == 0):  # Multiple de 3 et 5

14. print 'Fizz Buzz!',                # Affichage sans retour à la ligne

15. elif ((i % 3) == 0):                   # Multiple de 3 uniquement

16. print 'Fizz!',

17. elif ((i % 5) == 0):                   # Multiple de 5 uniquement

18. print 'Buzz!',

19. else:

20. print i,

$ python fizz.py

1 2 Fizz! 4 Buzz! Fizz! 7 8 Fizz! Buzz! 11 Fizz! 13 14 Fizz Buzz! 16...

Source : fizz.py

X-C - Algorithme d'Euclide

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2014-09-17 18:25:00 ycopin>

5. """

6. Calcul du PGCD de deux entiers

7. """

9. __author__ = "Adrien Licari <adrien.licari@ens-lyon.fr>"

11. # Entiers dont on calcule le PGCD

12. a = 306

13. b = 756

15. # Test usuels : les nombres sont bien positifs et a > b

16. assert a > 0

17. assert b > 0

18. if (a < b):

19. a, b = b, a                   # Interversion

21. a0, b0 = a, b                     # On garde une copie des valeurs originales

23. # On boucle jusqu'à ce que le reste soit nul, d'où la boucle while. Il

24. # faut être sûr que l'algorithme converge dans tous les cas !

25. while (b != 0):

26. # On remplace a par b et b par le reste de la division euclidienne

27. # de a par b

28. a, b = b, a % b

30. print 'Le PGCD de', a0, 'et', b0, 'vaut', a  # On affiche le résultat

31. # Vérifications

32. print a0 // a, 'x', a, '=', (a0 // a * a)  # a//b: division euclidienne

33. print b0 // a, 'x', a, '=', (b0 // a * a)

$ python pgcd.py

Le PGCD de 756 et 306 vaut 18
42 x 18 = 756
17 x 18 = 306

Source : pgcd.py

X-D - Crible d'Ératosthène

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. """

5. Crible d'Ératosthène.

7. Source: http://fr.wikibooks.org/wiki/Exemples_de_scripts_Python#Implémentation_du_crible_d'Ératosthène

8. """

10. __author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"

12. # start-sys

13. # Gestion simplifiée d'un argument entier sur la ligne de commande

14. import sys

16. if sys.argv[1:]:  # Présence d'au moins un argument sur la ligne de commande

17. try:

18. n = int(sys.argv[1])  # Essayer de lire le 1er argument comme un entier

19. except ValueError:

20. raise ValueError("L'argument '{}' n'est pas un entier"

21. .format(sys.argv[1]))

22. else:                        # Pas d'argument sur la ligne de commande

23. n = 101                  # Valeur par défaut

24. # end-sys

26. # Liste des entiers *potentiellement* premiers. Les nb non premiers

27. # seront étiquetés par 0 au fur et à mesure.

28. l = range(n + 1)                          # <0,...,n>, 0 n'est pas premier

29. l[1] = 0                                # 1 n'est pas premier

31. i = 2                                   # Entier à tester

32. while i ** 2 <= n:                        # Inutile de tester jusqu'à n

33. if l[i] != 0:                       # Si i n'est pas étiqueté (=0)...

34. # ...étiqueter tous les multiples de i

35. l[2 * i::i] = [0] * len(l[2 * i::i])

36. i += 1                              # Passer à l'entier à tester suivant

38. # Afficher la liste des entiers premiers (c-à-d non étiquetés)

39. print "Liste des entiers premiers <= {} :".format(n)

40. print [i for i in l if i != 0]

$ python crible.py

Liste des entiers premiers <= 101
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]

Source : crible.py

X-E - Carré magique

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2017-07-13 14:29 ycopin@lyonovae03.in2p3.fr>

5. """

6. Création et affichage d'un carré magique d'ordre impair.

7. """

9. __author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"

11. n = 5                                   # Ordre du carré magique

13. # On vérifie que l'ordre est bien impair

14. assert n % 2 == 1, "L'ordre {} n'est pas impair".format(n)

16. # Le carré sera stocké dans une liste de n listes de n entiers

17. # Initialisation du carré: liste de n listes de n zéros.

18. array = [[0 for j in range(n)] for i in range(n)]

20. # Initialisation de l'algorithme

21. i, j = n, (n + 1) / 2         # Indices de l'algo (1-indexation)

22. array[i - 1][j - 1] = 1       # Attention: python utilise une 0-indexation

24. # Boucle sur valeurs restant à inclure dans le carré (de 2 à n**2)

25. for k in range(2, n**2 + 1):

26. # Test de la case i+1, j+1 (modulo n)

27. i2 = (i + 1) % n

28. j2 = (j + 1) % n

29. if array[i2 - 1][j2 - 1] == 0:  # La case est vide: l'utiliser

30. i, j = i2, j2

31. # La case est déjà remplie : utiliser la case i-1, j

32. else:

33. i = (i - 1) % n

34. array[i - 1][j - 1] = k       # Remplissage de la case

36. # Affichage, avec vérification des sommes par ligne et par colonne

37. print "Carré magique d'ordre {}:".format(n)

38. for row in array:

39. print '  '.join('{:2d}'.format(k) for k in row), '=', sum(row)

40. print '  '.join('==' for k in row)

41. print '  '.join(str(sum(array[i][j] for i in range(n))) for j in range(n))

42. $ python carre.py

$ python carre.py

Carré magique d'ordre 5:
11  18  25   2   9 = 65
10  12  19  21   3 = 65
4   6  13  20  22 = 65
23   5   7  14  16 = 65
17  24   1   8  15 = 65
==  ==  ==  ==  ==
65  65  65  65  65

Source : carre.py

X-F - Suite de Syracuse

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2014-09-24 17:31:53 ycopin>

5. __author__ = "Adrien Licari <adrien.licari@ens-lyon.fr>; Yannick Copin <y.copin@ipnl.in2p3.fr>"

8. def suite_syracuse(n):

9. """

10. Retourne la suite de Syracuse pour l'entier n.

12. >>> suite_syracuse(15)

13. [15, 46, 23, 70, 35, 106, 53, 160, 80, 40, 20, 10, 5, 16, 8, 4, 2, 1]

14. """

16. seq = [n]                      # La suite de Syracuse sera complétée...

17. while (seq[-1] != 1):         # ...jusqu'à tomber sur 1

18. if seq[-1] % 2 == 0:      # u_n est pair

19. seq.append(seq[-1] // 2)  # Division euclidienne par 2

20. else:                     # u_n est impair

21. seq.append(seq[-1] * 3 + 1)

23. return seq

26. def temps_syracuse(n, altitude=False):

27. """

28. Calcule le temps de vol (éventuellement en altitude) de la suite

29. de Syracuse pour l'entier n.

31. >>> temps_syracuse(15)

32. 17

33. >>> temps_syracuse(15, altitude=True)

34. 10

35. """

37. seq = suite_syracuse(n)

38. if not altitude:            # Temps de vol total

39. return len(seq) - 1

40. else:                       # Temps de vol en altitude

41. # Construction de la séquence en altitude

42. alt = []

43. for i in seq:

44. if i >= n:

45. alt.append(i)

46. else:

47. break

48. return len(alt) - 1

50. if __name__ == '__main__':

52. n = 15

53. print "Suite de Syracuse pour n =", n

54. print suite_syracuse(n)

55. print "Temps de vol total:      ", temps_syracuse(n)

56. print "Temps de vol en altitude:", temps_syracuse(n, altitude=True)

Suite de Syracuse pour n = 15
[15, 46, 23, 70, 35, 106, 53, 160, 80, 40, 20, 10, 5, 16, 8, 4, 2, 1]
Temps de vol total:       17
Temps de vol en altitude: 10

Source : syracuse.py

X-G - Flocon de Koch

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. from __future__ import division  # Pas de division euclidienne par défaut

6. """

7. Tracé (via 'turtle') d'un flocon de Koch d'ordre arbitraire.

9. Dans le même genre:

11. - courbe de Peano (http://fr.wikipedia.org/wiki/Courbe_de_Peano)

12. - courbe de Hilbert (http://fr.wikipedia.org/wiki/Courbe_de_Hilbert)

13. - île de Gosper (http://fr.wikipedia.org/wiki/Île_de_Gosper)

15. Voir également:

17. - L-système: http://fr.wikipedia.org/wiki/L-système

18. - Autres exemples: http://natesoares.com/tutorials/python-fractals/

19. """

21. __version__ = "Time-stamp: <2013-01-14 00:49 ycopin@lyopc469>"

22. __author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"

24. import turtle as T

27. def koch(niveau=3, iter=0, taille=100, delta=0):

28. """

29. Tracé du flocon de Koch de niveau 'niveau', de taille 'taille'

30. (px).

32. Cette fonction récursive permet d'initialiser le flocon (iter=0,

33. par défaut), de tracer les branches fractales (0<iter<=niveau) ou

34. bien juste de tracer un segment (iter>niveau).

35. """

37. if iter == 0:                         # Dessine le triangle de niveau 0

38. T.title("Flocon de Koch - niveau {}".format(niveau))

39. koch(iter=1, niveau=niveau, taille=taille, delta=delta)

40. T.right(120)

41. koch(iter=1, niveau=niveau, taille=taille, delta=delta)

42. T.right(120)

43. koch(iter=1, niveau=niveau, taille=taille, delta=delta)

44. elif iter <= niveau:                  # Trace une section _/\_ du flocon

45. koch(iter=iter + 1, niveau=niveau, taille=taille, delta=delta)

46. T.left(60 + delta)

47. koch(iter=iter + 1, niveau=niveau, taille=taille, delta=delta)

48. T.right(120 + 2 * delta)

49. koch(iter=iter + 1, niveau=niveau, taille=taille, delta=delta)

50. T.left(60 + delta)

51. koch(iter=iter + 1, niveau=niveau, taille=taille, delta=delta)

52. else:                               # Trace le segment de dernier niveau

53. T.forward(taille / 3 ** (niveau + 1))

55. if __name__ == '__main__':

57. # start-argparse

58. # Exemple d'utilisation de la bibliothèque de gestion d'arguments 'argparse'

59. import argparse

61. desc = u"Tracé (via 'turtle') d'un flocon de Koch d'ordre arbitraire."

63. # Définition des options

64. parser = argparse.ArgumentParser(description=desc)

65. parser.add_argument('ordre', nargs='?', type=int,

66. help="Ordre du flocon, >0 [%(default)s]",

67. default=3)

68. parser.add_argument('-t', '--taille', type=int,

69. help="Taille de la figure, >0 [%(default)s px]",

70. default=500)

71. parser.add_argument('-d', '--delta', type=float,

72. help="Delta [%(default)s deg]",

73. default=0.)

74. parser.add_argument('-f', '--figure', type=str,

75. help="Nom de la figure de sortie (format EPS)")

76. parser.add_argument('-T', '--turbo',

77. action="store_true", default=False,

78. help="Mode Turbo")

80. # Déchiffrage des options et arguments

81. args = parser.parse_args()

83. # Quelques tests sur les args et options

84. if not args.ordre > 0:

85. parser.error("Ordre requis '{}' invalide".format(args.ordre))

87. if not args.taille > 0:

88. parser.error("La taille de la figure doit être positive")

89. # end-argparse

91. if args.turbo:

92. T.hideturtle()

93. T.speed(0)

95. # Tracé du flocon de Koch de niveau 3

96. koch(niveau=args.ordre, taille=args.taille, delta=args.delta)

97. if args.figure:

98. # Sauvegarde de l'image

99. print "Sauvegarde de la figure dans '{}'".format(args.figure)

100. T.getscreen().getcanvas().postscript(file=args.figure)

102. T.exitonclick()

Source : koch.py

X-H - Jeu du plus ou moins

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. import random

6. nmin, nmax = 1, 100

7. nsol = random.randint(nmin, nmax)

9. print "Vous devez deviner un nombre entre {} et {}.".format(nmin, nmax)

10. ncoups = 0

11. try:

12. while True:

13. proposition = raw_input("Votre proposition: ")

14. try:

15. ntest = int(proposition)

16. if not nmin <= ntest <= nmax:

17. raise ValueError("Proposition invalide")

18. except ValueError:

19. print "Votre proposition {!r} n'est pas un entier " \

20. "compris entre {} et {}.".format(proposition, nmin, nmax)

21. continue

22. ncoups += 1

23. if nsol > ntest:

24. print "C'est plus."

25. elif nsol < ntest:

26. print "C'est moins."

27. else:

28. print "Vous avez trouvé en {} coup{}!".format(

29. ncoups, 's' if ncoups > 1 else '')

30. break  # Interrompt la boucle while

31. except (KeyboardInterrupt, EOFError):  # Interception de Ctrl-C et Ctrl-D

32. print "\nVous abandonnez après seulement {} coup{}!".format(

33. ncoups, 's' if ncoups > 1 else '')

Source : pm.py

X-I - Animaux

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. import pytest                   # Module (non-- standard) de tests

7. class Animal(object):  # *object* est la classe dont dérivent toutes les autres

9. """

10. Classe définissant un `Animal`, caractérisé par son nom et son

11. poids.

12. """

14. def __init__(self, nom, masse):

15. """

16. Méthode d'instanciation à partir d'un nom (str) et d'un poids (float).

17. """

19. # Ici, convertir les paramètres pour être sûr qu'il ont le bon

20. # type. On utilisera `str` et `float`

21. self.nom = str(nom)

22. self.masse = float(masse)

24. self.vivant = True       # Les animaux sont vivants à l'instanciation

25. self.empoisonne = False  # Animal empoisonné ?

27. def __str__(self):

28. """

29. Surcharge de l'opérateur `str`: l'affichage *informel* de l'objet

30. dans l'interpréteur, p.ex. `print self` sera résolu comme

31. `self.__str__()`

33. Retourne une chaîne de caractères.

34. """

36. return "{a.nom} ({a.masse:.1f} kg)".format(a=self)

38. def estVivant(self):

39. """Méthode booléenne, vraie si l'animal est vivant."""

41. return self.vivant

43. def mourir(self):

44. """Change l'état interne de l'objet (ne retourne rien)."""

46. self.vivant = False

48. def __cmp__(self, other):

49. """

50. Surcharge des opérateurs de comparaison, sur la base de la masse

51. des animaux.

52. """

54. return cmp(self.masse, other.masse)

56. def __call__(self, other):

57. """

58. Surcharge de l'opérateur '()' pour manger un autre animal (qui

59. meurt s'il est vivant) et prendre du poids (mais pas plus que

60. la masse de l'autre ou 10 % de son propre poids).  Attention aux

61. animaux empoisonnés !

63. L'instruction `self(other)` sera résolue comme

64. `self.__call__(other).

65. """

67. other.mourir()

68. poids = min(other.masse, self.masse * 0.1)

69. self.masse += poids

70. other.masse -= poids

71. if other.empoisonne:

72. self.mourir()

75. class Chien(Animal):

77. """

78. Un `Chien` hérite de `Animal` avec des méthodes additionnelles

79. (p.ex. l'aboiement et l'odorat).

80. """

82. def __init__(self, nom, masse=20, odorat=0.5):

83. """Définit un chien plus ou moins fin limier."""

85. # Initialisation de la classe parente

86. Animal.__init__(self, nom, masse)

88. # Attribut propre à la classe dérivée

89. self.odorat = float(odorat)

91. def __str__(self):

93. return "{a.nom} (Chien, {a.masse:.1f} kg)".format(a=self)

95. def aboyer(self):

96. """Une méthode bien spécifique aux chiens."""

98. print("Ouaf ! Ouaf !")

100. def estVivant(self):

101. """Quand on vérifie qu'un chien est vivant, il aboie."""

103. vivant = Animal.estVivant(self)

105. if vivant:

106. self.aboyer()

108. return vivant

110. # start-tests

113. def test_empty_init():

114. with pytest.raises(TypeError):

115. Animal()

118. def test_wrong_init():

119. with pytest.raises(ValueError):

120. Animal('Youki', 'lalala')

123. def test_init():

124. youki = Animal('Youki', 600)

125. assert youki.masse == 600

126. assert youki.vivant

127. assert youki.estVivant()

128. assert not youki.empoisonne

129. # end-tests

132. def test_str():

133. youki = Animal('Youki', 600)

134. assert str(youki) == 'Youki (600.0 kg)'

137. def test_mort():

138. youki = Animal('Youki', 600)

139. assert youki.estVivant()

140. youki.mourir()

141. assert not youki.estVivant()

144. def test_lt():

145. medor = Animal('Medor', 600)

146. kiki = Animal('Kiki', 20)

147. assert kiki < medor

148. with pytest.raises(AttributeError):

149. medor < 1

152. def test_mange():

153. medor = Animal('Medor', 600)

154. kiki = Animal('Kiki', 20)

155. medor(kiki)                 # Médor mange Kiki

156. assert medor.estVivant()

157. assert not kiki.estVivant()

158. assert kiki.masse == 0

159. assert medor.masse == 620

160. kiki = Animal("Kiki Jr.", 20)

161. kiki(medor)                 # Kiki Jr. mange Médor

162. assert not medor.estVivant()

163. assert kiki.estVivant()

164. assert kiki.masse == 22

165. assert medor.masse == 618   # Médor a perdu du poids en se faisant manger!

168. def test_init_chien():

169. medor = Chien('Medor', 600)

170. assert isinstance(medor, Animal)

171. assert isinstance(medor, Chien)

172. assert medor.odorat == 0.5

173. assert str(medor) == 'Medor (Chien, 600.0 kg)'

174. assert medor.estVivant()

=========================== test session starts ============================
platform linux2 -- Python 2.7.6 -- py-1.4.24 -- pytest-2.6.2
collected 8 items

animauxSol.py ........

========================= 8 passed in 0.04 seconds =========================

Source : animauxSol.py

X-J - Particules

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2014-09-04 16:56:47 alicari>

6. from __future__ import division  # division réelle de type python 3, admis

7. import pytest                    # pytest importé pour les tests unitaires

8. import math

10. """

11. Définition d'une classe point matériel, avec sa masse,

12. sa position et sa vitesse, et des méthodes pour le déplacer.

13. Le main test applique cela à un problème à force centrale

14. gravitationnel ou électrostatique.

15. Remarque : toutes les unités ont été choisies adimensionnées.

16. """

18. __author__ = "Adrien Licari <adrien.licari@ens-lyon.fr>"

21. # Un critère pour déterminer l'égalité entre réels

22. tolerance = 1e-8

25. #############################################################################

26. ### Définition de la classe Vector, utile pour la position et la vitesse. ###

27. #############################################################################

29. class Vector(object):

31. """

32. Une classe structure simple contenant 3 coordonnées.

33. Une méthode est disponible pour en calculer la norme et

34. une surcharge des opérateurs ==, !=, +, - et * est proposée.

35. """

37. def __init__(self, x=0, y=0, z=0):

38. """

39. Constructeur de la classe vector.

40. Par défaut, construit le vecteur nul.

42. Args:

43. x,y,z(float): Les composantes du vecteur à construire.

45. Raises:

46. TypeError en cas de composantes non réelles

47. """

48. try:

49. self.x = float(x)

50. self.y = float(y)

51. self.z = float(z)

52. except (ValueError, TypeError, AttributeError):

53. raise TypeError("The given coordinates must be numbers")

55. def __str__(self):

56. """

57. Surcharge de l'opérateur `str`.

59. Returns :

60. "(x,y,z)" with 2 decimals

61. """

62. return "({:.2f},{:.2f},{:.2f})".format(self.x, self.y, self.z)

64. def __eq__(self, other):

65. """

66. Surcharge de l'opérateur `==` pour tester l'égalité

67. entre deux vecteurs.

69. Args :

70. other(Vector): Un autre vecteur

72. Raises :

73. TypeError si other n'est pas un objet Vector

74. """

75. try:

76. return abs(self.x - other.x) < tolerance and \

77. abs(self.y - other.y) < tolerance and \

78. abs(self.z - other.z) < tolerance

79. except (ValueError, TypeError, AttributeError):

80. raise TypeError("Tried to compare Vector and non-Vector objects")

82. def __ne__(self, other):

83. """

84. Surcharge de l'opérateur `!=` pour tester l'inégalité

85. entre deux vecteurs.

87. Args :

88. other(Vector): Un autre vecteur

90. Raises :

91. TypeError si other n'est pas un objet Vector

92. """

93. return not self == other

95. def __add__(self, other):

96. """

97. Surcharge de l'opérateur `+` pour les vecteurs.

99. Args :

100. other(Vector): Un autre vecteur

102. Raises :

103. TypeError si other n'est pas un objet Vector

104. """

105. try:

106. return Vector(self.x + other.x, self.y + other.y, self.z + other.z)

107. except (ValueError, TypeError, AttributeError):

108. raise TypeError("Tried to add Vector and non-Vector objects")

110. def __sub__(self, other):

111. """

112. Surcharge de l'opérateur `-` pour les vecteurs.

114. Args :

115. other(Vector): Un autre vecteur

117. Raises :

118. TypeError si other n'est pas un objet Vector

119. """

120. try:

121. return Vector(self.x - other.x, self.y - other.y, self.z - other.z)

122. except (ValueError, TypeError, AttributeError):

123. raise TypeError("Tried to substract Vector and non-Vector objects")

125. def __mul__(self, number):

126. """

127. Surcharge de l'opérateur `*` pour la multiplication entre

128. un vecteur et un nombre.

130. Args :

131. number(float): Un nombre à multiplier par le Vector.

133. Raises :

134. TypeError si other n'est pas un nombre

135. """

136. try:

137. return Vector(number * self.x, number * self.y, number * self.z)

138. except (ValueError, TypeError, AttributeError):

139. raise TypeError("Tried to multiply Vector and non-number objects")

141. __rmul__ = __mul__  # Ligne pour autoriser la multiplication à droite

143. def norm(self):

144. """

145. Calcul de la norme 2 d'un vecteur.

147. Returns :

148. sqrt(x**2 + y**2 + z**2)

149. """

150. return (self.x ** 2 + self.y ** 2 + self.z ** 2) ** (1 / 2)

152. def clone(self):

153. """

154. Méthode pour construire un nouveau Vecteur, copie de self.

155. """

156. return Vector(self.x, self.y, self.z)

159. ###############################################

160. ##### Quelques test pour la classe Vector #####

161. ###############################################

163. def test_VectorInit():

164. with pytest.raises(TypeError):

165. vec = Vector('Test', 'avec', 'strings')

166. vec = Vector(Vector())

167. vec = Vector([])

168. vec = Vector(0, -53.76, math.pi)

169. assert vec.x == 0

170. assert vec.y == -53.76

171. assert vec.z == math.pi

174. def test_VectorStr():

175. vec = Vector(0, 600, -2)

176. assert str(vec) == '(0.00,600.00,-2.00)'

179. def test_VectorEq():  # teste aussi l'opérateur !=

180. vec = Vector(2, 3, -5)

181. vec2 = Vector(2, 3, -4)

182. assert vec != vec2

183. assert vec != Vector(0, 3, -5)

184. with pytest.raises(TypeError):

185. Vector(2, 1, 4) == "Testing strings"

186. Vector(2, 1, 4) == 42

187. Vector(2, 1, 4) == ['list']

190. def test_VectorAdd():

191. vec = Vector(2, 3, -5)

192. vec2 = Vector(2, -50, 5)

193. assert (vec + vec2) == Vector(4, -47, 0)

196. def test_VectorSub():

197. vec = Vector(1, -7, 9)

198. vec2 = Vector()

199. assert (vec - vec) == Vector()

200. assert (vec - vec2) == vec

203. def test_VectorMul():

204. vec = Vector(1, -7, 9) * 2

205. vec2 = 6 * Vector(1, -1, 2)

206. assert vec == Vector(2, -14, 18)

207. assert vec2 == Vector(6, -6, 12)

210. def test_VectorNorm():

211. assert Vector().norm() == 0

212. assert Vector(1, 0, 0).norm() == 1

213. assert Vector(2, -5, -4).norm() == 45 ** (1 / 2)

216. def test_VectorClone():

217. vec = Vector(3, 2, 9)

218. vec2 = vec.clone()

219. assert vec == vec2

220. vec2.x = 1

221. assert vec != vec2

224. ############################################################

225. ##### Une classe point matériel qui se gère en interne #####

226. ############################################################

228. class Particle(object):

230. """

231. La classe Particle représente un point matériel doté d'une masse,

232. d'une position et d'une vitesse. Elle possède également une méthode

233. pour calculer la force gravitationnelle exercée par une autre particule.

234. Enfin, la méthode update lui permet de mettre à jour sa position et

235. sa vitesse en fonction des forces subies.

236. """

238. def __init__(self, mass=1, position=Vector(), speed=Vector()):

239. """

240. Le constructeur de la classe Particle.

241. Définit un point matériel avec une position et une vitesse initiales.

243. Args :

244. mass(float): La masse de la particule (doit être

245. strictement positive)

246. position(Vector): La position initiale de la particule

247. speed(Vector): La vitesse initiale de la particule

249. Raises :

250. TypeError si la masse n'est pas un nombre, ou si la position ou

251. la vitesse ne sont pas des Vector

252. ValueError si la masse est négative ou nulle

253. """

254. try:

255. self.mass = float(mass)

256. self.position = position.clone()

257. self.speed = speed.clone()

258. except (ValueError, TypeError, AttributeError):

259. raise TypeError("The mass must be a positive float number. "

260. "The position and speed must Vector objects.")

261. try:

262. assert mass > 0  # une masse négative ou nulle pose des problèmes

263. except AssertionError:

264. raise ValueError("The mass must be strictly positive")

265. self.force = Vector()

267. def __str__(self):

268. """

269. Surcharge de l'opérateur `str`.

271. Returns :

272. "Particle with mass m, position (x,y,z) and speed (vx,vy,vz)"

273. with 2 decimals

274. """

275. return "Particle with mass {:.2f}, position {} " \

276. "and speed {}".format(self.mass, self.position, self.speed)

278. def computeForce(self, other):

279. """

280. Calcule la force gravitationnelle exercée par une Particule

281. other sur self.

283. Args :

284. other(Particle): Une autre particule, source de l'interaction

286. Raises :

287. TypeError si other n'est pas un objet Vector

288. """

289. try:

290. r = self.position - other.position

291. self.force = -self.mass * other.mass / r.norm() ** 3 * r

292. except AttributeError:

293. raise TypeError("Tried to compute the force created by "

294. "a non-Particle object")

296. def update(self, dt):

297. """

298. Mise à jour de la position et la vitesse au cours du temps.

300. Args :

301. dt(float): Pas de temps d'intégration.

302. """

303. try:

304. d = float(dt)

305. except (ValueError, TypeError, AttributeError):

306. raise TypeError("The integration timestep must be a number")

307. self.speed += self.force * dt * (1 / self.mass)

308. self.position += self.speed * dt

311. #############################################

312. ##### Des tests pour la classe Particle #####

313. #############################################

315. def test_ParticleInit():

316. with pytest.raises(TypeError):

317. p = Particle("blabla")

318. p = Particle(2, position='hum')  # on vérifie les erreurs sur Vector

319. p = Particle([])

320. p = Particle(3, Vector(2, 1, 4), Vector(-1, -1, -1))

321. assert p.mass == 3

322. assert p.position == Vector(2, 1, 4)

323. assert p.speed == Vector(-1, -1, -1)

324. assert p.force == Vector()

327. def test_ParticleStr():

328. p = Particle(3, Vector(1, 2, 3), Vector(-1, -2, -3))

329. assert str(p) == "Particle with mass 3.00, position (1.00,2.00,3.00) " \

330. "and speed (-1.00,-2.00,-3.00)"

333. def test_ParticleForce():

334. p = Particle(1, Vector(1, 0, 0))

335. p2 = Particle()

336. p.computeForce(p2)

337. assert p.force == Vector(-1, 0, 0)

338. p.position = Vector(2, -3, 6)

339. p.mass = 49

340. p.computeForce(p2)

341. assert p.force == Vector(-2 / 7, 3 / 7, -6 / 7)

344. def test_ParticleUpdate():

345. dt = 0.1

346. p = Particle(1, Vector(1, 0, 0), Vector())

347. p.computeForce(Particle())

348. p.update(dt)

349. assert p.speed == Vector(-0.1, 0, 0)

350. assert p.position == Vector(0.99, 0, 0)

353. #######################################################

354. ##### Une classe Ion qui hérite de point matériel #####

355. #######################################################

357. class Ion (Particle):

359. """

360. Un Ion est une particule ayant une charge en plus de sa masse et

361. intéragissant électrostatiquement plutôt que gravitationnellement.

362. La méthode computeForce remplace donc le calcul de la force

363. gravitationnelle de Newton par celui de la force

364. électrostatique de Coulomb.

365. """

367. def __init__(self, mass=1, charge=1, position=Vector(), speed=Vector()):

368. """

369. Le constructeur de la classe Ion.

371. Args :

372. mass(float): La masse de l'ion (doit être strictement positive)

373. charge(float): La charge de l'ion (doit être entière et

374. strictement positive)

375. position(Vector): La position initiale de la particule

376. speed(Vector): La vitesse initiale de la particule

378. Raises :

379. ValueError si charge < 0

380. TypeError si la masse n'est pas un réel,

381. si la charge n'est pas un entier,

382. si position ou speed ne sont pas des Vector

383. """

384. Particle.__init__(self, mass, position, speed)

385. try:

386. self.charge = int(charge)

387. except (ValueError, AttributeError, TypeError):

388. raise TypeError("The charge must be an integer.")

389. try:

390. assert self.charge > 0

391. except AssertionError:

392. raise ValueError("The charge must be positive.")

394. def __str__(self):

395. """

396. Surcharge de l'opérateur `str`.

398. Returns :

399. "Ion with mass m, charge q, position (x,y,z)

400. and speed (vx,vy,vz)" avec q entier et le reste à 2 décimales

401. """

402. return "Ion with mass {:.2f}, charge {:d}, position {} " \

403. "and speed {}".format(self.mass, self.charge,

404. self.position, self.speed)

406. def computeForce(self, other):

407. """

408. Calcule la force électrostatique de Coulomb exercée par un Ion other

409. sur self. Masque la méthode de Particle.

411. Args :

412. other(Ion): Un autre Ion, source de l'interaction.

413. Raises :

414. TypeError si other n'est pas un objet Ion

415. """

416. try:

417. r = self.position - other.position

418. self.force = self.charge * other.charge / r.norm() ** 3 * r

419. except (AttributeError, TypeError, ValueError):

420. raise TypeError("Tried to compute the force created by "

421. "a non-Ion object")

424. #######################################

425. ##### Des test pour la classe Ion #####

426. #######################################

428. def test_IonInit():

429. with pytest.raises(TypeError):

430. ion = Ion("blabla")

431. ion = Ion(2, position='hum')  # on vérifie une erreur sur Vector

432. ion = Ion(2, 'hum')           # on vérifie une erreur sur la charge

433. ion = Ion(2, 3, Vector(2, 1, 4), Vector(-1, -1, -1))

434. assert ion.mass == 2

435. assert ion.charge == 3

436. assert ion.position == Vector(2, 1, 4)

437. assert ion.speed == Vector(-1, -1, -1)

438. assert ion.force == Vector()

441. def test_IonStr():

442. ion = Ion(3, 2, Vector(1, 2, 3), Vector(-1, -2, -3))

443. assert str(ion) == "Ion with mass 3.00, charge 2, " \

444. "position (1.00,2.00,3.00) and speed (-1.00,-2.00,-3.00)"

447. def test_IonForce():

448. ion = Ion(mass=1, charge=1, position=Vector(1, 0, 0))

449. ion2 = Ion(charge=3)

450. ion.computeForce(ion2)

451. assert ion.force == Vector(3, 0, 0)

452. ion = Ion(charge=49, position=Vector(2, -3, 6))

453. ion.computeForce(ion2)

454. assert ion.force == Vector(6 / 7, -9 / 7, 18 / 7)

457. ###########################

458. ##### Un main de test #####

459. ###########################

461. if __name__ == '__main__':

463. # On lance tous les tests en bloc pour commencer

464. print " Test functions ".center(50, "*")

465. print "Testing Vector class...",

466. test_VectorInit()

467. test_VectorStr()

468. test_VectorEq()

469. test_VectorAdd()

470. test_VectorSub()

471. test_VectorMul()

472. test_VectorNorm()

473. test_VectorClone()

474. print "ok"

475. print "Testing Particle class...",

476. test_ParticleInit()

477. test_ParticleStr()

478. test_ParticleForce()

479. test_ParticleUpdate()

480. print "ok"

481. print "Testing Ion class...",

482. test_IonInit()

483. test_IonStr()

484. test_IonForce()

485. print "ok"

486. print " Test end ".center(50, "*"), "\n"

488. # Un petit calcul physique

489. print " Physical computations ".center(50, "*")

490. dt = 0.0001

492. # problème à force centrale gravitationnelles, cas circulaire

493. ntimesteps = int(10000 * math.pi)  # durée pour parcourir pi

494. center = Particle()

495. M = Particle(mass=1, position=Vector(1, 0, 0), speed=Vector(0, 1, 0))

496. print "** Gravitationnal computation of central-force motion for a {}" \

497. .format(str(M))

498. for i in range(ntimesteps):

499. M.computeForce(center)

500. M.update(dt)

501. print "\t => Final system : {}".format(str(M))

503. # problème à force centrale électrostatique, cas rectiligne

504. center = Ion()

505. M = Ion(charge=4, position=Vector(0, 0, 1), speed=Vector(0, 0, -1))

506. print "** Electrostatic computation of central-force motion for a {}" \

507. .format(str(M))

508. for i in range(ntimesteps):

509. M.computeForce(center)

510. M.update(dt)

511. print "\t => Final system : {}".format(str(M))

513. print " Physical computations end ".center(50, "*")

Source : particleSol.py

X-K - Jeu de la vie

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2017-07-09 19:14 ycopin@lyonovae03.in2p3.fr>

5. import random

8. class Life(object):

10. cells = {False: ".", True: "#"}  # Dead and living cell representations

12. def __init__(self, h, w, periodic=False):

13. """

14. Create a 2D-list (the game grid *G*) with the wanted size (*h*

15. rows, *w* columns) and initialize it with random booleans

16. (dead/alive). The world is periodic if *periodic* is True.

17. """

19. self.h = int(h)

20. self.w = int(w)

21. assert self.h > 0 and self.w > 0

22. # Random initialization of a h×w world

23. self.world = [[random.choice([True, False])

24. for j in range(self.w)]

25. for i in range(self.h)]  # h rows of w elements

26. self.periodic = periodic

28. def get(self, i, j):

29. """

30. This method returns the state of cell (*i*,*j*) safely, even

31. if the (*i*,*j*) is outside the grid.

32. """

34. if self.periodic:

35. return self.world[i % self.h][j % self.w]  # Periodic conditions

36. else:

37. if (0 <= i < self.h) and (0 <= j < self.w):  # Inside grid

38. return self.world[i][j]

39. else:                       # Outside grid

40. return False            # There's nobody out there...

42. def __str__(self):

43. """

44. Convert the grid to a visually handy string.

45. """

47. return '\n'.join([''.join([self.cells[val] for val in row])

48. for row in self.world])

50. def evolve_cell(self, i, j):

51. """

52. Tells if cell (*i*,*j*) will survive during game iteration,

53. depending on the number of living neighboring cells.

54. """

56. alive = self.get(i, j)           # Current cell status

57. # Count living cells *around* current one (excluding current one)

58. count = sum(self.get(i + ii, j + jj)

59. for ii in [-1, 0, 1]

60. for jj in [-1, 0, 1]

61. if (ii, jj) != (0, 0))

63. if count == 3:

64. # A cell w/ 3 neighbors will either stay alive or resuscitate

65. future = True

66. elif count < 2 or count > 3:

67. # A cell w/ too few or too many neighbors will die

68. future = False

69. else:

70. # A cell w/ 2 or 3 neighbors will stay as it is (dead or alive)

71. future = alive             # Current status

73. return future

75. def evolve(self):

76. """

77. Evolve the game grid by one step.

78. """

80. # Update the grid

81. self.world = [[self.evolve_cell(i, j)

82. for j in range(self.w)]

83. for i in range(self.h)]

85. if __name__ == "__main__":

87. import time

89. h, w = (20, 60)                       # (nrows,ncolumns)

91. # Instanciation (including initialization)

92. life = Life(h, w, periodic=True)

94. while True:                         # Infinite loop! (Ctrl-C to break)

95. print life                      # Print current world

96. print "\n"

97. time.sleep(0.1)                 # Pause a bit

98. life.evolve()                   # Evolve world

Source : life.py

X-L - Median Absolute Deviation

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. import numpy as N

7. def mad(a, axis=None):

8. """

9. Compute *Median Absolute Deviation* of an array along given axis.

10. """

12. # Median along given axis, but *keeping* the reduced axis so that

13. # result can still broadcast against a.

14. med = N.median(a, axis=axis, keepdims=True)

15. mad = N.median(N.absolute(a - med), axis=axis)  # MAD along given axis

17. return mad

19. if __name__ == '__main__':

21. x = N.arange(4 * 5, dtype=float).reshape(4, 5)

23. print "x =\n", x

24. print "MAD(x, axis=None)   =", mad(x)

25. print "MAD(x, axis=0)      =", mad(x, axis=0)

26. print "MAD(x, axis=1)      =", mad(x, axis=1)

27. print "MAD(x, axis=(0, 1)) =", mad(x, axis=(0, 1))

Source : mad.py

X-M - Distribution du pull

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. import numpy as N

7. def pull(x, dx):

8. """

9. Compute the pull from x, dx.

11. * Input data: x = [x_i], error dx = [s_i] Optimal

12. * (variance-weighted) mean: E = sum(x_i/s_i²)/sum(1/s_i²) Variance

13. * on weighted mean: var(E) = 1/sum(1/s_i²) Pull: p_i = (x_i -

14. * E_i)/sqrt(var(E_i) + s_i²) where E_i and var(E_i) are computed

15. * without point i.

17. If errors s_i are correct, the pull distribution is centered on 0

18. with standard deviation of 1.

19. """

21. assert x.ndim == dx.ndim == 1, "pull works on 1D-arrays only."

22. assert len(x) == len(dx), "dx should be the same length as x."

24. n = len(x)

26. i = N.resize(N.arange(n), n * n)  # 0,...,n-1,0,...n-1,..., n times (n²,)

27. i[::n + 1] = -1                  # Mark successively 0,1,2,...,n-1

28. # Remove marked indices & reshape (n,n-1)

29. j = i[i >= 0].reshape((n, n - 1))

31. v = dx ** 2                      # Variance

32. w = 1 / v                        # Variance (optimal) weighting

34. Ei = N.average(x[j], weights=w[j], axis=1)  # Weighted mean (n,)

35. vEi = 1 / N.sum(w[j], axis=1)                # Variance of mean (n,)

37. p = (x - Ei) / N.sqrt(vEi + v)               # Pull (n,)

39. return p

41. if __name__ == '__main__':

43. import matplotlib.pyplot as P

44. import scipy.stats as SS

46. n = 1000

47. mu = 1.

48. sig = 2.

50. # Normally distributed random sample of size n, with mean=mu and std=sig

51. x = N.random.normal(loc=mu, scale=sig, size=n)

52. dx = N.full_like(x, sig)              # Formal (true) errors

54. p = pull(x, dx)                       # Pull computation

56. m, s = p.mean(), p.std(ddof=1)

57. print "Pull ({} entries): mean={:.2f}, std={:.2f}".format(n, m, s)

59. fig, ax = P.subplots()

60. _, bins, _ = ax.hist(p, bins='auto', normed=True,

61. histtype='stepfilled',

62. label=u"#={}, µ={:.3f}, σ={:.3f}".format(n, m, s))

63. y = N.linspace(-3, 3)

64. ax.plot(y, SS.norm.pdf(y), label=ur"$\mathcal{N}$(µ=0, σ²=1)")

65. ax.set(title='Pull distribution', xlabel='Pull')

66. ax.legend(loc='upper left')

68. P.show()

Source : pull.py

X-N - Calculs numériques (numerique.ipynb)

numerique.ipynb

In [1]

import numpy as N
import matplotlib.pyplot as P
# Insert figures within notebook
%matplotlib inline

X-N-1 - Quadrature

Calcul de l'intégrale

In [3] :import scipy.integrate as SI

In [3]

def integrand(x):

return x**3 / (N.exp(x) - 1)

In [4]

q, dq = SI.quad(integrand, 0, N.inf)
print "Intégrale:", q
print "Erreur estimée:", dq
print "Erreur absolue:", (q - (N.pi ** 4 / 15))

Intégrale : 6.49393940227

Erreur estimée : 2.62847076684e-09

Erreur absolue : 1.7763568394e-15

/data/ycopin/Softs/local/lib/python2.7/site-packages/ipykernel/__main__.py:3: RuntimeWarning: overflow encountered in exp
app.launch_new_instance()

X-N-2 - Zéro d'une fonction

Résolution numérique de l'équation

In [5]

def func(x):

return x * N.exp(x) / (N.exp(x) - 1) - 5

Il faut d'abord déterminer un intervalle contenant la solution, c.-à-d. le zéro de func. Puisque et , il est intéressant de tracer l'allure de la courbe sur ce domaine :

In [6]

x = N.logspace(-1, 1)  # 50 points logarithmiquement espacés de 10**-1 = 0.1 à 10**1 = 10
P.plot(x, func(x));

In [7] :import scipy.optimize as SO

In [8]

zero = SO.brentq(func, 1, 10)
print "Solution:", zero

Solution : 4.96511423174

X-O - Quartet d'Anscombe

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

4. import numpy as N

5. import scipy.stats as SS

6. import matplotlib.pyplot as P

9. def printStats(x, y):

10. """

11. Print out means and variances for x and y, as well as correlation

12. coeff. (Pearson) and linear regression for y vs. x.

13. """

15. assert N.shape(x) == N.shape(y), "Incompatible input arrays"

17. print "x: mean={:.2f}, variance={:.2f}".format(N.mean(x), N.var(x))

18. print "y: mean={:.2f}, variance={:.2f}".format(N.mean(y), N.var(y))

19. print "y vs. x: corrcoeff={:.2f}".format(SS.pearsonr(x, y)[0])

20. # slope, intercept, r_value, p_value, std_err

21. a, b, _, _, _ = SS.linregress(x, y)

22. print "y vs. x: y = {:.2f} x + {:.2f}".format(a, b)

25. def plotStats(ax, x, y, title='', fancy=True):

26. """

27. Plot y vs. x, and linear regression.

28. """

30. assert N.shape(x) == N.shape(y), "Incompatible input arrays"

32. # slope, intercept, r_value, p_value, std_err

33. a, b, r, _, _ = SS.linregress(x, y)

35. # Data + corrcoeff

36. ax.plot(x, y, 'bo', label="r = {:.2f}".format(r))

38. # Linear regression

39. xx = N.array([0, 20])

40. yy = a * xx + b

41. ax.plot(xx, yy, 'r-', label="y = {:.2f} x + {:.2f}".format(a, b))

43. leg = ax.legend(loc='upper left', fontsize='small')

45. if fancy:                   # Additional stuff

46. # Add mean line ± stddev

47. m = N.mean(x)

48. s = N.std(x, ddof=1)

49. ax.axvline(m, color='g', ls='--', label='_')  # Mean

50. ax.axvspan(m - s, m + s, color='g', alpha=0.2, label='_')  # Std-dev

52. m = N.mean(y)

53. s = N.std(y, ddof=1)

54. ax.axhline(m, color='g', ls='--', label='_')  # Mean

55. ax.axhspan(m - s, m + s, color='g', alpha=0.2, label='_')  # Std-dev

57. # Title and labels

58. ax.set_title(title)

59. if ax.is_last_row():

60. ax.set_xlabel("x")

61. if ax.is_first_col():

62. ax.set_ylabel("y")

65. if __name__ == '__main__':

67. quartet = N.genfromtxt("anscombe.dat")  # Read Anscombe's Quartet

69. fig = P.figure()

71. for i in range(4):                  # Loop over quartet sets x,y

72. ax = fig.add_subplot(2, 2, i + 1)

73. print "Dataset #{} ".format(i + 1) + "=" * 20

74. x, y = quartet[:, 2 * i:2 * i + 2].T

75. printStats(x, y)                # Print main statistics

76. plotStats(ax, x, y, title=str(i + 1))  # Plots

78. fig.suptitle("Anscombe's Quartet", fontsize='x-large')

80. P.show()

$ python anscombe.py

Dataset 1 ====================
x: mean=9.00, variance=10.00
y: mean=7.50, variance=3.75
y vs. x: corrcoeff=0.82
y vs. x: y = 0.50 x + 3.00
Dataset 2 ====================
x: mean=9.00, variance=10.00
y: mean=7.50, variance=3.75
y vs. x: corrcoeff=0.82
y vs. x: y = 0.50 x + 3.00
Dataset 3 ====================
x: mean=9.00, variance=10.00
y: mean=7.50, variance=3.75
y vs. x: corrcoeff=0.82
y vs. x: y = 0.50 x + 3.00
Dataset 4 ====================
x: mean=9.00, variance=10.00
y: mean=7.50, variance=3.75
y vs. x: corrcoeff=0.82
y vs. x: y = 0.50 x + 3.00

Source : anscombe.py

X-P - Suite logistique

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2012-09-05 02:37 ycopin@lyopc469>

5. import numpy as np

6. import random

7. import matplotlib.pyplot as plt

10. def iteration(r, niter=100):

12. x = random.uniform(0, 1)

13. i = 0

14. while i < niter and x < 1:

15. x = r * x * (1 - x)

16. i += 1

18. return x if x < 1 else -1

21. def generate_diagram(r, ntrials=50):

22. """

23. Cette fonction retourne (jusqu'à) *ntrials* valeurs d'équilibre

24. pour les *r* d'entrée.  Elle renvoie un tuple :

26. + le premier élément est la liste des valeurs prises par le paramètre *r*

27. + le second est la liste des points d'équilibre correspondants

28. """

30. r_v = []

31. x_v = []

32. for rr in r:

33. j = 0

34. while j < ntrials:

35. xx = iteration(rr)

36. if xx > 0:  # Convergence : il s'agit d'une valeur d'équilibre

37. r_v.append(rr)

38. x_v.append(xx)

39. j += 1                      # Nouvel essai

41. return r_v, x_v

43. r = np.linspace(0, 4, 1000)

44. x, y = generate_diagram(r)

46. plt.plot(x, y, 'r,')

47. plt.xlabel('r')

48. plt.ylabel('x')

49. plt.show()

Source : logistique.py

X-Q - Ensemble de Julia

1. #!/usr/bin/env python

2. # -*- coding: utf-8 -*-

3. # Time-stamp: <2012-09-05 02:10 ycopin@lyopc469>

5. """Visualisation de l'ensemble de julia

6. <http://fr.wikipedia.org/wiki/Ensemble_de_Julia>`_.

8. Exercice: proposer des solutions pour accélérer le calcul.

9. """

11. import numpy as np

12. import matplotlib.pyplot as plt

14. c = complex(0.284, 0.0122)           # Constante

16. xlim = 1.5                          # [-xlim,xlim] × i[-xlim,xlim]

17. nx = 1000                           # Nb de pixels

18. niter = 100                         # Nb d'itérations

20. x = np.linspace(-xlim, xlim, nx)    # nx valeurs de -xlim à +xlim

21. xx, yy = np.meshgrid(x, x)           # Tableaux 2D

22. z = xx + 1j * yy                      # Portion du plan complexe

23. for i in range(niter):              # Itération: z(n+1) = z(n)**2 + c

24. z = z ** 2 + c

26. # Visualisation

27. plt.imshow(np.abs(z), extent=[-xlim, xlim, -xlim, xlim], aspect='equal')

28. plt.title(c)

29. plt.show()

Source : julia.py

X-R - Trajectoire d'un boulet de canon (canon.ipynb)

Nous allons intégrer les équations du mouvement pour un boulet de canon soumis à des forces de frottement « turbulentes » (non linéaires) :

Cette équation différentielle non linéaire du 2d ordre doit être réécrite sous la forme de deux équations différentielles couplées du 1er ordre:

Il s'agit donc de résoudre une seule équation différentielle du 1er ordre en

In [1] :

1. %matplotlib inline

3. import numpy as N

4. import scipy.integrate as SI

5. import matplotlib.pyplot as P

Valeurs numériques pour un boulet de canon de 36 livres :

In [2] :

1. g = 9.81 # Pesanteur [m/s2]

2. cx = 0.45 # Coefficient de frottement d'une sphère

3. rhoAir = 1.2 # Masse volumique de l'air [kg/m3] au niveau de la mer, T=20°C

4. rad = 0.1748/2 # Rayon du boulet [m]

5. rho = 6.23e3 # Masse volumique du boulet [kg/m3]

6. mass = 4./3.*N.pi*rad**3 * rho # Masse du boulet [kg]

7. alpha = 0.5*cx*rhoAir*N.pi*rad**2 / mass # Coefficient de frottement par unité de masse

8. print "Masse du boulet: {:.2f} kg".format(mass)

9. print "Coefficient de frottement par unité de masse: {} S.I.".format(alpha)

Out [2] :

Masse du boulet : 17.42 kg
Coefficient de frottement par unité de masse : 0.000371899460424 S.I.

Conditions initiales :

In [3] :

1. v0 = 450. # Vitesse initiale [m/s]

2. alt = 45. # Inclinaison du canon [deg]

3. alt *= N.pi/180. # Inclinaison [rad]

4. z0 = (0.,0.,v0*N.cos(alt),v0*N.sin(alt)) # (x0,y0,vx0,vy0)

Temps caractéristique du système : (durée du régime transitoire). L'intégration des équations se fera sur un temps caractéristique, avec des pas de temps significativement plus petits.

In [4] :

1. tc = N.sqrt(mass/(g * alpha))

2. print "Temps caractéristique: {:.1f} s".format(tc)

3. t = N.linspace(0, tc, 100)

Temps caractéristique : 69.1 s

Définition de la fonction avec

In [5] :

1. def zdot(z, t):

2. """Calcul de la dérivée de z=(x,y,vx,vy) à l'instant t."""

4. x,y,vx,vy = z

5. alphav = alpha * N.hypot(vx, vy)

7. return (vx,vy,-alphav*vx,-g-alphav*vy) # dz/dt = (vx,vy,x..,y..)

Intégration numérique des équations du mouvement à l'aide de la fonction scipy.integrate.odeint :

In [6] : zs = SI.odeint(zdot, z0, t)

Le tableau zs contient les valeurs de à chaque instant : il est donc de taille (len(t),4).

In [7] :

1. ypos = zs[:,1]>=0 # y>0?

2. print "temps de coll. t(y~0) = {:.0f} s".format(t[ypos][-1]) # Dernier instant pour lequel y>0

3. print "portée x(y~0) = {:.0f} m".format(zs[ypos,0][-1]) # Portée approximative du canon

4. #print "y(y~0) = {:.0f} m".format(zs[ypos,1][-1]) # ~0

5. print "vitesse(y~0): {:.0f} m/s".format(N.hypot(zs[ypos,2][-1],zs[ypos,3][-1]))

temps de coll. t(y~0) = 36 s
portée x(y~0) = 3966 m
vitesse(y~0): 140 m/s

In [8] :

1. fig,ax = P.subplots()

2. ax.plot(zs[ypos,0],zs[ypos,1], 'b.-')

3. ax.set_xlabel("x [m]")

4. ax.set_ylabel("y [m]")

5. ax.set_title("Trajectoire du boulet")

Out[8] : <matplotlib.text.Text at 0x7fe8d8b06950>

canon.ipynb

XI - Note de la rédaction de Developpez.com

Nous tenons à remercier Yannick Copin qui nous a aimablement autorisés à publier son tutoriel : Analyse scientifique avec Python. Nous remercions également Claude Leloup pour la relecture orthographique.

L'auteur, Yannick Copin, est maître de conférences au sein de l'équipe Cosmologie observationnelle/EUCLID de l'Institut de Physique Nucléaire de Lyon de l'Université Claude Bernard (Lyon 1).

(1)ou from __future__ import braces :-)

(2)La notion de « portée » est plus complexe, je simplifie…

(3)Par opposition au LBYL du C/C++, basé sur une série exhaustive de tests a priori.

(4)De nombreuses distributions Linux sont encore basées sur Python 2.7 par défaut.

(5)Python dispose également d'une bibliothèque de calcul formel, sympy, et d'un environnement de calcul mathématique, sage.

(6)Les structures pandas.Panel (de rang 3), pandas.Panel4D (de rang 4) et pandas.PanelND (de rang arbitraire) sont considérées comme dépréciées et seront retirées dans une version ultérieure. Utiliser une indexation hiérarchique ou xarray.

(7)Dont seulement 19 ont été écrits.

(8)« Don't reinvent the wheel, unless you plan on learning more about wheels » - Jeff Atwood

(9)« Premature optimization is the root of all evil » - Donald Knuth

(10)pytest ne fait pas partie de la bibliothèque standard. Il vous faudra donc l'installer indépendamment si vous voulez l'utiliser.

(11)Et maintenant du code Windows !