[wiki: Calcul]

#!comment http://matplotlib.org/faq/usage_faq.html#matplotlib-pylab-and-pyplot-how-are-they-related

VisPy ? http://reseau-loops.github.io/journee_2014_06.html

Module `math`

Ce module fournit un ensemble de fonctions mathématiques pour les réels.

pi
sqrt
cos, sin, tan, acos, …
…

#!html
<br clear=right>

NumPy

Contenu

Tableaux multidimensionnels
Arithmétique et fonctions mathématiques sur les tableaux
Algèbre linéaire (LAPACK)
Transformée de Fourier (FFTPACK)
Nombres aléatoires
Outils pour intégrer du code Fortran, C/C++
Outils d’installation avec support tests unitaires

Remarques

Fonctions très rapides, proches des temps d’exécution en C sur les gros tableaux.
Code optimisé, de nombreux sous-modules disponibles.
Effort de vectorisation des algorithmes.
Tableaux multidimensionnels intégrés directement comme type (jusqu’à maintenant, seul le Fortran proposait cela).
Les matrices se comportent comme des tableaux, à part dans certaines opérations comme la multiplication, par exemple.
Distribueé au sein de !SciPy (voir ci-dessous).

Tableaux multidimensionnels

tous les élements du même type (entier, réel, …)
stockage compact des données, compatible C/Fortran
opérations efficaces
arithmétique
indexation flexible

Exemples d’opérations simples :

#!div style="float:left;width:48%"


{{{
#!python



author    = 'Simon CHOLLET (simon.chollet@ens.fr)'
date      = '18/08/2013'



import numpy



# Declaration des tableaux
int_array     = numpy.array([1, 2, 4, 8, 16])
float_array   = numpy.array([1.2, 2.4, 4.8, 8.16, 16.32])
complex_array = numpy.array(1, 2], [3, 4, dtype=complex)
array_2D = numpy.array(1, 2, 3],[4, 5, 6)
array_3D = numpy.array([1, 2, 3], [4, 5, 6, \
  4, 5, 6], [7, 8, 9,7, 8, 9], [1, 2, 3])



# Opérations sur tous les éléments
print int_array ** 2
print 10.0 * numpy.sin(float_array)



# Recuperation du type de valeurs
print int_array.dtype
print float_array.dtype
print complex_array.dtype



# Les dimensions de tableau
print "Dim(complex_array) = %d" % (complex_array.ndim)
print "Dim(array_2D)      = %d" % (array_2D.ndim)
print "Dim(array_3D)      = %d" % (array_3D.ndim)



# Le nombre d'elements
print "complex_array = %d elements." % (complex_array.size)
print "array_2D      = %d elements." % (array_2D.size)
print "array_3D      = %d elements." % (array_3D.size)



# Les tailles en octets des elements
print "Element int_array     = %d octets." % (int_array.itemsize)
print "Element float_array   = %d octets." % (float_array.itemsize)
print "Element complex_array = %d octets." % (complex_array.itemsize)



# Creation d'un tableau contenant que des 0
print numpy.zeros(10)
print numpy.zeros(10, dtype=numpy.int64)
print numpy.zeros((3, 5))



# Creation de tableaux contenant des 1
print numpy.ones(10)
print numpy.ones(5, dtype=numpy.complex)
print numpy.ones((4, 2))



# Creation de tableau vide
print numpy.empty(4)
print numpy.empty(4, dtype=numpy.int32)
print numpy.empty((2, 3))



# Creation d'un tableau allant de 5 en 5
print numpy.arange(10, 50, 5)



# Creation d'un tableau avec des valeurs regulierement espacees
print numpy.linspace(0, 9, 10)



# Combinaison de tableaux
print int_array - int_array
print int_array + int_array
print float_array * float_array
print float_array / float_array



# Multiplication de matrices
mat1 = numpy.array(1, 1], [0, 1)
mat2 = numpy.array(2, 0], [3, 4)
print mat1 * mat2
print numpy.dot(mat1, mat2)

}}}

#!div style="float:right;width:48%"


