Analyse scientifique avec Python

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2014-09-17 19:30:18 ycopin>

# Division réelle de type Python 3 - ligne admise
from __future__ import division

"""
Définition de la fonction carré et calcul de son intégrale entre 0
et 1 dans le main par la méthode des rectangles (subdivision en 100
pas)
"""

__author__ = "Adrien Licari <adrien.licari@ens-lyon.fr>"

# Définition de la fonction sq, admise au stade du TD 1


def sq(x):
    return x * x

# Début du programme principal (main)

# On définit les bornes d'intégration a et b, le nombre de pas n
a = 0
b = 1
n = 100

# Largeur des rectangles dx
h = (b - a) / n      # Division réelle!

integ = 0         # Cette variable accumulera les aires des rectangles

# On sait déjà que l'on va calculer n aires de rectangles, donc une
# boucle for est appropriée
for i in range(n):            # Boucle de 0 à n-1
    x = a + (i + 0.5) * h         # Abscisse du rectangle
    integ += sq(x) * h        # On ajoute au total l'aire du rectangle

print "Intégrale de x**2 entre a =", a, "et b =", b, "avec n =", n
# On affiche le résultat numérique
print "Résultat numérique: ", integ
theorie = (b ** 3 - a ** 3) / 3                     # Résultat analytique
# On affiche le résultat analytique
print "Résultat analytique:", theorie
print "Erreur relative:", (integ / theorie - 1)  # On affiche l'erreur

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2017-05-28 21:42 ycopin@lyonovae03.in2p3.fr>

"""
Jeu du Fizz Buzz
"""

__author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"


for i in range(1, 100):                    # Entiers de 1 à 99
    if ((i % 3) == 0) and ((i % 5) == 0):  # Multiple de 3 et 5
        print 'Fizz Buzz!',                # Affichage sans retour à la ligne
    elif ((i % 3) == 0):                   # Multiple de 3 uniquement
        print 'Fizz!',
    elif ((i % 5) == 0):                   # Multiple de 5 uniquement
        print 'Buzz!',
    else:
        print i,

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2014-09-17 18:25:00 ycopin>

"""
Calcul du PGCD de deux entiers
"""

__author__ = "Adrien Licari <adrien.licari@ens-lyon.fr>"

# Entiers dont on calcule le PGCD
a = 306
b = 756

# Test usuels : les nombres sont bien positifs et a > b
assert a > 0
assert b > 0
if (a < b):
    a, b = b, a                   # Interversion

a0, b0 = a, b                     # On garde une copie des valeurs originales

# On boucle jusqu'à ce que le reste soit nul, d'où la boucle while. Il
# faut être sûr que l'algorithme converge dans tous les cas !
while (b != 0):
    # On remplace a par b et b par le reste de la division euclidienne
    # de a par b
    a, b = b, a % b

print 'Le PGCD de', a0, 'et', b0, 'vaut', a  # On affiche le résultat
# Vérifications
print a0 // a, 'x', a, '=', (a0 // a * a)  # a//b: division euclidienne
print b0 // a, 'x', a, '=', (b0 // a * a)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Crible d'Ératosthène.

Source: http://fr.wikibooks.org/wiki/Exemples_de_scripts_Python#Implémentation_du_crible_d'Ératosthène
"""

__author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"

# start-sys
# Gestion simplifiée d'un argument entier sur la ligne de commande
import sys

if sys.argv[1:]:  # Présence d'au moins un argument sur la ligne de commande
    try:
        n = int(sys.argv[1])  # Essayer de lire le 1er argument comme un entier
    except ValueError:
        raise ValueError("L'argument '{}' n'est pas un entier"
                         .format(sys.argv[1]))
else:                        # Pas d'argument sur la ligne de commande
    n = 101                  # Valeur par défaut
# end-sys

# Liste des entiers *potentiellement* premiers. Les nb non premiers
# seront étiquetés par 0 au fur et à mesure.
l = range(n + 1)                          # <0,...,n>, 0 n'est pas premier
l[1] = 0                                # 1 n'est pas premier

i = 2                                   # Entier à tester
while i ** 2 <= n:                        # Inutile de tester jusqu'à n
    if l[i] != 0:                       # Si i n'est pas étiqueté (=0)...
        # ...étiqueter tous les multiples de i
        l[2 * i::i] = [0] * len(l[2 * i::i])
    i += 1                              # Passer à l'entier à tester suivant

# Afficher la liste des entiers premiers (c-à-d non étiquetés)
print "Liste des entiers premiers <= {} :".format(n)
print [i for i in l if i != 0]

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2017-07-13 14:29 ycopin@lyonovae03.in2p3.fr>

"""
Création et affichage d'un carré magique d'ordre impair.
"""

__author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"

n = 5                                   # Ordre du carré magique

# On vérifie que l'ordre est bien impair
assert n % 2 == 1, "L'ordre {} n'est pas impair".format(n)

# Le carré sera stocké dans une liste de n listes de n entiers
# Initialisation du carré: liste de n listes de n zéros.
array = [[0 for j in range(n)] for i in range(n)]

# Initialisation de l'algorithme
i, j = n, (n + 1) / 2         # Indices de l'algo (1-indexation)
array[i - 1][j - 1] = 1       # Attention: python utilise une 0-indexation

# Boucle sur valeurs restant à inclure dans le carré (de 2 à n**2)
for k in range(2, n**2 + 1):
    # Test de la case i+1, j+1 (modulo n)
    i2 = (i + 1) % n
    j2 = (j + 1) % n
    if array[i2 - 1][j2 - 1] == 0:  # La case est vide: l'utiliser
        i, j = i2, j2
    # La case est déjà remplie : utiliser la case i-1, j
    else:
        i = (i - 1) % n
    array[i - 1][j - 1] = k       # Remplissage de la case

# Affichage, avec vérification des sommes par ligne et par colonne
print "Carré magique d'ordre {}:".format(n)
for row in array:
    print '  '.join('{:2d}'.format(k) for k in row), '=', sum(row)
print '  '.join('==' for k in row)
print '  '.join(str(sum(array[i][j] for i in range(n))) for j in range(n))
$ python carre.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2014-09-24 17:31:53 ycopin>

__author__ = "Adrien Licari <adrien.licari@ens-lyon.fr>; Yannick Copin <y.copin@ipnl.in2p3.fr>"


def suite_syracuse(n):
    """
    Retourne la suite de Syracuse pour l'entier n.

    >>> suite_syracuse(15)
    [15, 46, 23, 70, 35, 106, 53, 160, 80, 40, 20, 10, 5, 16, 8, 4, 2, 1]
    """

    seq = [n]                      # La suite de Syracuse sera complétée...
    while (seq[-1] != 1):         # ...jusqu'à tomber sur 1
        if seq[-1] % 2 == 0:      # u_n est pair
            seq.append(seq[-1] // 2)  # Division euclidienne par 2
        else:                     # u_n est impair
            seq.append(seq[-1] * 3 + 1)

    return seq


def temps_syracuse(n, altitude=False):
    """
    Calcule le temps de vol (éventuellement en altitude) de la suite
    de Syracuse pour l'entier n.

    >>> temps_syracuse(15)
    17
    >>> temps_syracuse(15, altitude=True)
    10
    """

    seq = suite_syracuse(n)
    if not altitude:            # Temps de vol total
        return len(seq) - 1
    else:                       # Temps de vol en altitude
        # Construction de la séquence en altitude
        alt = []
        for i in seq:
            if i >= n:
                alt.append(i)
            else:
                break
        return len(alt) - 1

if __name__ == '__main__':

    n = 15
    print "Suite de Syracuse pour n =", n
    print suite_syracuse(n)
    print "Temps de vol total:      ", temps_syracuse(n)
    print "Temps de vol en altitude:", temps_syracuse(n, altitude=True)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division  # Pas de division euclidienne par défaut

"""
Tracé (via 'turtle') d'un flocon de Koch d'ordre arbitraire.

