Notes de cours scientifiques Python

SciPy est composée de sous-modules spéciaux :

Ils reposent tous sur NumPyhttp://docs.scipy.org/doc/numpy/reference/index.html#numpy, mais sont complètement indépendants les uns des autres. La manière standard pour importer la librairie NumPy et la librairie SciPy est la suivante :

>>> import numpy as np
>>> from scipy import stats  # same for other sub-modules

L'import principal SciPy contient le plus souvent des fonctions identiques aux fonctions NumPy (comparer scipy.cos et np.cos). Celles-ci sont fournies uniquement pour des raisons historiques, il n'y a donc aucune raison d'effectuer un import scipy directement dans votre code.

• Chargement et enregistrement de fichiers Matlab :

>>> from scipy import io as spio
>>> a = np.ones((3, 3))
>>> spio.savemat('file.mat', {'a': a}) # savemat expects a dictionary
>>> data = spio.loadmat('file.mat', struct_as_record=True)
>>> data['a']
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])

• Lecture d'images :

>>> from scipy import misc
>>> misc.imread('fname.png')
>>> # Matplotlib also has a similar function
>>> import matplotlib.pyplot as plt
>>> plt.imread('fname.png')

Voir aussi :

• charger/enregistrer des fichiers texte : numpy.loadtxt()http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html#numpy.loadtxt/numpy.savetxt()http://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html#numpy.savetxt ;
• chargement « intelligent » de fichiers TEXT ou CSV : numpy.genfromtxt()http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html#numpy.genfromtxt/numpy.recfromcsv() ;
• rapide et efficace, mais spécifique à NumPy, pour le format binaire : numpy.save()http://docs.scipy.org/doc/numpy/reference/generated/numpy.save.html#numpy.save/numpy.load()http://docs.scipy.org/doc/numpy/reference/generated/numpy.load.html#numpy.load.

Comme leurs noms l'indiquent, les fonctions spéciales ne sont pas des fonctions ordinaires. Les docstrings (ou « commentaires dans le code ») du module scipy.specialhttp://docs.scipy.org/doc/scipy/reference/special.html#scipy.special sont suffisamment bien détaillées pour ne pas en faire la liste ici.

Les fonctions les plus souvent utilisées sont :

• fonction de Bessel, comme scipy.special.jn() (fonction de Bessel d'ordre n) ;
• fonction elliptique (scipy.special.ellipj() pour la fonction Jacobienne elliptique) ;
• fonction Gamma : scipy.special.gamma(), aussi notée scipy.special.gammaln() qui donnera le logarithme du Gamma pour une précision numérique plus grande ;
• Erf, la surface sous une courbe de distribution gaussienne : scipy.special.erf().

Le module scipy.linalghttp://docs.scipy.org/doc/scipy/reference/linalg.html#scipy.linalg propose tout un tas d'opérations d'algèbre linéaire standards, reposant sur une implémentation qui a déjà fait ses preuves (BLAS, LAPACK).

• La fonction scipy.linalg.det()http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.det.html#scipy.linalg.det calcule le déterminant de la matrice carrée :

>>> from scipy import linalg
>>> arr = np.array([[1, 2],
...                 [3, 4]])
>>> linalg.det(arr)
-2.0
>>> arr = np.array([[3, 2],
...                 [6, 4]])
>>> linalg.det(arr)
0.0
>>> linalg.det(np.ones((3, 4)))
Traceback (most recent call last):
...
ValueError: expected square matrix

• La fonction scipy.linalg.inv()http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.inv.html#scipy.linalg.inv calcule l'inverse de la matrice carrée :

>>> arr = np.array([[1, 2],
...                 [3, 4]])
>>> iarr = linalg.inv(arr)
>>> iarr
array([[-2. ,  1. ],
       [ 1.5, -0.5]])
>>> np.allclose(np.dot(arr, iarr), np.eye(2))
True

Enfin, calculer l'inverse d'une matrice singulière (ou non inversible, car son déterminant est zéro) remontera l'erreur LinAlgError :

>>> arr = np.array([[3, 2],
...                 [6, 4]])
>>> linalg.inv(arr)
Traceback (most recent call last):
...
LinAlgError: singular matrix

• Des opérations plus complexes sont disponibles, par exemple la méthode de décomposition en valeurs singulières (SVD) :

>>> arr = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1])
>>> uarr, spec, vharr = linalg.svd(arr)

Le résultat de la décomposition est :

>>> spec    
array([ 14.88982544,   0.45294236,   0.29654967])

La matrice originale peut être recomposée en réalisant une multiplication matricielle des sorties de svd avec np.dot :

>>> sarr = np.diag(spec)
>>> svd_mat = uarr.dot(sarr).dot(vharr)
>>> np.allclose(svd_mat, arr)
True

SVD est très souvent utilisée en statistiques et en traitement du signal. Beaucoup d'autres méthodes de décomposition standard (QR, LU, Cholesky, Schur), ainsi que des méthodes de résolution de systèmes linéaires sont également disponibles dans scipy.linalghttp://docs.scipy.org/doc/scipy/reference/linalg.html#scipy.linalg.