{{{
#!html
<pre style="margin-right:20px;border: 1px solid palegoldenrod;padding:4px;background-color:lightyellow">



[ 1 4 16 64 256 ]
[ 9.32039086 6.75463181 -9.96164609 9.5354065 -5.74535674 ]



int64
float64
complex128



Dim(complex_array) = 2
Dim(array_2D)      = 2
Dim(array_3D)      = 3



complex_array = 4 elements.
array_2D      = 6 elements.
array_3D      = 18 elements.



Element int_array     = 8 octets.
Element float_array   = 8 octets.
Element complex_array = 16 octets.



[ 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.



[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. ]
[ 1.+0.j 1.+0.j 1.+0.j 1.+0.j 1.+0.j ]
1. 1. ] ↵ [ 1. 1. ] ↵ [ 1. 1. ] ↵ [ 1. 1.



[ 6.93537601e-310 1.97731475e-316 1.75700415e-316 1.75700415e-316 ]
[ 39083840 0 39316480 0 ]
[[ 6.93537601e-310 1.94233015e-316 6.93537694e-310 ]
 [ 1.80609056e-316 1.92870896e-316 6.93537601e-310 ]]



[ 10 15 20 25 30 35 40 45 ]



[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. ]



[ 0 0 0 0 0 ]
[ 2 4 8 16 32 ]
[ 1.44 5.76 23.04 66.5856 266.3424 ]
[ 1. 1. 1. 1. 1. ]



2 0 ] ↵ [ 0 4
5 4 ] ↵ [ 3 4
</pre>

}}}

Quelques manipulations avancées :

#!div style="float:left;width:48%"


{{{
#!python



author    = 'Simon CHOLLET (simon.chollet@ens.fr)'
date      = '18/08/2013'



import numpy
import pylab



# 1 si nous utilisons la methode du module 'matplot' pour l'histogramme
MATPLOT_HISTOGRAM = 0



# Declaration des tableaux
int_array     = numpy.array([1, 2, 4, 8, 16])
float_array   = numpy.array([1.2, 2.4, 4.8, 8.16, 16.32])
complex_array = numpy.array(1, 2], [3, 4, dtype=complex)
array_2D = numpy.array(1, 2, 3],[4, 5, 6)
array_3D = numpy.array([1, 2, 3], [4, 5, 6, \
  4, 5, 6], [7, 8, 9,7, 8, 9], [1, 2, 3])



# Aplatissement des elements du tableau
for item in array_3D.flat:
    print item,
print



# Retrouve les dimensions d'un tableau
print array_3D.shape
print array_2D.shape



# Decoupage de tableau en plusieurs
rand_array = numpy.floor(10 * numpy.random.random((2, 12)))
pritn rand_array
print numpy.hsplit(rand_array, 3)
print numpy.vsplit(rand_array, 2)



# Donne la reference du tableau
print id(rand_array)



# Copie du tableau
int_array2 = int_array.copy()
print int_array2



# Teste les elements
bool_array = int_array > 5
print bool_array
print int_array[bool_array]



# Construction d'un vecteur de 10000 valeurs avec
# une variance de 0.5² et une moyenne de 2
mu, sigma = 2, 0.5
v = numpy.random.normal(mu, sigma, 10000)



# Trace un histogramme ou une courbe
if MATPLOT_HISTOGRAM:
    # Trace un histogramme
    pylab.hist(v, bins=50, normed=1)
    pylab.show()
else:
    # Trace une courbe
    (n, bins) = numpy.histogram(v, bins=50, normed=True)
    pylab.plot(.5 * (bins[1:] + bins[:-1]), n)
    pylab.show()

}}}

#!div style="float:right;width:48%"


{{{
#!html
<pre style="margin-right:20px;border: 1px solid palegoldenrod;padding:4px;background-color:lightyellow">



1 2 3 4 5 6 4 5 6 7 8 9 7 8 9 1 2 3



(3L, 2L, 3L)
(2L, 3L)



[[ 3., 7., 1., 1., 2., 5., 3., 9., 5., 1., 0., 1. ]
 [ 4., 0., 8., 5., 5., 7., 1., 8., 0., 6., 5., 6. ]]



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



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



[ False False False  True  True ]
[ 8 16 ]

}}}

#!html
<br clear=right>

Un exemple plus complexe :

author    = 'Simon CHOLLET (simon.chollet@ens.fr)'
date      = '18/08/2013'


import numpy



def mandelbrot(h, w, maxit=50):
    '''
    Fonction qui retourne une image de fractale de Mandelbrot
    de taille (h, w).
    -------------------------------------------------------------------------
    Arguments :
        - h     : Hauteur de l'image a generer.
        - w     : Largeur de l'image a generer.
        - maxit : Finesse.
    -------------------------------------------------------------------------
    Retour :
        Image fractale.
    '''
    y, x = numpy.ogrid[-1.4:1.4:h * 1j, -2:0.8:w * 1j]
    c = x + y * 1j
    z = c
    divtime = maxit + numpy.zeros(z.shape, dtype=numpy.int)



for i in xrange(maxit):
        z  = z ** 2 + c
        diverge = z * numpy.conj(z) > 2 ** 2
        div_now = diverge & (divtime == maxit)
        divtime[div_now] = i
        z[diverge] = 2



return divtime



if name == "main":



import pylab



MANDELBROT_NB_POINTS = 50
    # 0 = version vectorisee
    # 1 = version non vectorisee
    # 2 = version par fonction
    MANDELBROT_VERSION = 0



# Calculs preliminaires
    divergence = numpy.zeros((MANDELBROT_NB_POINTS, MANDELBROT_NB_POINTS))
    if MADELBROT_VERSION == 0:
        # Version non vectorisee
        for c_x in numpy.linspace(-1.5, 0.5, MANDELBROT_NB_POINTS):
            for c_y in numpy.linspace(-1, 1, MANDELBROT_NB_POINTS):
                c = c_x + c_y * 1j
                z = 0
                for i in range(50):
                    z = z ** 2 + c
                    if numpy.abs(z) > 10:
                        divergence[(c_x + 1.49) * MANDELBROT_NB_POINTS / 2,
                                   (c_y + 0.99) * MANDELBROT_NB_POINTS / 2] = 50 - i
                        break
    elif MANDELBROT_VERSION == 1:
        # Version vectorisee
        c_x, c_y = numpy.ogrid[-1.5:0.5:MANDELBROT_NB_POINTS * 1j,
                               -1:1:MANDELBROT_NB_POINTS * 1j]
        c = c_x + c_y * 1j
        divergence = numpy.zeros((MANDELBROT_NB_POINTS, MANDELBROT_NB_POINTS))
        z = numpy.zeros((MANDELBROT_NB_POINTS, MANDELBROT_NB_POINTS))
        masque = numpy.ones((MANDELBROT_NB_POINTS, MANDELBROT_NB_POINTS),
                            dtype=numpy.bool)
        for i in range(50):
            z[masque] = z[masque] ** 2 + c[masque]
            masque = (numpy.abs(z) > 10)
            divergence += masque



# Trace du resultat
    if MANDELBROT_VERSION == 0 or MANDELBROT_VERSION == 1:
        pylab.imshow(divergence, cmap=pylab.cm.spectral, extent=(-1, 1, -1, 1))
    else:
        pylab.imshow(mandelbrot(MANDELBROT_NB_POINTS, MANDELBROT_NB_POINTS))
    pylab.show()

#!comment
Comment optimiser !NumPy :


utiliser l’API C directement,
utiliser Swig,
utiliser f2py,
utiliser cython,...

#!html
<br clear=right>

Matplotlib

Bibliothèque de visualisation scientifique, conçue comme une alternative à Matlab.

Contenu

Plots 2D
Plots 3D avec l’extension mplot3d
Affichage à l’écran, EPS, PDF, …
Intégrable dans un interface graphique

PyLab : distribué avec Matplotlib, ce module intègre dans le même espace de nom tout numpy et tout matplotlib.pyplot, afin d’offrir à l’utilisateur un environnement le plus proche possible de MatLab. Voir la FAQ.

#!html
<br clear=right>

SymPy

Pour le calcul symbolique.

#!html
<br clear=right>

SciPy

#!comment

PySide ? http://ipython.org/ipython-doc/stable/interactive/qtconsole.html ? Matplotlib ? http://www.pyqtgraph.org/ ?

Description :

Comme numpy, fournit des fonctions mathématiques de base sur tableaux et matrices.
Intègre des algorithmes numériques : résolution d’équations, calculs statistiques, etc.
Il est vaste et grandit régulièrement.
Point d’entrée « unique » pour les scientifiques.
Utilise quasiment les mêmes méthodes que numpy.
Aussi optimisé en C.

Les composants principaux de la Distribution !SciPy :

Python
!NumPy
Bibliothèque !SciPy
Matplotlib
Pandas
!SimPy
IPython
Nose

Quelques paquets de la Bibliothèque !SciPy :

fft : Transformées de Fourrier.
integrate : Fonctions d’intégration numérique.
interpolate : Fonctions d’interpolation, linéaire, cubique, etc.
io : Fonctions d’entrées / sorties vers les fichiers.
linalg : Fonctions sur les fonctions d’algèbre linéaire.
ndimage : Fonctions de traitement d’images.
optimize : Fonctions d’optimisation de fonction.
sparse : Fonctions sur les matrices.
special : Série de fonctions spéciales.
stats : Fonctions pour les calculs statistiques.
tests : Fonctions qui permettent de tester la librairie, elle-même.
utils : Fonctions utilitaires diverses.

SciKits

scikit-aero Aeronautical engineering calculations in Python.
scikit-commpy Digital Communication Algorithms with Python
scikit-fmm An extension module implimenting the fast marching method
scikit-image Image processing routines for !SciPy
scikit-learn A set of python modules for machine learning and data mining
scikit-monaco Python modules for Monte Carlo integration
scikit-nano Python toolkit for generating and analyzing nanostructure data
scikit-rf Open Source RF Engineering
scikit-tensor Python module for multilinear algebra and tensor factorizations
scikit-vi Scikit providing Virtual Instruments
scikits-image Image processing routines for SciPy
ann Approximate Nearest Neighbor library wrapper for Numpy
audiolab A python module to make noise from numpy arrays
bootstrap Bootstrap confidence interval estimation routines for !SciPy
bvp1lg Boundary value problem (legacy) solvers for ODEs
bvp_solver Python package for solving two-point boundary value problems
cuda Python interface to GPU-powered libraries
datasmooth Scikits data smoothing package
eartho Earth Observation routines for !SciPy
example Scikits example package
fitting Framework for fitting functions to data with !SciPy
hydroclimpy Environmental time series manipulation
learn A set of python modules for machine learning and data mining
odes A python module for ordinary differential equation anddifferential algebraic equation solvers
optimization A python module for numerical optimization
samplerate A python module for high quality audio resampling
scattpy Light Scattering methods for Python
sparse Scikits sparse matrix package
statsmodels Statistical computations and models for use with !SciPy
talkbox Talkbox, a set of python modules for speech/signal processing
timeseries Time series manipulation
vectorplot Vector fields plotting algorithms.

Exemples :

#!div style="float:left;width:48%"


{{{
#!python



author    = 'Simon CHOLLET (simon.chollet@ens.fr)'
date      = '18/08/2013'



import scipy
import scipy.special
import scipy.interpolate
import scipy.linalg
import pylab



# Declaration des tableaux
int_array   = scipy.array([1, 2, 4, 8, 16])
float_array = scipy.array([1.2, 2.4, 4.8, 8.16, 16.32])
array_2D = scipy.array([[1, 2, 3],
                        [4, 5, 6]])
array_3D = scipy.array([1, 2, 3], [4, 5, 6,
                        4, 5, 6], [7, 8, 9,
                        7, 8, 9], [1, 2, 3])



tab = scipy.array(1, 2], [3, 4)
mat = scipy.mat(tab)
print "* Valeurs propres mat = ", scipy.linalg.eigvals(mat)
print "* Determinant mat = ", scipy.linalg.det(mat)
print "* Matrice inverse = \n", mat.I



# Calcul de FFT
rate = 1000.0
t = scipy.r_[0:0.6:1 / rate]
nbT = len(t)



# Le signal a analyser
s = scipy.sin(2.0 * scipy.pi * 50.0 * t) + \
    scipy.sin(2.0 * scipy.pi * 70.0 * t + scipy.pi / 4.0)



# Calcul de la FFT
fftS = scipy.fft(s)



# Calcul des limites de graphiques
f = rate * scipy.r_[0:(nbT / 2)] / nbT
n = len(f)
pylab.plot(f, abs(fftS[0:n]) / nbT)
pylab.show()

}}}

#!div style="float:right;width:48%"


{{{
#!html
<pre style="margin-right:20px;border: 1px solid palegoldenrod;padding:4px;background-color:lightyellow">



Valeurs propres mat =  [-0.37228132+0.j  5.37228132+0.j]
Determinant mat =  -2.0
[[-2.   1. ]
 [ 1.5 -0.5]]

}}}

#!html
<br clear=right>

Intégration et interpolation :

#!div style="float:left;width:48%"


{{{
#!python



author    = 'Simon CHOLLET (simon.chollet@ens.fr)'
date      = '18/08/2013'



import scipy
import scipy.special
import scipy.integrate
import scipy.interpolate
import scipy.linalg
import pylab



def f(x):
    return x  4 + x  3 - x ** 2 - 10



# Calcul de l'intergrale de x^5 de 0 a 1
value, err = scipy.integrate.quad(func=pow, a=0., b=1., args=(5,))
print "* Integrale de x^5, entre 0 et 1 = %f" % (value)



# Calcul de l'intergrale de f(x) de 0 a 1
value, err = scipy.integrate.quad(func=f, a=0., b=10., args=())
print "* Integrale de f(x), entre 0 et 10 = %f" % (value)



# Interpolation
# Axe X
x = scipy.linspace(0, 1, 10)
# Axe Y
y = scipy.sin(2 * scipy.pi * x)
# Interpolation lineaire
linear_interp = scipy.interpolate.interp1d(x, y)
# Interpolation cubique
cubic_interp = scipy.interpolate.interp1d(x, y, kind='cubic')



# Affichage des courbes
fine_x = scipy.linspace(0, 1, 50)
y_fromlinearinterp = linear_interp(fine_x)
y_fromcubicinterp = cubic_interp(fine_x)
pylab.plot(fine_x, scipy.sin(2 * scipy.pi * fine_x), 'b-',
           fine_x, y_fromlinearinterp, 'g-^',
           fine_x, y_fromcubicinterp, 'r-.')
pylab.show()

}}}

#!div style="float:right;width:48%"


{{{
#!html
<pre style="margin-right:20px;border: 1px solid palegoldenrod;padding:4px;background-color:lightyellow">



Integrale de x^5 = 0.166667



Integrale de x^5 = 0.166667

}}}

#!html
<br clear=right>

Traitement d’image :

#!python


author    = 'Simon CHOLLET (simon.chollet@ens.fr)'
date      = '18/08/2013'



import scipy
import scipy.ndimage
import pylab



img = scipy.misc.lena()
imgFloue  = scipy.ndimage.gaussian_filter(img, 1)
imgRotate = scipy.ndimage.rotate(img, 45)
	
#pylab.imshow(img)
#pylab.imshow(imgFloue)
pylab.imshow(imgRotate)
pylab.show()

VPython

VPython regroupe le langage de programmation Python et le module graphique 3D visual, développé par David Scherer. Il repose également sur !NumPy et wxPython (interface utilisateur).
Exemple de quelques lignes : faisons tourner une boule grise sur un anneau, qui serait dessiné autour d’une plus grosse boule bleue :

#!python


author    = 'Simon CHOLLET (simon.chollet@ens.fr)'
date      = '17/08/2013'



import visual
import math



# Creation d'une grosse sphere
s1 = visual.sphere()
# Changement de sa couleur : bleue
s1.color = (0, 0, 1)
# Creation d'un anneau
r = visual.ring(radius=2)
# Creation d'une petite sphere de couleur verte
s2 = visual.sphere(pos=(0, 0, 2), color=(0, 1, 0), radius=0.7)
# Mise en scene des objets
visual.scene.autoscale = 0
# Animation ... tout le temps
t = 0
while True:
    # Temporisation de 30 images par seconde
    visual.rate(30)
    # Changement de la position de la petite sphere
    s2.pos = (0,
        2 * math.cos(t * math.pi / 30),
        2 * math.sin(t * math.pi / 30))
    # Increment du temps
    t += 1

RPy

Interface Python pour le langage de programmation R, dédié aux traitements et analyses statistiques.

Implémentation reposant sur le module numpy.
Fonctionnalités intéressantes : traiter des séries de données, statistiques, dessiner des images au format PNG, etc.
Mais aussi : unique(), sort(), table(), barplot(), etc.
ATTENTION : ne fonctionne actuellement que sous Linux, un peu plus dur sous Windows …

Exemple de code :

#!python


import random
from rpy import r as R



#==========
# PARAMETRES
#==========



# Nombre d'elements a generer dans les exemples
RPY_NB_ITEMS = 500
# Fichier image a generer
RPY_OUTPUT_IMAGE_1 = '../output/code_rpy1.png'
RPY_OUTPUT_IMAGE_2 = '../output/code_rpy2.png'



#==========
# GENERATION ALEATOIRE D'ELEMENTS
#==========



# Construction d'une liste de plusieurs elements croissants
x = range(1, RPY_NB_ITEMS + 1)



# Construction d'une liste de plusieurs elements aleatoires
y = []
for i in x:
    y.append(random.random())



# Enregistrement du graphique dans un fichier PNG
print '* Creation du fichier Image : %s' % (RPY_OUTPUT_IMAGE_1)
R.png(RPY_OUTPUT_IMAGE_1)



# Dessine les points sur l'image
R.plot(x, y, xlab="Position", ylab="Coordonnées", col="black", pch=3)



# Fin du graphique
R.dev_off()



#==========
# GENERATION D'HISTOGRAMME
#==========



# La sequence de donnees a analyser
dataList = [1, 2, 3, 1, 2, 3, 4, 5, 5, 6, 7, 3, 1,
             8, 4, 6, 9, 9, 2, 1, 5, 6, 8, 9, 4, 3,
             3, 2, 7, 5, 9, 7, 5, 6, 4, 3, 3, 2, 1]
print '* Donnees a analyser :', dataList



# Enregistrement du graphique dans un fichier PNG
print '* Creation du fichier Image : %s' % (RPY_OUTPUT_IMAGE_2)
R.png(RPY_OUTPUT_IMAGE_2)



# Tri des donnees
dataSort = R.sort(dataList)
print '* Donnees triees :', dataSort



# Extrait les bases
bases = R.unique(dataSort)
print '* Bases pour histogramme :', bases



# Effectif de chaque base
effectifs = R.table(dataSort)
print '* Effectifs dans chaque classe :', effectifs



# Dessine l'histogramme et sauvegarde de la position des abscisses
coords = R.barplot(effectifs, ylab="Nombre")



# Ajout du texte pour l'axe des abscisses
R.text(coords, -0.5, bases, xpd=True, cex=1.5, font=2)



# Fin du graphique
R.dev_off()

#!comment


== ScientificPython ==

Conclusions

Le module numpy est largement utilisé.
Scipy est un très gros module.
Nous pouvons le compléter avec des ‘!SciKits’ …
Nous avons parcouru qu’une infime partie des possibilités de ces modules.
Il existe sûrement des algorithmes et exemples qui peuvent aider à résoudre certains de vos problèmes spécifiques.
De nombreux autres domaines sont traités : aéronautique, simulation, biologie, chimie, médical, etc.
Vous pouvez aussi y apporter votre pierre avec notamment l’interfaçage de vos propres bibliothèques de calculs.

#!html
<br clear=right>

Wiki Electronique

BibliosCalcul

[wiki: Calcul]

Module `math`

NumPy

Matplotlib

SymPy

SciPy

VPython

RPy

Conclusions

References

[wiki: Calcul]

Module math

NumPy

Matplotlib

SymPy

SciPy

VPython

RPy

Conclusions

References

Module `math`