Différentiation automatique et Python

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()