Le module http://docs.scipy.org/doc/scipy/reference/fftpack.html#scipy.fftpackhttp://docs.scipy.org/doc/scipy/reference/fftpack.html#scipy.fftpack permet de calculer des transformées de Fourier. Comme dans l'illustration, un signal d'entrée bruité pourrait ressembler à :

>>> time_step = 0.02
>>> period = 5.
>>> time_vec = np.arange(0, 20, time_step)
>>> sig = np.sin(2 * np.pi / period * time_vec) + \
...       0.5 * np.random.randn(time_vec.size)

L'observateur ne connaît pas la fréquence du signal, mais uniquement la fréquence d'échantillonnage (échantillonnage dans le domaine temporel) du signal sig. Le signal est supposé provenir d'une fonction réelle de telle sorte que la transformée de Fourier soit symétrique. La fonction scipy.fftpack.fftfreq()http://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.fftfreq.html#scipy.fftpack.fftfreq va générer les fréquences d'échantillonnage et scipy.fftpack.fft()http://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.fft.html#scipy.fftpack.fft va calculer la transformée de Fourier.

>>> from scipy import fftpack
>>> sample_freq = fftpack.fftfreq(sig.size, d=time_step)
>>> sig_fft = fftpack.fft(sig)

Parce que la puissance (power) résultante est symétrique, seule la partie positive du spectre suffira à trouver la fréquence :

>>> pidxs = np.where(sample_freq > 0)
>>> freqs = sample_freq[pidxs]
>>> power = np.abs(sig_fft)[pidxs]

import numpy as np
from scipy import fftpack
import pylab as pl


np.random.seed(1234)

time_step = 0.02
period = 5.

time_vec = np.arange(0, 20, time_step)
sig = np.sin(2 * np.pi / period * time_vec) + \
      0.5 * np.random.randn(time_vec.size)

sample_freq = fftpack.fftfreq(sig.size, d=time_step)
sig_fft = fftpack.fft(sig)
pidxs = np.where(sample_freq > 0)
freqs, power = sample_freq[pidxs], np.abs(sig_fft)[pidxs]
freq = freqs[power.argmax()]

pl.figure()
pl.plot(freqs, power)
pl.xlabel('Frequency [Hz]')
pl.ylabel('plower')
axes = pl.axes([0.3, 0.3, 0.5, 0.5])
pl.title('Peak frequency')
pl.plot(freqs[:8], power[:8])
pl.setp(axes, yticks=[])

La fréquence du signal peut être obtenue ainsi :

>>> freq = freqs[power.argmax()]
>>> np.allclose(freq, 1./period)  # check that correct freq is found
True

À présent, le bruit haute fréquence va être supprimé à partir du signal de la transformée de Fourier :

>>> sig_fft[np.abs(sample_freq) > freq] = 0

Le signal filtré résultant peut être obtenu par la fonction de Fourier inverse scipy.fftpack.ifft()http://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.ifft.html#scipy.fftpack.ifft :

>>> main_sig = fftpack.ifft(sig_fft)

Le résultat sera visible avec :

>>> import pylab as plt
>>> plt.figure()
>>> plt.plot(time_vec, sig)
>>> plt.plot(time_vec, main_sig, linewidth=3)
>>> plt.xlabel('Time [s]')
>>> plt.ylabel('Amplitude')

import numpy as np
from scipy import fftpack
import pylab as pl

np.random.seed(1234)

time_step = 0.02
period = 5.

time_vec = np.arange(0, 20, time_step)
sig = np.sin(2 * np.pi / period * time_vec) + \
      0.5 * np.random.randn(time_vec.size)

sample_freq = fftpack.fftfreq(sig.size, d=time_step)
sig_fft = fftpack.fft(sig)
pidxs = np.where(sample_freq > 0)
freqs, power = sample_freq[pidxs], np.abs(sig_fft)[pidxs]
freq = freqs[power.argmax()]

sig_fft[np.abs(sample_freq) > freq] = 0
main_sig = fftpack.ifft(sig_fft)

pl.figure()
pl.plot(time_vec, sig)
pl.plot(time_vec, main_sig, linewidth=3)
pl.xlabel('Time [s]')
pl.ylabel('Amplitude')

"""
Discover the periods in ../../data/populations.txt
"""
import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('../../data/populations.txt')
years = data[:, 0]
populations = data[:, 1:]

ft_populations = np.fft.fft(populations, axis=0)
frequencies = np.fft.fftfreq(populations.shape[0], years[1] - years[0])
periods = 1 / frequencies

plt.figure()
plt.plot(years, populations * 1e-3)
plt.xlabel('Year')
plt.ylabel('Population number ($\cdot10^3$)')
plt.legend(['hare', 'lynx', 'carrot'], loc=1)

plt.figure()
plt.plot(periods, abs(ft_populations) * 1e-3, 'o')
plt.xlim(0, 22)
plt.xlabel('Period')
plt.ylabel('Power ($\cdot10^3$)')

plt.show()

# There's probably a period of around 10 years (obvious from the
# plot), but for this crude a method, there's not enough data to say
# much more.

"""
Simple image blur by convolution with a Gaussian kernel
"""

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt

