Voici un bout de code Python qui illustre le concept de différentiation automatique par surcharge d’opérateur. Le programme ci-dessous comporte des méthodes de test. On peut donc directement le donner à l’interpréteur. Pour plus d’informations.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy
import pylab
class adarray():
def __init__(self,v):
self.v = v
def __repr__(self):
return str(adarray(self.v))
def __str__(self):
return str(self.v)
def __add__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
return adarray(self.v+other.v)
def __radd__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
return adarray(self.v+other.v)
def __sub__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
return adarray(self.v-other.v)
def __rsub__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
return adarray(other.v-self.v)
def __mul__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=self.v[:,0]*other.v[:,0]
for i in range(1,self.v.shape[1]):
result[:,i]=other.v[:,0]*self.v[:,i]+self.v[:,0]*other.v[:,i]
return adarray(result)
def __rmul__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=self.v[:,0]*other.v[:,0]
for i in range(1,self.v.shape[1]):
result[:,i]=other.v[:,0]*self.v[:,i]+self.v[:,0]*other.v[:,i]
return adarray(result)
def __div__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=self.v[:,0]/other.v[:,0]
o2=other.v[:,0]*other.v[:,0]
for i in range(1,self.v.shape[1]):
result[:,i]=(other.v[:,0]*self.v[:,i]-self.v[:,0]*other.v[:,i])/o2
return adarray(result)
def __rdiv__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=other.v[:,0]/self.v[:,0]
o2=self.v[:,0]*self.v[:,0]
for i in range(1,self.v.shape[1]):
result[:,i]=(other.v[:,i]*self.v[:,0]-self.v[:,i]*other.v[:,0])/o2
return adarray(result)
def __pow__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=pow(self.v[:,0],other.v[:,0])
for i in range(1,self.v.shape[1]):
result[:,i]=(other.v[:,i]*log(self.v[:,i])+other.v[:,0]*self.v[:,i]/self.v[:,0])*result[:,0]
return adarray(result)
def __rpow__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=pow(other.v[:,0],self.v[:,0])
for i in range(1,self.v.shape[1]):
result[:,i]=(self.v[:,i]*log(other.v[:,i])+self.v[:,0]*other.v[:,i]/other.v[:,0])*result[:,0]
return adarray(result)
def sin(x):
if not isinstance(x,adarray) :
return numpy.sin(x)
else :
result=numpy.zeros(x.v.shape)
result[:,0]=numpy.sin(x.v[:,0])
for i in range(1,x.v.shape[1]):
result[:,i]=x.v[:,i]*numpy.cos(x.v[:,0])
return adarray(result)
def cos(x):
if not isinstance(x,adarray) :
return numpy.sin(x)
else :
result=numpy.zeros(x.v.shape)
result[:,0]=numpy.cos(x.v[:,0])
for i in range(1,x.v.shape[1]):
result[:,i]=-x.v[:,i]*numpy.sin(x.v[:,0])
return adarray(result)
def exp(x):
if not isinstance(x,adarray) :
return numpy.exp(x)
else :
result=numpy.zeros(x.v.shape)
result[:,0]=numpy.exp(x.v[:,0])
for i in range(1,x.v.shape[1]):
result[:,i]=x.v[:,i]*numpy.exp(x.v[:,0])
return adarray(result)
def log(x):
if not isinstance(x,adarray) :
return numpy.log(x)
else :
result=numpy.zeros(x.v.shape)
result[:,0]=numpy.log(x.v[:,0])
for i in range(1,x.v.shape[1]):
result[:,i]=x.v[:,i]/x.v[:,0]
return adarray(result)
def CreateVar(base,calrange):
nbvar=calrange.shape[1]
for i in range(nbvar):
exec ("global %s%s" % (base,i))
exec ("globals()['%s%s'] = adarray(numpy.zeros([calrange.shape[0],calrange.shape[1]+1]))" % (base,i))
exec ("%s%s.v[:,0] = calrange[:,i]" % (base,i))
exec ("%s%s.v[:,i+1] = numpy.ones(calrange.shape[0])" % (base,i))
#exec ("print %s%s" % (base,i))
def testMonoDim():
### Cas monodimensionnel
plagex0=numpy.arange(0.001, 4., 0.1)
plage=numpy.array([plagex0]).T
CreateVar('x',plage)
## Opérations avec les scalaires
print "Opérations avec les scalaires :\n"
print "Matrice initiale : x = \n"+str(x0)
# Addition avec un scalaire après
print "Addition avec un scalaire après, \t x + 5 = \n%s" % (x0+5)
# Addition avec un scalaire avant
print "Addition avec un scalaire avant, \t 5 + x = \n%s" % (5+x0)
# Soustraction avec un scalaire après
print "Soustraction avec un scalaire après, \t x - 5 = \n%s" % (x0-5)
# Soustraction avec un scalaire avant
print "Soustraction avec un scalaire avant, \t 5 - x = \n%s" % (5-x0)
# Multiplication avec un scalaire après
print "Multiplication avec un scalaire après, \t x * 5 = \n%s" % (x0*5)
# Multiplication avec un scalaire avant
print "Multiplication avec un scalaire avant, \t 5 * x = \n%s" % (5*x0)
# Division avec un scalaire après
print "Division avec un scalaire après, \t x / 3 = \n%s" % (x0/3)
# Division avec un scalaire avant
print "Division avec un scalaire avant, \t 3 / x = \n%s" % (3/x0)
# Puissance avec un scalaire
print "Puissance avec un scalaire, \t x^2 = \n%s" % (pow(x0,2))
# Puissance avec un scalaire
print "Puissance avec un scalaire, \t x^3 = \n%s" % (pow(x0,3))
# Puissance avec un scalaire
print "Puissance avec un scalaire, \t 2.^x = \n%s" % (pow(2,x0))
# Puissance avec un scalaire
print "Puissance avec un scalaire, \t 3.^x = \n%s" % (pow(3,x0))
## Opérations avec la variable
print "\nOpérations avec la variable :\n"
print "Vecteur initial : x = \n"+str(x0)
# Addition avec un pair
print "Addition,\t x + x = \n%s" % (x0+x0)
# Soustraction avec un pair
print "Soustraction,\t x - x = \n%s" % (x0-x0)
# Multiplication avec un pair
print "Multiplication,\t x * x = \n%s" % (x0*x0)
# Division avec un pair
print "Division,\t x / x = \n%s" % (x0/x0)
# Cosinus avec un pair
print "Cosinus,\t cos(x) = \n%s" % (cos(x0))
# Sinus avec un pair
print "Sinus, \t sin(x) = \n%s" % (sin(x0))
# Log avec un pair
print "Log, \t log(x) = \n%s" % (log(x0))
# Exp avec un pair
print "Exp, \t exp(x) = \n%s" % (exp(x0))
# Puissance avec un pair
print "Puissance avec un pair, \t x^x = \n%s" % (pow(x0,x0))
# Puissance avec un pair
print "Puissance avec un pair, \t x^(3*x) = \n%s" % (pow(x0,3*x0))
## Opérations complexes
print "\nOpérations complexes :\n"
print "Vecteur initial : x = \n"+str(x0)
# Addition
print "Addition,\t 10*x + 3*x = \n%s" % (10*x0+3*x0)
# Soustraction
print "Soustraction,\t 4*x - 5*x = \n%s" % (4*x0-5*x0)
# Multiplication
print "Multiplication,\t 3*x*4*x = \n%s" % (3*x0*4*x0)
# Division
print "Division,\t x/(x+1) = \n%s" % (x0/(x0+1))
# Cosinus
print "Cosinus,\t cos(pi*x/2-pi/2) = \n%s" % (cos(numpy.pi*x0/2-numpy.pi/2))
# Sinus
print "Sinus, \t sin(pi*x/2-pi/2) = \n%s" % (sin(numpy.pi*x0/2-numpy.pi/2))
# Log
print "Log, \t log(10*x)/log(10) = \n%s" % (log(10*x0)/log(10))
# Exp
print "Exp, \t x*exp(x/2) = \n%s" % (x0*exp(x0/2))
# Formule
print "Exp, \t x^5*exp(x/2)*sin(x0) = \n%s" % (pow(x0,5)*sin(x0)*exp(x0/2))
def testMonoDimGraph():
## Graphique
min=0.001
max=numpy.pi
nbpts=300
step=(max-min)/nbpts
print "Le pas vaut %s" %step
t=numpy.arange(min,max,step)
plage=numpy.array([t]).T
res=numpy.zeros(nbpts)
dres=numpy.zeros(nbpts)
CreateVar('x',plage)
y=sin(x0)*5*x0-1/(1+x0*x0)
pylab.figure(1)
pylab.plot(t,y.v[:,0])
pylab.plot(t,y.v[:,1])
pylab.show()
def testBiDim():
### Cas bidimensionnel
nbpts=6.
minx=1.
maxx=3.
rx=numpy.arange(minx, maxx, (maxx-minx)/nbpts)
miny=0.5
maxy=2.33
ry=numpy.arange(miny, maxy, (maxy-miny)/nbpts)
plage=numpy.array([rx,ry]).T
CreateVar('x',plage)
## Opérations avec les scalaires
print "Opérations avec les scalaires :\n"
print "Matrice initiale : x = \n"+str(x0)
print "Matrice initiale : y = \n"+str(x1)
# Addition avec un scalaire après
print "Addition avec un scalaire après, \t x + y + 5 = \n%s" % (x0+x1+5)
# Addition avec un scalaire avant
print "Addition avec un scalaire avant, \t 5 + x + y = \n%s" % (5+x0+x1)
# Soustraction avec un scalaire après
print "Soustraction avec un scalaire après, \t x + y - 5 = \n%s" % (x0+x1-5)
# Soustraction avec un scalaire avant
print "Soustraction avec un scalaire avant, \t 5 - x -y = \n%s" % (5-x0-x1)
# Multiplication avec un scalaire après
print "Multiplication avec un scalaire après, \t x * y * 5 = \n%s" % (x0*x1*5)
# Multiplication avec un scalaire avant
print "Multiplication avec un scalaire avant, \t 5 * x * y = \n%s" % (5*x0*x1)
# Division avec un scalaire après
print "Division avec un scalaire après, \t x * y / 3 = \n%s" % (x0*x1/3)
# Division avec un scalaire avant
print "Division avec un scalaire avant, \t 3 / (x * y) = \n%s" % (3/(x0*x1))
## Opérations complexes
print "\nOpérations complexes :\n"
print "Vecteur initial : x = \n"+str(x0)
print "Vecteur initial : y = \n"+str(x1)
# Addition
print "Addition,\t 10*x + 3*y = \n%s" % (10*x0+3*x1)
# Soustraction
print "Soustraction,\t 4*x - 5*y = \n%s" % (4*x0-5*x1)
# Multiplication
print "Multiplication,\t 3*x*4*y = \n%s" % (3*x0*4*x1)
# Division
print "Division,\t x/(y+1) = \n%s" % (x0/(x1+1))
# Cosinus
print "Cosinus,\t cos(pi*x*y/2) = \n%s" % (cos(numpy.pi*x0*x1/2))
# Sinus
print "Sinus, \t sin(pi*x*y/2) = \n%s" % (sin(numpy.pi*x0*x1/2))
# Log
print "Log, \t log(10*x*y)/log(10) = \n%s" % (log(10*x0*x1)/log(10))
# Exp
print "Exp, \t x*y*exp(x*y/2) = \n%s" % (x0*x1*exp(x0*x1/2))
# Complexe
print "Formule \t x*y*cos(x)-y*log(x*y) = \n%s" % (x0*x1*cos(x0)-x1*log(x0*x1))
testMonoDim()
testMonoDimGraph()
testBiDim()