Dans le même genre:

- courbe de Peano (http://fr.wikipedia.org/wiki/Courbe_de_Peano)
- courbe de Hilbert (http://fr.wikipedia.org/wiki/Courbe_de_Hilbert)
- île de Gosper (http://fr.wikipedia.org/wiki/Île_de_Gosper)

Voir également:

- L-système: http://fr.wikipedia.org/wiki/L-système
- Autres exemples: http://natesoares.com/tutorials/python-fractals/
"""

__version__ = "Time-stamp: <2013-01-14 00:49 ycopin@lyopc469>"
__author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"

import turtle as T


def koch(niveau=3, iter=0, taille=100, delta=0):
    """
    Tracé du flocon de Koch de niveau 'niveau', de taille 'taille'
    (px).

    Cette fonction récursive permet d'initialiser le flocon (iter=0,
    par défaut), de tracer les branches fractales (0<iter<=niveau) ou
    bien juste de tracer un segment (iter>niveau).
    """

    if iter == 0:                         # Dessine le triangle de niveau 0
        T.title("Flocon de Koch - niveau {}".format(niveau))
        koch(iter=1, niveau=niveau, taille=taille, delta=delta)
        T.right(120)
        koch(iter=1, niveau=niveau, taille=taille, delta=delta)
        T.right(120)
        koch(iter=1, niveau=niveau, taille=taille, delta=delta)
    elif iter <= niveau:                  # Trace une section _/\_ du flocon
        koch(iter=iter + 1, niveau=niveau, taille=taille, delta=delta)
        T.left(60 + delta)
        koch(iter=iter + 1, niveau=niveau, taille=taille, delta=delta)
        T.right(120 + 2 * delta)
        koch(iter=iter + 1, niveau=niveau, taille=taille, delta=delta)
        T.left(60 + delta)
        koch(iter=iter + 1, niveau=niveau, taille=taille, delta=delta)
    else:                               # Trace le segment de dernier niveau
        T.forward(taille / 3 ** (niveau + 1))

if __name__ == '__main__':

    # start-argparse
    # Exemple d'utilisation de la bibliothèque de gestion d'arguments 'argparse'
    import argparse

    desc = u"Tracé (via 'turtle') d'un flocon de Koch d'ordre arbitraire."

    # Définition des options
    parser = argparse.ArgumentParser(description=desc)
    parser.add_argument('ordre', nargs='?', type=int,
                        help="Ordre du flocon, >0 [%(default)s]",
                        default=3)
    parser.add_argument('-t', '--taille', type=int,
                        help="Taille de la figure, >0 [%(default)s px]",
                        default=500)
    parser.add_argument('-d', '--delta', type=float,
                        help="Delta [%(default)s deg]",
                        default=0.)
    parser.add_argument('-f', '--figure', type=str,
                        help="Nom de la figure de sortie (format EPS)")
    parser.add_argument('-T', '--turbo',
                        action="store_true", default=False,
                        help="Mode Turbo")

    # Déchiffrage des options et arguments
    args = parser.parse_args()

    # Quelques tests sur les args et options
    if not args.ordre > 0:
        parser.error("Ordre requis '{}' invalide".format(args.ordre))

    if not args.taille > 0:
        parser.error("La taille de la figure doit être positive")
    # end-argparse

    if args.turbo:
        T.hideturtle()
        T.speed(0)

    # Tracé du flocon de Koch de niveau 3
    koch(niveau=args.ordre, taille=args.taille, delta=args.delta)
    if args.figure:
        # Sauvegarde de l'image
        print "Sauvegarde de la figure dans '{}'".format(args.figure)
        T.getscreen().getcanvas().postscript(file=args.figure)

    T.exitonclick()

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import random

nmin, nmax = 1, 100
nsol = random.randint(nmin, nmax)

print "Vous devez deviner un nombre entre {} et {}.".format(nmin, nmax)
ncoups = 0
try:
    while True:
        proposition = raw_input("Votre proposition: ")
        try:
            ntest = int(proposition)
            if not nmin <= ntest <= nmax:
                raise ValueError("Proposition invalide")
        except ValueError:
            print "Votre proposition {!r} n'est pas un entier " \
                "compris entre {} et {}.".format(proposition, nmin, nmax)
            continue
        ncoups += 1
        if nsol > ntest:
            print "C'est plus."
        elif nsol < ntest:
            print "C'est moins."
        else:
            print "Vous avez trouvé en {} coup{}!".format(
                ncoups, 's' if ncoups > 1 else '')
            break  # Interrompt la boucle while
except (KeyboardInterrupt, EOFError):  # Interception de Ctrl-C et Ctrl-D
    print "\nVous abandonnez après seulement {} coup{}!".format(
        ncoups, 's' if ncoups > 1 else '')

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pytest                   # Module (non-- standard) de tests


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

    """
    Classe définissant un `Animal`, caractérisé par son nom et son
    poids.
    """

    def __init__(self, nom, masse):
        """
        Méthode d'instanciation à partir d'un nom (str) et d'un poids (float).
        """

        # Ici, convertir les paramètres pour être sûr qu'il ont le bon
        # type. On utilisera `str` et `float`
        self.nom = str(nom)
        self.masse = float(masse)

        self.vivant = True       # Les animaux sont vivants à l'instanciation
        self.empoisonne = False  # Animal empoisonné ?

    def __str__(self):
        """
        Surcharge de l'opérateur `str`: l'affichage *informel* de l'objet
        dans l'interpréteur, p.ex. `print self` sera résolu comme
        `self.__str__()`

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

        return "{a.nom} ({a.masse:.1f} kg)".format(a=self)

    def estVivant(self):
        """Méthode booléenne, vraie si l'animal est vivant."""

        return self.vivant

    def mourir(self):
        """Change l'état interne de l'objet (ne retourne rien)."""

        self.vivant = False

    def __cmp__(self, other):
        """
        Surcharge des opérateurs de comparaison, sur la base de la masse
        des animaux.
        """

        return cmp(self.masse, other.masse)

    def __call__(self, other):
        """
        Surcharge de l'opérateur '()' pour manger un autre animal (qui
        meurt s'il est vivant) et prendre du poids (mais pas plus que
        la masse de l'autre ou 10 % de son propre poids).  Attention aux
        animaux empoisonnés !

        L'instruction `self(other)` sera résolue comme
        `self.__call__(other).
        """

        other.mourir()
        poids = min(other.masse, self.masse * 0.1)
        self.masse += poids
        other.masse -= poids
        if other.empoisonne:
            self.mourir()


class Chien(Animal):

    """
    Un `Chien` hérite de `Animal` avec des méthodes additionnelles
    (p.ex. l'aboiement et l'odorat).
    """

    def __init__(self, nom, masse=20, odorat=0.5):
        """Définit un chien plus ou moins fin limier."""

        # Initialisation de la classe parente
        Animal.__init__(self, nom, masse)

        # Attribut propre à la classe dérivée
        self.odorat = float(odorat)

    def __str__(self):

        return "{a.nom} (Chien, {a.masse:.1f} kg)".format(a=self)

    def aboyer(self):
        """Une méthode bien spécifique aux chiens."""

        print("Ouaf ! Ouaf !")

    def estVivant(self):
        """Quand on vérifie qu'un chien est vivant, il aboie."""

        vivant = Animal.estVivant(self)

        if vivant:
            self.aboyer()

        return vivant

# start-tests


def test_empty_init():
    with pytest.raises(TypeError):
        Animal()


def test_wrong_init():
    with pytest.raises(ValueError):
        Animal('Youki', 'lalala')


def test_init():
    youki = Animal('Youki', 600)
    assert youki.masse == 600
    assert youki.vivant
    assert youki.estVivant()
    assert not youki.empoisonne
# end-tests


def test_str():
    youki = Animal('Youki', 600)
    assert str(youki) == 'Youki (600.0 kg)'


def test_mort():
    youki = Animal('Youki', 600)
    assert youki.estVivant()
    youki.mourir()
    assert not youki.estVivant()


def test_lt():
    medor = Animal('Medor', 600)
    kiki = Animal('Kiki', 20)
    assert kiki < medor
    with pytest.raises(AttributeError):
        medor < 1


def test_mange():
    medor = Animal('Medor', 600)
    kiki = Animal('Kiki', 20)
    medor(kiki)                 # Médor mange Kiki
    assert medor.estVivant()
    assert not kiki.estVivant()
    assert kiki.masse == 0
    assert medor.masse == 620
    kiki = Animal("Kiki Jr.", 20)
    kiki(medor)                 # Kiki Jr. mange Médor
    assert not medor.estVivant()
    assert kiki.estVivant()
    assert kiki.masse == 22
    assert medor.masse == 618   # Médor a perdu du poids en se faisant manger!


def test_init_chien():
    medor = Chien('Medor', 600)
    assert isinstance(medor, Animal)
    assert isinstance(medor, Chien)
    assert medor.odorat == 0.5
    assert str(medor) == 'Medor (Chien, 600.0 kg)'
    assert medor.estVivant()

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2014-09-04 16:56:47 alicari>


from __future__ import division  # division réelle de type python 3, admis
import pytest                    # pytest importé pour les tests unitaires
import math

"""
Définition d'une classe point matériel, avec sa masse,
sa position et sa vitesse, et des méthodes pour le déplacer.
Le main test applique cela à un problème à force centrale
gravitationnel ou électrostatique.
Remarque : toutes les unités ont été choisies adimensionnées.
"""

__author__ = "Adrien Licari <adrien.licari@ens-lyon.fr>"


# Un critère pour déterminer l'égalité entre réels
tolerance = 1e-8


#############################################################################
### Définition de la classe Vector, utile pour la position et la vitesse. ###
#############################################################################

class Vector(object):

    """
    Une classe structure simple contenant 3 coordonnées.
    Une méthode est disponible pour en calculer la norme et
    une surcharge des opérateurs ==, !=, +, - et * est proposée.
    """

    def __init__(self, x=0, y=0, z=0):
        """
        Constructeur de la classe vector.
        Par défaut, construit le vecteur nul.

        Args:
                x,y,z(float): Les composantes du vecteur à construire.

        Raises:
                TypeError en cas de composantes non réelles
        """
        try:
            self.x = float(x)
            self.y = float(y)
            self.z = float(z)
        except (ValueError, TypeError, AttributeError):
            raise TypeError("The given coordinates must be numbers")

    def __str__(self):
        """
        Surcharge de l'opérateur `str`.

        Returns :
                "(x,y,z)" with 2 decimals
        """
        return "({:.2f},{:.2f},{:.2f})".format(self.x, self.y, self.z)

    def __eq__(self, other):
        """
        Surcharge de l'opérateur `==` pour tester l'égalité
        entre deux vecteurs.

        Args :
                other(Vector): Un autre vecteur

        Raises :
                TypeError si other n'est pas un objet Vector
        """
        try:
            return abs(self.x - other.x) < tolerance and \
                abs(self.y - other.y) < tolerance and \
                abs(self.z - other.z) < tolerance
        except (ValueError, TypeError, AttributeError):
            raise TypeError("Tried to compare Vector and non-Vector objects")

    def __ne__(self, other):
        """
        Surcharge de l'opérateur `!=` pour tester l'inégalité
        entre deux vecteurs.

        Args :
                other(Vector): Un autre vecteur

        Raises :
                TypeError si other n'est pas un objet Vector
        """
        return not self == other

    def __add__(self, other):
        """
        Surcharge de l'opérateur `+` pour les vecteurs.

        Args :
                other(Vector): Un autre vecteur

        Raises :
                TypeError si other n'est pas un objet Vector
        """
        try:
            return Vector(self.x + other.x, self.y + other.y, self.z + other.z)
        except (ValueError, TypeError, AttributeError):
            raise TypeError("Tried to add Vector and non-Vector objects")

    def __sub__(self, other):
        """
        Surcharge de l'opérateur `-` pour les vecteurs.

        Args :
                other(Vector): Un autre vecteur

        Raises :
                TypeError si other n'est pas un objet Vector
        """
        try:
            return Vector(self.x - other.x, self.y - other.y, self.z - other.z)
        except (ValueError, TypeError, AttributeError):
            raise TypeError("Tried to substract Vector and non-Vector objects")

    def __mul__(self, number):
        """
        Surcharge de l'opérateur `*` pour la multiplication entre
        un vecteur et un nombre.

        Args :
                number(float): Un nombre à multiplier par le Vector.

        Raises :
                TypeError si other n'est pas un nombre
        """
        try:
            return Vector(number * self.x, number * self.y, number * self.z)
        except (ValueError, TypeError, AttributeError):
            raise TypeError("Tried to multiply Vector and non-number objects")

    __rmul__ = __mul__  # Ligne pour autoriser la multiplication à droite

    def norm(self):
        """
        Calcul de la norme 2 d'un vecteur.

        Returns :
                sqrt(x**2 + y**2 + z**2)
        """
        return (self.x ** 2 + self.y ** 2 + self.z ** 2) ** (1 / 2)

    def clone(self):
        """
        Méthode pour construire un nouveau Vecteur, copie de self.
        """
        return Vector(self.x, self.y, self.z)


###############################################
##### Quelques test pour la classe Vector #####
###############################################

def test_VectorInit():
    with pytest.raises(TypeError):
        vec = Vector('Test', 'avec', 'strings')
        vec = Vector(Vector())
        vec = Vector([])
    vec = Vector(0, -53.76, math.pi)
    assert vec.x == 0
    assert vec.y == -53.76
    assert vec.z == math.pi


def test_VectorStr():
    vec = Vector(0, 600, -2)
    assert str(vec) == '(0.00,600.00,-2.00)'


def test_VectorEq():  # teste aussi l'opérateur !=
    vec = Vector(2, 3, -5)
    vec2 = Vector(2, 3, -4)
    assert vec != vec2
    assert vec != Vector(0, 3, -5)
    with pytest.raises(TypeError):
        Vector(2, 1, 4) == "Testing strings"
        Vector(2, 1, 4) == 42
        Vector(2, 1, 4) == ['list']


def test_VectorAdd():
    vec = Vector(2, 3, -5)
    vec2 = Vector(2, -50, 5)
    assert (vec + vec2) == Vector(4, -47, 0)


def test_VectorSub():
    vec = Vector(1, -7, 9)
    vec2 = Vector()
    assert (vec - vec) == Vector()
    assert (vec - vec2) == vec


def test_VectorMul():
    vec = Vector(1, -7, 9) * 2
    vec2 = 6 * Vector(1, -1, 2)
    assert vec == Vector(2, -14, 18)
    assert vec2 == Vector(6, -6, 12)


def test_VectorNorm():
    assert Vector().norm() == 0
    assert Vector(1, 0, 0).norm() == 1
    assert Vector(2, -5, -4).norm() == 45 ** (1 / 2)


def test_VectorClone():
    vec = Vector(3, 2, 9)
    vec2 = vec.clone()
    assert vec == vec2
    vec2.x = 1
    assert vec != vec2


############################################################
##### Une classe point matériel qui se gère en interne #####
############################################################

class Particle(object):

    """
    La classe Particle représente un point matériel doté d'une masse,
    d'une position et d'une vitesse. Elle possède également une méthode
    pour calculer la force gravitationnelle exercée par une autre particule.
    Enfin, la méthode update lui permet de mettre à jour sa position et
    sa vitesse en fonction des forces subies.
    """

    def __init__(self, mass=1, position=Vector(), speed=Vector()):
        """
        Le constructeur de la classe Particle.
        Définit un point matériel avec une position et une vitesse initiales.

        Args :
                mass(float): La masse de la particule (doit être
                        strictement positive)
                position(Vector): La position initiale de la particule
                speed(Vector): La vitesse initiale de la particule

        Raises :
                TypeError si la masse n'est pas un nombre, ou si la position ou
                        la vitesse ne sont pas des Vector
                ValueError si la masse est négative ou nulle
        """
        try:
            self.mass = float(mass)
            self.position = position.clone()
            self.speed = speed.clone()
        except (ValueError, TypeError, AttributeError):
            raise TypeError("The mass must be a positive float number. "
                            "The position and speed must Vector objects.")
        try:
            assert mass > 0  # une masse négative ou nulle pose des problèmes
        except AssertionError:
            raise ValueError("The mass must be strictly positive")
        self.force = Vector()

    def __str__(self):
        """
        Surcharge de l'opérateur `str`.

        Returns :
                "Particle with mass m, position (x,y,z) and speed (vx,vy,vz)"
                        with 2 decimals
        """
        return "Particle with mass {:.2f}, position {} " \
            "and speed {}".format(self.mass, self.position, self.speed)

    def computeForce(self, other):
        """
        Calcule la force gravitationnelle exercée par une Particule
        other sur self.

        Args :
                other(Particle): Une autre particule, source de l'interaction

        Raises :
                TypeError si other n'est pas un objet Vector
        """
        try:
            r = self.position - other.position
            self.force = -self.mass * other.mass / r.norm() ** 3 * r
        except AttributeError:
            raise TypeError("Tried to compute the force created by "
                            "a non-Particle object")

    def update(self, dt):
        """
        Mise à jour de la position et la vitesse au cours du temps.

        Args :
                dt(float): Pas de temps d'intégration.
        """
        try:
            d = float(dt)
        except (ValueError, TypeError, AttributeError):
            raise TypeError("The integration timestep must be a number")
        self.speed += self.force * dt * (1 / self.mass)
        self.position += self.speed * dt


#############################################
##### Des tests pour la classe Particle #####
#############################################

def test_ParticleInit():
    with pytest.raises(TypeError):
        p = Particle("blabla")
        p = Particle(2, position='hum')  # on vérifie les erreurs sur Vector
        p = Particle([])
    p = Particle(3, Vector(2, 1, 4), Vector(-1, -1, -1))
    assert p.mass == 3
    assert p.position == Vector(2, 1, 4)
    assert p.speed == Vector(-1, -1, -1)
    assert p.force == Vector()


def test_ParticleStr():
    p = Particle(3, Vector(1, 2, 3), Vector(-1, -2, -3))
    assert str(p) == "Particle with mass 3.00, position (1.00,2.00,3.00) " \
        "and speed (-1.00,-2.00,-3.00)"


def test_ParticleForce():
    p = Particle(1, Vector(1, 0, 0))
    p2 = Particle()
    p.computeForce(p2)
    assert p.force == Vector(-1, 0, 0)
    p.position = Vector(2, -3, 6)
    p.mass = 49
    p.computeForce(p2)
    assert p.force == Vector(-2 / 7, 3 / 7, -6 / 7)


def test_ParticleUpdate():
    dt = 0.1
    p = Particle(1, Vector(1, 0, 0), Vector())
    p.computeForce(Particle())
    p.update(dt)
    assert p.speed == Vector(-0.1, 0, 0)
    assert p.position == Vector(0.99, 0, 0)


#######################################################
##### Une classe Ion qui hérite de point matériel #####
#######################################################

class Ion (Particle):

    """
    Un Ion est une particule ayant une charge en plus de sa masse et
    intéragissant électrostatiquement plutôt que gravitationnellement.
    La méthode computeForce remplace donc le calcul de la force
    gravitationnelle de Newton par celui de la force
    électrostatique de Coulomb.
    """

    def __init__(self, mass=1, charge=1, position=Vector(), speed=Vector()):
        """
        Le constructeur de la classe Ion.

        Args :
                mass(float): La masse de l'ion (doit être strictement positive)
                charge(float): La charge de l'ion (doit être entière et
                        strictement positive)
                position(Vector): La position initiale de la particule
                speed(Vector): La vitesse initiale de la particule

        Raises :
                ValueError si charge < 0
                TypeError si la masse n'est pas un réel,
                        si la charge n'est pas un entier,
                        si position ou speed ne sont pas des Vector
        """
        Particle.__init__(self, mass, position, speed)
        try:
            self.charge = int(charge)
        except (ValueError, AttributeError, TypeError):
            raise TypeError("The charge must be an integer.")
        try:
            assert self.charge > 0
        except AssertionError:
            raise ValueError("The charge must be positive.")

    def __str__(self):
        """
        Surcharge de l'opérateur `str`.

        Returns :
                "Ion with mass m, charge q, position (x,y,z)
                and speed (vx,vy,vz)" avec q entier et le reste à 2 décimales
        """
        return "Ion with mass {:.2f}, charge {:d}, position {} " \
            "and speed {}".format(self.mass, self.charge,
                                  self.position, self.speed)

    def computeForce(self, other):
        """
        Calcule la force électrostatique de Coulomb exercée par un Ion other
        sur self. Masque la méthode de Particle.

        Args :
                other(Ion): Un autre Ion, source de l'interaction.
        Raises :
                TypeError si other n'est pas un objet Ion
        """
        try:
            r = self.position - other.position
            self.force = self.charge * other.charge / r.norm() ** 3 * r
        except (AttributeError, TypeError, ValueError):
            raise TypeError("Tried to compute the force created by "
                            "a non-Ion object")


#######################################
##### Des test pour la classe Ion #####
#######################################

def test_IonInit():
    with pytest.raises(TypeError):
        ion = Ion("blabla")
        ion = Ion(2, position='hum')  # on vérifie une erreur sur Vector
        ion = Ion(2, 'hum')           # on vérifie une erreur sur la charge
    ion = Ion(2, 3, Vector(2, 1, 4), Vector(-1, -1, -1))
    assert ion.mass == 2
    assert ion.charge == 3
    assert ion.position == Vector(2, 1, 4)
    assert ion.speed == Vector(-1, -1, -1)
    assert ion.force == Vector()


def test_IonStr():
    ion = Ion(3, 2, Vector(1, 2, 3), Vector(-1, -2, -3))
    assert str(ion) == "Ion with mass 3.00, charge 2, " \
        "position (1.00,2.00,3.00) and speed (-1.00,-2.00,-3.00)"


def test_IonForce():
    ion = Ion(mass=1, charge=1, position=Vector(1, 0, 0))
    ion2 = Ion(charge=3)
    ion.computeForce(ion2)
    assert ion.force == Vector(3, 0, 0)
    ion = Ion(charge=49, position=Vector(2, -3, 6))
    ion.computeForce(ion2)
    assert ion.force == Vector(6 / 7, -9 / 7, 18 / 7)


###########################
##### Un main de test #####
###########################

if __name__ == '__main__':

    # On lance tous les tests en bloc pour commencer
    print " Test functions ".center(50, "*")
    print "Testing Vector class...",
    test_VectorInit()
    test_VectorStr()
    test_VectorEq()
    test_VectorAdd()
    test_VectorSub()
    test_VectorMul()
    test_VectorNorm()
    test_VectorClone()
    print "ok"
    print "Testing Particle class...",
    test_ParticleInit()
    test_ParticleStr()
    test_ParticleForce()
    test_ParticleUpdate()
    print "ok"
    print "Testing Ion class...",
    test_IonInit()
    test_IonStr()
    test_IonForce()
    print "ok"
    print " Test end ".center(50, "*"), "\n"

    # Un petit calcul physique
    print " Physical computations ".center(50, "*")
    dt = 0.0001

    # problème à force centrale gravitationnelles, cas circulaire
    ntimesteps = int(10000 * math.pi)  # durée pour parcourir pi
    center = Particle()
    M = Particle(mass=1, position=Vector(1, 0, 0), speed=Vector(0, 1, 0))
    print "** Gravitationnal computation of central-force motion for a {}" \
        .format(str(M))
    for i in range(ntimesteps):
        M.computeForce(center)
        M.update(dt)
    print "\t => Final system : {}".format(str(M))

    # problème à force centrale électrostatique, cas rectiligne
    center = Ion()
    M = Ion(charge=4, position=Vector(0, 0, 1), speed=Vector(0, 0, -1))
    print "** Electrostatic computation of central-force motion for a {}" \
        .format(str(M))
    for i in range(ntimesteps):
        M.computeForce(center)
        M.update(dt)
    print "\t => Final system : {}".format(str(M))

    print " Physical computations end ".center(50, "*")

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2017-07-09 19:14 ycopin@lyonovae03.in2p3.fr>

import random


class Life(object):

    cells = {False: ".", True: "#"}  # Dead and living cell representations

    def __init__(self, h, w, periodic=False):
        """
        Create a 2D-list (the game grid *G*) with the wanted size (*h*
        rows, *w* columns) and initialize it with random booleans
        (dead/alive). The world is periodic if *periodic* is True.
        """

        self.h = int(h)
        self.w = int(w)
        assert self.h > 0 and self.w > 0
        # Random initialization of a h×w world
        self.world = [[random.choice([True, False])
                       for j in range(self.w)]
                      for i in range(self.h)]  # h rows of w elements
        self.periodic = periodic

    def get(self, i, j):
        """
        This method returns the state of cell (*i*,*j*) safely, even
        if the (*i*,*j*) is outside the grid.
        """

        if self.periodic:
            return self.world[i % self.h][j % self.w]  # Periodic conditions
        else:
            if (0 <= i < self.h) and (0 <= j < self.w):  # Inside grid
                return self.world[i][j]
            else:                       # Outside grid
                return False            # There's nobody out there...

    def __str__(self):
        """
        Convert the grid to a visually handy string.
        """

        return '\n'.join([''.join([self.cells[val] for val in row])
                          for row in self.world])

    def evolve_cell(self, i, j):
        """
        Tells if cell (*i*,*j*) will survive during game iteration,
        depending on the number of living neighboring cells.
        """

        alive = self.get(i, j)           # Current cell status
        # Count living cells *around* current one (excluding current one)
        count = sum(self.get(i + ii, j + jj)
                    for ii in [-1, 0, 1]
                    for jj in [-1, 0, 1]
                    if (ii, jj) != (0, 0))

        if count == 3:
            # A cell w/ 3 neighbors will either stay alive or resuscitate
            future = True
        elif count < 2 or count > 3:
            # A cell w/ too few or too many neighbors will die
            future = False
        else:
            # A cell w/ 2 or 3 neighbors will stay as it is (dead or alive)
            future = alive             # Current status

        return future

    def evolve(self):
        """
        Evolve the game grid by one step.
        """

        # Update the grid
        self.world = [[self.evolve_cell(i, j)
                       for j in range(self.w)]
                      for i in range(self.h)]

if __name__ == "__main__":

    import time

    h, w = (20, 60)                       # (nrows,ncolumns)

    # Instanciation (including initialization)
    life = Life(h, w, periodic=True)

    while True:                         # Infinite loop! (Ctrl-C to break)
        print life                      # Print current world
        print "\n"
        time.sleep(0.1)                 # Pause a bit
        life.evolve()                   # Evolve world

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as N


def mad(a, axis=None):
    """
    Compute *Median Absolute Deviation* of an array along given axis.
    """

    # Median along given axis, but *keeping* the reduced axis so that
    # result can still broadcast against a.
    med = N.median(a, axis=axis, keepdims=True)
    mad = N.median(N.absolute(a - med), axis=axis)  # MAD along given axis

    return mad

if __name__ == '__main__':

    x = N.arange(4 * 5, dtype=float).reshape(4, 5)

    print "x =\n", x
    print "MAD(x, axis=None)   =", mad(x)
    print "MAD(x, axis=0)      =", mad(x, axis=0)
    print "MAD(x, axis=1)      =", mad(x, axis=1)
    print "MAD(x, axis=(0, 1)) =", mad(x, axis=(0, 1))

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as N


def pull(x, dx):
    """
    Compute the pull from x, dx.

    * Input data: x = [x_i], error dx = [s_i] Optimal
    * (variance-weighted) mean: E = sum(x_i/s_i²)/sum(1/s_i²) Variance
    * on weighted mean: var(E) = 1/sum(1/s_i²) Pull: p_i = (x_i -
    * E_i)/sqrt(var(E_i) + s_i²) where E_i and var(E_i) are computed
    * without point i.

    If errors s_i are correct, the pull distribution is centered on 0
    with standard deviation of 1.
    """

    assert x.ndim == dx.ndim == 1, "pull works on 1D-arrays only."
    assert len(x) == len(dx), "dx should be the same length as x."

    n = len(x)

    i = N.resize(N.arange(n), n * n)  # 0,...,n-1,0,...n-1,..., n times (n²,)
    i[::n + 1] = -1                  # Mark successively 0,1,2,...,n-1
    # Remove marked indices & reshape (n,n-1)
    j = i[i >= 0].reshape((n, n - 1))

    v = dx ** 2                      # Variance
    w = 1 / v                        # Variance (optimal) weighting

    Ei = N.average(x[j], weights=w[j], axis=1)  # Weighted mean (n,)
    vEi = 1 / N.sum(w[j], axis=1)                # Variance of mean (n,)

    p = (x - Ei) / N.sqrt(vEi + v)               # Pull (n,)

    return p

if __name__ == '__main__':

    import matplotlib.pyplot as P
    import scipy.stats as SS

    n = 1000
    mu = 1.
    sig = 2.

    # Normally distributed random sample of size n, with mean=mu and std=sig
    x = N.random.normal(loc=mu, scale=sig, size=n)
    dx = N.full_like(x, sig)              # Formal (true) errors

    p = pull(x, dx)                       # Pull computation

    m, s = p.mean(), p.std(ddof=1)
    print "Pull ({} entries): mean={:.2f}, std={:.2f}".format(n, m, s)

    fig, ax = P.subplots()
    _, bins, _ = ax.hist(p, bins='auto', normed=True,
                         histtype='stepfilled',
                         label=u"#={}, µ={:.3f}, σ={:.3f}".format(n, m, s))
    y = N.linspace(-3, 3)
    ax.plot(y, SS.norm.pdf(y), label=ur"$\mathcal{N}$(µ=0, σ²=1)")
    ax.set(title='Pull distribution', xlabel='Pull')
    ax.legend(loc='upper left')

    P.show()

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as N
import scipy.stats as SS
import matplotlib.pyplot as P


def printStats(x, y):
    """
    Print out means and variances for x and y, as well as correlation
    coeff. (Pearson) and linear regression for y vs. x.
    """

    assert N.shape(x) == N.shape(y), "Incompatible input arrays"

    print "x: mean={:.2f}, variance={:.2f}".format(N.mean(x), N.var(x))
    print "y: mean={:.2f}, variance={:.2f}".format(N.mean(y), N.var(y))
    print "y vs. x: corrcoeff={:.2f}".format(SS.pearsonr(x, y)[0])
    # slope, intercept, r_value, p_value, std_err
    a, b, _, _, _ = SS.linregress(x, y)
    print "y vs. x: y = {:.2f} x + {:.2f}".format(a, b)


def plotStats(ax, x, y, title='', fancy=True):
    """
    Plot y vs. x, and linear regression.
    """

    assert N.shape(x) == N.shape(y), "Incompatible input arrays"

    # slope, intercept, r_value, p_value, std_err
    a, b, r, _, _ = SS.linregress(x, y)

    # Data + corrcoeff
    ax.plot(x, y, 'bo', label="r = {:.2f}".format(r))

    # Linear regression
    xx = N.array([0, 20])
    yy = a * xx + b
    ax.plot(xx, yy, 'r-', label="y = {:.2f} x + {:.2f}".format(a, b))

    leg = ax.legend(loc='upper left', fontsize='small')

    if fancy:                   # Additional stuff
        # Add mean line ± stddev
        m = N.mean(x)
        s = N.std(x, ddof=1)
        ax.axvline(m, color='g', ls='--', label='_')  # Mean
        ax.axvspan(m - s, m + s, color='g', alpha=0.2, label='_')  # Std-dev

        m = N.mean(y)
        s = N.std(y, ddof=1)
        ax.axhline(m, color='g', ls='--', label='_')  # Mean
        ax.axhspan(m - s, m + s, color='g', alpha=0.2, label='_')  # Std-dev

        # Title and labels
        ax.set_title(title)
        if ax.is_last_row():
            ax.set_xlabel("x")
        if ax.is_first_col():
            ax.set_ylabel("y")


if __name__ == '__main__':

    quartet = N.genfromtxt("anscombe.dat")  # Read Anscombe's Quartet

    fig = P.figure()

    for i in range(4):                  # Loop over quartet sets x,y
        ax = fig.add_subplot(2, 2, i + 1)
        print "Dataset #{} ".format(i + 1) + "=" * 20
        x, y = quartet[:, 2 * i:2 * i + 2].T
        printStats(x, y)                # Print main statistics
        plotStats(ax, x, y, title=str(i + 1))  # Plots

    fig.suptitle("Anscombe's Quartet", fontsize='x-large')

    P.show()

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2012-09-05 02:37 ycopin@lyopc469>

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


def iteration(r, niter=100):

    x = random.uniform(0, 1)
    i = 0
    while i < niter and x < 1:
        x = r * x * (1 - x)
        i += 1

    return x if x < 1 else -1


def generate_diagram(r, ntrials=50):
    """
    Cette fonction retourne (jusqu'à) *ntrials* valeurs d'équilibre
    pour les *r* d'entrée.  Elle renvoie un tuple :

    + le premier élément est la liste des valeurs prises par le paramètre *r*
    + le second est la liste des points d'équilibre correspondants
    """

    r_v = []
    x_v = []
    for rr in r:
        j = 0
        while j < ntrials:
            xx = iteration(rr)
            if xx > 0:  # Convergence : il s'agit d'une valeur d'équilibre
                r_v.append(rr)
                x_v.append(xx)
            j += 1                      # Nouvel essai

    return r_v, x_v

r = np.linspace(0, 4, 1000)
x, y = generate_diagram(r)

plt.plot(x, y, 'r,')
plt.xlabel('r')
plt.ylabel('x')
plt.show()

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Time-stamp: <2012-09-05 02:10 ycopin@lyopc469>

"""Visualisation de l'ensemble de julia
<http://fr.wikipedia.org/wiki/Ensemble_de_Julia>`_.

Exercice: proposer des solutions pour accélérer le calcul.
"""

import numpy as np
import matplotlib.pyplot as plt

c = complex(0.284, 0.0122)           # Constante

xlim = 1.5                          # [-xlim,xlim] × i[-xlim,xlim]
nx = 1000                           # Nb de pixels
niter = 100                         # Nb d'itérations

x = np.linspace(-xlim, xlim, nx)    # nx valeurs de -xlim à +xlim
xx, yy = np.meshgrid(x, x)           # Tableaux 2D
z = xx + 1j * yy                      # Portion du plan complexe
for i in range(niter):              # Itération: z(n+1) = z(n)**2 + c
    z = z ** 2 + c

# Visualisation
plt.imshow(np.abs(z), extent=[-xlim, xlim, -xlim, xlim], aspect='equal')
plt.title(c)
plt.show()

%matplotlib inline

import numpy as N
import scipy.integrate as SI
import matplotlib.pyplot as P

g = 9.81 # Pesanteur [m/s2]
cx = 0.45 # Coefficient de frottement d'une sphère
rhoAir = 1.2 # Masse volumique de l'air [kg/m3] au niveau de la mer, T=20°C
rad = 0.1748/2 # Rayon du boulet [m]
rho = 6.23e3 # Masse volumique du boulet [kg/m3]
mass = 4./3.*N.pi*rad**3 * rho # Masse du boulet [kg]
alpha = 0.5*cx*rhoAir*N.pi*rad**2 / mass # Coefficient de frottement par unité de masse
print "Masse du boulet: {:.2f} kg".format(mass)
print "Coefficient de frottement par unité de masse: {} S.I.".format(alpha)

v0 = 450. # Vitesse initiale [m/s]
alt = 45. # Inclinaison du canon [deg]
alt *= N.pi/180. # Inclinaison [rad]
z0 = (0.,0.,v0*N.cos(alt),v0*N.sin(alt)) # (x0,y0,vx0,vy0)

tc = N.sqrt(mass/(g * alpha))
print "Temps caractéristique: {:.1f} s".format(tc)
t = N.linspace(0, tc, 100)

def zdot(z, t):
    """Calcul de la dérivée de z=(x,y,vx,vy) à l'instant t."""

    x,y,vx,vy = z
    alphav = alpha * N.hypot(vx, vy)

    return (vx,vy,-alphav*vx,-g-alphav*vy) # dz/dt = (vx,vy,x..,y..)

ypos = zs[:,1]>=0 # y>0?
print "temps de coll. t(y~0) = {:.0f} s".format(t[ypos][-1]) # Dernier instant pour lequel y>0
print "portée x(y~0) = {:.0f} m".format(zs[ypos,0][-1]) # Portée approximative du canon
#print "y(y~0) = {:.0f} m".format(zs[ypos,1][-1]) # ~0
print "vitesse(y~0): {:.0f} m/s".format(N.hypot(zs[ypos,2][-1],zs[ypos,3][-1]))

fig,ax = P.subplots()
ax.plot(zs[ypos,0],zs[ypos,1], 'b.-')
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_title("Trajectoire du boulet")