# read image
img = plt.imread('../../data/elephant.png')

# prepare an 1-D Gaussian convolution kernel
t = np.linspace(-10, 10, 30)
bump = np.exp(-0.1*t**2)
bump /= np.trapz(bump) # normalize the integral to 1

# make a 2-D kernel out of it
kernel = bump[:, np.newaxis] * bump[np.newaxis, :]

# padded fourier transform, with the same shape as the image
kernel_ft = fftpack.fft2(kernel, shape=img.shape[:2], axes=(0, 1))

# convolve
img_ft = fftpack.fft2(img, axes=(0, 1))
img2_ft = kernel_ft[:, :, np.newaxis] * img_ft
img2 = fftpack.ifft2(img2_ft, axes=(0, 1)).real

# clip values to range
img2 = np.clip(img2, 0, 1)

# plot output
plt.imshow(img2)
plt.show()

# Further exercise (only if you are familiar with this stuff):
#
# A "wrapped border" appears in the upper left and top edges of the
# image. This is because the padding is not done correctly, and does
# not take the kernel size into account (so the convolution "flows out
# of bounds of the image").  Try to remove this artifact.

L'optimisation est le procédé numérique de recherche d'un extremum ou d'un point d'inflexion d'une fonction (c.-à-d. les points où le gradient de la fonction est nul).

Le module scipy.optimizehttp://docs.scipy.org/doc/scipy/reference/optimize.html#scipy.optimize propose de nombreux algorithmes très utiles pour la recherche de minimum de fonctions (scalaires ou multidimensionnelles), l'approximation de courbes et la recherche des racines.

>>> from scipy import optimize

Trouver le minimum d'une fonction scalaire

On définit la fonction suivante :

>>> def f(x):
...     return x**2 + 10*np.sin(x)

et on l'affiche :

>>> x = np.arange(-10, 10, 0.1)
>>> plt.plot(x, f(x)) 
>>> plt.show()

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return x**2 + 10*np.sin(x)


x = np.arange(-10, 10, 0.1)
plt.plot(x, f(x))

Cette fonction possède un minimum global autour de -1.3 et un minimum local autour de 3.8.

La manière la plus commune et performante pour trouver un minimum de cette fonction est de réaliser une descente de gradient en commençant à partir d'un point initial donné. L'algorithme BFGS est une bonne façon pour réaliser cela :

>>> optimize.fmin_bfgs(f, 0)
Optimization terminated successfully.
         Current function value: -7.945823
         Iterations: 5
         Function evaluations: 24
         Gradient evaluations: 8
array([-1.30644003])

Toutefois, avec cette approche, on peut se retrouver face à un problème. Si la fonction possède des minimums locaux, l'algorithme peut très bien trouver ces minimums locaux à la place du minimum global en fonction du point initial donné :

>>> optimize.fmin_bfgs(f, 3, disp=0)
array([ 3.83746663])

Si nous ne connaissons pas le voisinage du minimum global afin de choisir un point initial cohérent, nous aurons besoin de recourir à une optimisation globale plus coûteuse en temps de calcul. Pour trouver le minimum global, l'algorithme le plus simple est le « brute force », la fonction est ainsi évaluée en chaque point de la grille donnée :

>>> grid = (-10, 10, 0.1)
>>> xmin_global = optimize.brute(f, (grid,))
>>> xmin_global
array([-1.30641113])

Pour des grilles de tailles plus importantes, scipy.optimize.brute()http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brute.html#scipy.optimize.brute devient assez lent. scipy.optimize.anneal()http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.anneal.html#scipy.optimize.anneal propose une alternative utilisant le « recuit simulé » (simulated annealing). Des algorithmes plus performants existent pour différents types de problèmes d'optimisation globale, mais cela est hors du cadre de cette présentation du module SciPy. Néanmoins, voici quelques librairies utiles consacrées à l'optimisation globale : OpenOpthttp://openopt.org/Welcome, IPOPThttps://github.com/xuy/pyipopt, PyGMOhttp://pagmo.sourceforge.net/pygmo/index.html et PyEvolvehttp://pyevolve.sourceforge.net/.

Afin de trouver le minimum local de la fonction, on peut aussi contraindre le domaine d'étude dans l'intervalle (0, 10) en utilisant scipy.optimize.fminbound()http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fminbound.html#scipy.optimize.fminbound :

>>> xmin_local = optimize.fminbound(f, 0, 10)    
>>> xmin_local
3.8374671...

Rechercher les racines d'une fonction scalaire

Afin de trouver une racine de la fonction f ci-dessous, c.-à-d. un point où f(x)=0, nous pouvons utiliser par exemple la fonction scipy.optimize.fsolve()http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve :

>>> root = optimize.fsolve(f, 1)  # our initial guess is 1
>>> root
array([ 0.])

Notez qu'ici, une seule racine a été trouvée. En inspectant le graphe de f, on remarque qu'il y a une autre racine autour de -2.5. Nous retrouvons la valeur exacte de cette racine en ajustant notre estimation initiale :

>>> root2 = optimize.fsolve(f, -2.5)
>>> root2
array([-2.47948183])

Approximation de courbes

Supposons que nous ayons des données échantillonnées de la fonction f auxquelles nous ajoutons un bruit.

>>> xdata = np.linspace(-10, 10, num=20)
>>> ydata = f(xdata) + np.random.randn(xdata.size)

À présent, si nous connaissons la fonction à partir de laquelle les échantillons ont été tirés (x^2 + sin(x) dans ce cas), mais pas les amplitudes associées, nous pouvons toutefois les trouver en approchant la courbe par la méthode des moindres carrés. Pour commencer, nous devons définir la fonction à approximer :

>>> def f2(x, a, b):
...     return a*x**2 + b*np.sin(x)

Ensuite, nous allons utiliser scipy.optimize.curve_fit()http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit pour trouver a et b :

>>> guess = [2, 2]
>>> params, params_covariance = optimize.curve_fit(f2, xdata, ydata, guess)
>>> params
array([ 0.99925147,  9.76065551])

Maintenant que nous avons trouvé les minimums et les racines de f et que nous avons utilisé l'approximation de courbe sur f, nous allons mettre tous les résultats ensemble sur le même graphe :

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

x = np.arange(-10, 10, 0.1)
def f(x):
    return x**2 + 10*np.sin(x)


grid = (-10, 10, 0.1)
xmin_global = optimize.brute(f, (grid,))
xmin_local = optimize.fminbound(f, 0, 10)
root = optimize.fsolve(f, 1)  # our initial guess is 1
root2 = optimize.fsolve(f, -2.5)

xdata = np.linspace(-10, 10, num=20)
np.random.seed(1234)
ydata = f(xdata) + np.random.randn(xdata.size)

def f2(x, a, b):
    return a*x**2 + b*np.sin(x)

guess = [2, 2]
params, params_covariance = optimize.curve_fit(f2, xdata, ydata, guess)


fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, f(x), 'b-', label="f(x)")
ax.plot(x, f2(x, *params), 'r--', label="Curve fit result")
xmins = np.array([xmin_global[0], xmin_local])
ax.plot(xmins, f(xmins), 'go', label="Minima")
roots = np.array([root, root2])
ax.plot(roots, f(roots), 'kv', label="Roots")
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('f(x)')

Vous trouverez les algorithmes avec les mêmes fonctionnalités pour le cas multidimensionnel dans scipy.optimizehttp://docs.scipy.org/doc/scipy/reference/optimize.html#scipy.optimize.

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def sixhump(x):
    return (4 - 2.1*x[0]**2 + x[0]**4 / 3.) * x[0]**2 + x[0] * x[1] + (-4 + \
        4*x[1]**2) * x[1] **2

x = np.linspace(-2, 2)
y = np.linspace(-1, 1)
xg, yg = np.meshgrid(x, y)

#plt.figure()  # simple visualization for use in tutorial
#plt.imshow(sixhump([xg, yg]))
#plt.colorbar()

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(xg, yg, sixhump([xg, yg]), rstride=1, cstride=1,
                       cmap=plt.cm.jet, linewidth=0, antialiased=False)

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
ax.set_title('Six-hump Camelback function')

Si vous souhaitez voir un exemple plus avancé, consultez l'exercice sur l'approximation de courbes par la méthode des moindres carrés : application pour extraire des points de mesure sur des données topographiques LiDAR.http://scipy-lectures.github.io/intro/summary-exercises/optimize-fit.html#summary-exercise-optimize

Le module scipy.statshttp://docs.scipy.org/doc/scipy/reference/stats.html#scipy.stats contient des outils statistiques et des descriptions probabilistes sur les nombres aléatoires. Les générateurs de nombres aléatoires de nombreux processus aléatoires sont disponibles dans numpy.random.

VII-A. Histogramme et fonction de densité de probabilité▲

Les observations d'un processus aléatoire peuvent être décrites par leur histogramme qui est un estimateur de la densité de probabilité PDF (probability density function) du process aléatoire :

 

Sélectionnez

>>> a = np.random.normal(size=1000)
>>> bins = np.arange(-4, 5)
>>> bins
array([-4, -3, -2, -1,  0,  1,  2,  3,  4])
>>> histogram = np.histogram(a, bins=bins, normed=True)[0]
>>> bins = 0.5*(bins[1:] + bins[:-1])
>>> bins
array([-3.5, -2.5, -1.5, -0.5,  0.5,  1.5,  2.5,  3.5])
>>> from scipy import stats
>>> b = stats.norm.pdf(bins)  # norm is a distribution

 

Sélectionnez

In [1]: pl.plot(bins, histogram)
In [2]: pl.plot(bins, b)

normal_distribution.py

Cacher/Afficher le codeSélectionnez

from scipy import stats
import numpy as np
import pylab as pl
a = np.random.normal(size=10000)
bins = np.linspace(-5, 5, 30)
histogram, bins = np.histogram(a, bins=bins, normed=True)
bins = 0.5*(bins[1:] + bins[:-1])
from scipy import stats
b = stats.norm.pdf(bins)
pl.plot(bins, histogram)
pl.plot(bins, b)

Si on connaît le processus aléatoire tel que le processus normal (gaussien), on peut réaliser une approximation par maximum de vraisemblance des observations pour estimer les paramètres de la distribution sous-jacente du processus aléatoire. Dans cet exemple, on réalise une approximation des données observées par un processus normal (gaussien) :

 

Sélectionnez

>>> loc, std = stats.norm.fit(a)
>>> loc     
-0.045256707490...
>>> std     
0.9870331586690...

Exercice : distributions de probabilité

Générez 1000 nombres aléatoires à partir d'une distribution gamma avec un paramètre de forme égal à 1, puis tracez l'histogramme à partir de cet échantillon de données. Pouvez-vous tracer la PDF par-dessus (cela devrait correspondre, n'est-ce pas ?)

Supplément : les distributions possèdent de nombreuses méthodes utiles. Examinez-les en lisant la docstring (commentaire dans le code) ou en utilisant la complétion IPython. Pouvez-vous retrouver le paramètre de forme = 1 en utilisant la méthode d'approximation fit sur votre jeu d'essai  ?

VII-B. Centiles (Percentiles en anglais souvent utilisé abusivement en français)▲

>>> np.median(a)     
-0.058028034...

>>> stats.scoreatpercentile(a, 50)     
-0.0580280347...

>>> stats.scoreatpercentile(a, 90)     
1.231593551...

VII-C. Tests statistiques▲

>>> a = np.random.normal(0, 1, size=100)
>>> b = np.random.normal(1, 1, size=10)
>>> stats.ttest_ind(a, b)   
(-3.75832707..., 0.00027786...)

Le module scipy.interpolatehttp://docs.scipy.org/doc/scipy/reference/interpolate.html#scipy.interpolate est pratique pour approximer une fonction à partir de données expérimentales et donc d'évaluer les points où il n'existe pas de mesure réelle. Ce module est basé sur la routine FITPACK Fortranhttp://www.netlib.org/dierckx/index.html du projet netlib.

Imaginons des données expérimentales proches de la fonction sinus :

>>> measured_time = np.linspace(0, 1, 10)
>>> noise = (np.random.random(10)*2 - 1) * 1e-1
>>> measures = np.sin(2 * np.pi * measured_time) + noise

La classe scipy.interpolate.interp1dhttp://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d permet de construire une fonction d'interpolation linéaire :

>>> from scipy.interpolate import interp1d
>>> linear_interp = interp1d(measured_time, measures)

L'instance scipy.interpolate.linear_interp a ensuite besoin d'être évaluée dans l'intervalle de temps considéré :

>>> computed_time = np.linspace(0, 1, 50)
>>> linear_results = linear_interp(computed_time)

On peut également réaliser une interpolation cubique à l'aide de mots-clés passés en argument :

>>> cubic_interp = interp1d(measured_time, measures, kind='cubic')
>>> cubic_results = cubic_interp(computed_time)

L'ensemble de ces résultats est visible sur la figure Matplotlib suivante :

"""Generate the interpolation.png image for the interpolate
section of the Scipy tutorial
"""

import numpy as np
from scipy.interpolate import interp1d
import pylab as pl

measured_time = np.linspace(0, 1, 10)
noise = (np.random.random(10)*2 - 1) * 1e-1
measures = np.sin(2 * np.pi * measured_time) + noise

linear_interp = interp1d(measured_time, measures)
computed_time = np.linspace(0, 1, 50)
linear_results = linear_interp(computed_time)
cubic_interp = interp1d(measured_time, measures, kind='cubic')
cubic_results = cubic_interp(computed_time)

pl.plot(measured_time, measures, 'o', ms=6, label='measures')
pl.plot(computed_time, linear_results, label='linear interp')
pl.plot(computed_time, cubic_results, label='cubic interp')
pl.legend()

scipy.interpolate.interp2dhttp://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html#scipy.interpolate.interp2d est l'équivalent au module scipy.interpolate.interp1dhttp://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d pour les données en 2D. Pour les interpolations, le domaine calculé (interpolé) doit rester dans le domaine mesuré. Voir l'exercice sur la Prédiction de la vitesse maximale du venthttp://scipy-lectures.github.io/intro/summary-exercises/stats-interpolate.html#summary-exercise-stat-interp de la station de Sprogø pour avoir un exemple sur la méthode d'interpolation Spline (le mot anglais reste d'usage !)

La routine d'intégration la plus commune est scipy.integrate.quad()http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad :

>>> from scipy.integrate import quad
>>> res, err = quad(np.sin, 0, np.pi/2)
>>> np.allclose(res, 1)
True
>>> np.allclose(err, 1 - res)
True

D'autres schémas d'intégration sont disponibles comme fixed_quad, quadrature, romberg.

Le module scipy.integratehttp://docs.scipy.org/doc/scipy/reference/integrate.html#scipy.integrate propose également des routines pour intégrer des Équations Différentielles Ordinaires (ODE en anglais). En particulier, scipy.integrate.odeint()http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint qui permet de réaliser des intégrations en utilisant la méthode LSODA (Livermore Solver for Ordinary Differential equations avec un aiguillage automatique pour résoudre soit des équations différentielles dites « raides », soit « non raides »).

scipy.integrate.odeint() résout des systèmes d'ODE de la forme :

dy/dt = rhs(y1, y2, .., t0,...)

Par exemple, essayons de résoudre l'ODE dy/dt = -2y dans l'intervalle t = 0..4 avec une condition initiale y(t=0) = 1. En premier lieu, la fonction dérivée dy/dt (variation de y en fonction du temps) doit être définie :

>>> def calc_derivative(ypos, time, counter_arr):
...     counter_arr += 1
...     return -2 * ypos
...

Une variable supplémentaire counter_arr a été ajoutée pour montrer que la fonction peut être appelée plusieurs fois de suite dans le même pas de temps, et cela, jusqu'à la convergence de la solution.

Le compteur (c'est un tableau) est défini ainsi :

>>> counter = np.zeros((1,), dtype=np.uint16)

La trajectoire sera calculée de la manière suivante :

>>> from scipy.integrate import odeint
>>> time_vec = np.linspace(0, 4, 40)
>>> yvec, info = odeint(calc_derivative, 1, time_vec,
...                     args=(counter,), full_output=True)

Par conséquent, la fonction dérivée a été appelée plus de 40 fois (ce qui était initialement le nombre de pas de temps) :

>>> counter
array([129], dtype=uint16)

Le nombre cumulé d'itérations pour chacun des 10 premiers pas de temps peut être obtenu par :

>>> info['nfe'][:10]
array([31, 35, 43, 49, 53, 57, 59, 63, 65, 69], dtype=int32)

À noter que la méthode de résolution a besoin de plus d'itérations lors du premier pas de temps. La solution yvec pour la trajectoire est tracée ci-dessous :

"""Solve the ODE dy/dt = -2y between t = 0..4, with the
initial condition y(t=0) = 1.
"""

import numpy as np
from scipy.integrate import odeint
import pylab as pl

def calc_derivative(ypos, time):
    return -2*ypos

time_vec = np.linspace(0, 4, 40)
yvec = odeint(calc_derivative, 1, time_vec)

pl.plot(time_vec, yvec)
pl.xlabel('Time [s]')
pl.ylabel('y position [m]')

Un autre exemple d'utilisation de scipy.integrate.odeint()http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint pourrait être celui d'un oscillateur harmonique, c.-à-d. oscillateur masse-ressort (oscillateur du 2^d ordre). La position d'une masse suspendue à un ressort obéit à une ODE du 2^d degré du type : y'' + 2 eps wo y' + wo^2 y = 0 avec wo^2 = k/m, k représentant la constante de raideur du ressort, m la masse et eps=c/(2 m wo) avec c le coefficient d'amortissement du fluide environnant. Pour cet exemple, nous avons choisi les paramètres suivants :

>>> mass = 0.5  # kg
>>> kspring = 4  # N/m
>>> cviscous = 0.4  # N s/m

Avec ces constantes, le système oscillant sera amorti peu à peu jusqu'à l'équilibre à cause de :

>>> eps = cviscous / (2 * mass * np.sqrt(kspring/mass))
>>> eps < 1
True

Avec la méthode de résolution scipy.integrate.odeint()http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint, l'équation du 2^d degré doit être transformée en un système de deux équations du 1^er degré avec le vecteur Y=(y, y').

Il sera alors pratique de définir nu = 2 eps * wo = c / m et om = wo^2 = k/m :

>>> nu_coef = cviscous / mass
>>> om_coef = kspring / mass

La fonction suivante va calculer la vitesse et l'accélération ainsi :

>>> def calc_deri(yvec, time, nuc, omc):
...     return (yvec[1], -nuc * yvec[1] - omc * yvec[0])
...
>>> time_vec = np.linspace(0, 10, 100)
>>> yarr = odeint(calc_deri, (1, 0), time_vec, args=(nu_coef, om_coef))

La position finale du ressort et la vitesse d'oscillation sont visibles sur la figure Matplotlib suivante :

"""Damped spring-mass oscillator
"""

import numpy as np
from scipy.integrate import odeint
import pylab as pl

mass = 0.5
kspring = 4
cviscous = 0.4

nu_coef = cviscous / mass
om_coef = kspring / mass

def calc_deri(yvec, time, nuc, omc):
    return (yvec[1], -nuc * yvec[1] - omc * yvec[0])

time_vec = np.linspace(0, 10, 100)
yarr = odeint(calc_deri, (1, 0), time_vec, args=(nu_coef, om_coef))

pl.plot(time_vec, yarr[:, 0], label='y')
pl.plot(time_vec, yarr[:, 1], label="y'")
pl.legend()

N. B. Il n'y a pas de méthode de résolution d'équations différentielles partielles (PDE) dans le module SciPy. Plusieurs librairies Python sont disponibles pour résoudre les PDE, tels que fipyhttp://www.ctcms.nist.gov/fipy/ ou SfePyhttp://code.google.com/p/sfepy/.

>>> from scipy import signal

scipy.signal.detrend()http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.detrend.html#scipy.signal.detrend supprime les composantes linéaires du signal :

t = np.linspace(0, 5, 100)
x = t + np.random.normal(size=100)

pl.plot(t, x, linewidth=3)
pl.plot(t, signal.detrend(x), linewidth=3)

import numpy as np
import pylab as pl
from scipy import signal
t = np.linspace(0, 5, 100)
x = t + np.random.normal(size=100)

pl.plot(t, x, linewidth=3)
pl.plot(t, signal.detrend(x), linewidth=3)

scipy.signal.resample()http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html#scipy.signal.resample, ré-échantillonnage du signal en n points à l'aide de la FFT.

t = np.linspace(0, 5, 100)
x = np.sin(t)

pl.plot(t, x, linewidth=3)
pl.plot(t[::2], signal.resample(x, 50), 'ko')

import numpy as np
import pylab as pl
from scipy import signal
t = np.linspace(0, 5, 100)
x = np.sin(t)

pl.plot(t, x, linewidth=3)
pl.plot(t[::2], signal.resample(x, 50), 'ko')

scipy.signalhttp://docs.scipy.org/doc/scipy/reference/signal.html#scipy.signal possède beaucoup de fonctions de fenêtrage comme scipy.signal.hamming()http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hamming.html#scipy.signal.hamming, scipy.signal.bartlett()http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.bartlett.html#scipy.signal.bartlett, scipy.signal.blackman()http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.blackman.html#scipy.signal.blackman…

scipy.signalhttp://docs.scipy.org/doc/scipy/reference/signal.html#scipy.signal possède également des fonctions filtres (filtre médian scipy.signal.medfilt()http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.medfilt.html#scipy.signal.medfilt, filtre Wiener scipy.signal.wiener()http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.wiener.html#scipy.signal.wiener). Ces fonctions filtres seront décrites dans le chapitre consacré au traitement d'images plus loin.

Le sous-module consacré au traitement d'images dans SciPy est scipy.ndimagehttp://docs.scipy.org/doc/scipy/reference/ndimage.html#scipy.ndimage.

>>> from scipy import ndimage

Les routines de traitement d'images sont classées par catégorie de traitement qu'elles réalisent.

XI-A. Transformations géométriques sur les images▲

Pour modifier l'orientation, la résolution…

 

Sélectionnez

>>> from scipy import misc
>>> lena = misc.lena()
>>> shifted_lena = ndimage.shift(lena, (50, 50))
>>> shifted_lena2 = ndimage.shift(lena, (50, 50), mode='nearest')
>>> rotated_lena = ndimage.rotate(lena, 30)
>>> cropped_lena = lena[50:-50, 50:-50]
>>> zoomed_lena = ndimage.zoom(lena, 2)
>>> zoomed_lena.shape
(1024, 1024)

 

Sélectionnez

In [35]: subplot(151)
Out[35]: <matplotlib.axes.AxesSubplot object at 0x925f46c>

In [36]: pl.imshow(shifted_lena, cmap=cm.gray)
Out[36]: <matplotlib.image.AxesImage object at 0x9593f6c>

In [37]: axis('off')
Out[37]: (-0.5, 511.5, 511.5, -0.5)

In [39]: # etc.

XI-B. Filtrage d'images▲

 

Sélectionnez

>>> from scipy import misc
>>> lena = misc.lena()
>>> import numpy as np
>>> noisy_lena = np.copy(lena).astype(np.float)
>>> noisy_lena += lena.std()*0.5*np.random.standard_normal(lena.shape)
>>> blurred_lena = ndimage.gaussian_filter(noisy_lena, sigma=3)
>>> median_lena = ndimage.median_filter(blurred_lena, size=5)
>>> from scipy import signal
>>> wiener_lena = signal.wiener(blurred_lena, (5,5))

De nombreux autres filtres peuvent être appliqués aux images et sont disponibles dans scipy.ndimage.filtershttp://docs.scipy.org/doc/scipy/reference/ndimage.html#scipy.ndimage.filters et scipy.signalhttp://docs.scipy.org/doc/scipy/reference/signal.html#scipy.signal.

Exercice

Comparez les histogrammes sur une image filtrée par différents filtres.

XI-C. Morphologie mathématique▲

La morphologie en mathématique est une théorie mathématique qui découle d'une série de théories. Elle permet de caractériser et de transformer des structures géométriques. Les images binaires (noir et blanc) en particulier, peuvent être transformées en utilisant cette théorie : la série de pixels à transformer correspond à celle des pixels de valeur non nulle au voisinage de cette série. La théorie s'étend également aux images en niveau de gris.

Les opérations élémentaires de morphologie en mathématique utilisent une « structure » (structuring element) afin de modifier d'autres structures géométriques.

Pour commencer, générons une structure élémentaire :

 

Sélectionnez

>>> el = ndimage.generate_binary_structure(2, 1)
>>> el
array([[False, True, False],
       [True, True, True],
       [False, True, False]], dtype=bool)
>>> el.astype(np.int)
array([[0, 1, 0],
       [1, 1, 1],
       [0, 1, 0]])

Érosion

 

Sélectionnez

>>> a = np.zeros((7,7), dtype=np.int)
>>> a[1:6, 2:5] = 1
>>> a
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 0]])
>>> ndimage.binary_erosion(a).astype(a.dtype)
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0]])
>>> #Erosion removes objects smaller than the structure
>>> ndimage.binary_erosion(a, structure=np.ones((5,5))).astype(a.dtype)
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0]])

Dilatation

 

Sélectionnez

>>> a = np.zeros((5, 5))
>>> a[2, 2] = 1
>>> a
array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])
>>> ndimage.binary_dilation(a).astype(a.dtype)
array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.],
       [ 0.,  1.,  1.,  1.,  0.],
       [ 0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])

Ouverture

 

Sélectionnez

>>> a = np.zeros((5,5), dtype=np.int)
>>> a[1:4, 1:4] = 1; a[4, 4] = 1
>>> a
array([[0, 0, 0, 0, 0],
       [0, 1, 1, 1, 0],
       [0, 1, 1, 1, 0],
       [0, 1, 1, 1, 0],
       [0, 0, 0, 0, 1]])
>>> # Opening removes small objects
>>> ndimage.binary_opening(a, structure=np.ones((3,3))).astype(np.int)
array([[0, 0, 0, 0, 0],
       [0, 1, 1, 1, 0],
       [0, 1, 1, 1, 0],
       [0, 1, 1, 1, 0],
       [0, 0, 0, 0, 0]])
>>> # Opening can also smooth corners
>>> ndimage.binary_opening(a).astype(np.int)
array([[0, 0, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 1, 1, 1, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 0, 0]])

Fermeture : ndimage.binary_closing

Exercice

Vérifier que l'ouverture correspond à l'érosion puis à la dilatation.

Une opération d'ouverture supprime les petites structures de l'image, alors que l'opération de fermeture va remplir les petits trous. De telles opérations peuvent donc être utilisées pour « nettoyer » une image.

 

Sélectionnez

>>> a = np.zeros((50, 50))
>>> a[10:-10, 10:-10] = 1
>>> a += 0.25*np.random.standard_normal(a.shape)
>>> mask = a>=0.5
>>> opened_mask = ndimage.binary_opening(mask)
>>> closed_mask = ndimage.binary_closing(opened_mask)

Exercice

Vérifier que la surface du carré reconstruit est plus petite que la surface du carré original.

(Le contraire serait arrivé si nous avions utilisé l'opération de fermeture avant celle d'ouverture).

Pour les images en niveau de gris, l'érosion (respectivement la dilatation) se traduit par le remplacement d'un pixel par la valeur minimale (respectivement maximale) parmi les pixels couverts par la structure élémentaire centrée en ce pixel.

 

Sélectionnez

>>> a = np.zeros((7,7), dtype=np.int)
>>> a[1:6, 1:6] = 3
>>> a[4,4] = 2; a[2,3] = 1
>>> a
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 3, 3, 3, 3, 3, 0],
       [0, 3, 3, 1, 3, 3, 0],
       [0, 3, 3, 3, 3, 3, 0],
       [0, 3, 3, 3, 2, 3, 0],
       [0, 3, 3, 3, 3, 3, 0],
       [0, 0, 0, 0, 0, 0, 0]])
>>> ndimage.grey_erosion(a, size=(3,3))
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 3, 2, 2, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0]])

XI-D. Mesures sur les images▲

Commençons par générer une jolie image de synthèse binaire.

 

Sélectionnez

>>> x, y = np.indices((100, 100))
>>> sig = np.sin(2*np.pi*x/50.)*np.sin(2*np.pi*y/50.)*(1+x*y/50.**2)**2
>>> mask = sig > 1

À présent, recherchons quelques informations sur l'image :

 

Sélectionnez

>>> labels, nb = ndimage.label(mask)
>>> nb
8
>>> areas = ndimage.sum(mask, labels, xrange(1, labels.max()+1))
>>> areas
array([ 190.,   45.,  424.,  278.,  459.,  190.,  549.,  424.])
>>> maxima = ndimage.maximum(sig, labels, xrange(1, labels.max()+1))
>>> maxima
array([  1.80238238,   1.13527605,   5.51954079,   2.49611818,
         6.71673619,   1.80238238,  16.76547217,   5.51954079])
>>> ndimage.find_objects(labels==4)
[(slice(30L, 48L, None), slice(30L, 48L, None))]
>>> sl = ndimage.find_objects(labels==4)
>>> import pylab as pl
>>> pl.imshow(sig[sl[0]])   
<matplotlib.image.AxesImage object at ...>

Pour en savoir plus, voir l'exercice d'application en traitement d'images : comptage de bulles et de grains.

Les exercices suivants utilisent principalement NumPy, SciPy et Matplotlib. Ils proposent quelques exemples concrets de calculs scientifiques exécutés avec Python. Maintenant que vous vous êtes familiarisé avec NumPy et SciPy, nous vous invitons à les tester.

Exercices

• Prédiction de la vitesse maximale du vent de la station de Sprogø.
• Approximation de courbes non linéaires par la méthode des moindres carrés : application pour extraire des points de mesure sur des données topographiques LiDAR.
• Application en traitement d'images : comptage de bulles et de grains.
• Exemple de solution pour l'exercice de traitement d'images : comptage de grains dans un verre.

Merci aux personnes suivantes pour leur aide durant la traduction :

• sweetasnz
• Deusyss
• f-leb