src/daComposant/daNumerics/ApproximatedDerivatives.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2012 EDF R&D
   4 #
   5 #  This library is free software; you can redistribute it and/or
   6 #  modify it under the terms of the GNU Lesser General Public
   7 #  License as published by the Free Software Foundation; either
   8 #  version 2.1 of the License.
   9 #
  10 #  This library is distributed in the hope that it will be useful,
  11 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 #  Lesser General Public License for more details.
  14 #
  15 #  You should have received a copy of the GNU Lesser General Public
  16 #  License along with this library; if not, write to the Free Software
  17 #  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 #  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 #  Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 __doc__ = """
  24     Définit les versions approximées des opérateurs tangents et adjoints.
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import os, numpy, time
  29 import logging
  30 # logging.getLogger().setLevel(logging.DEBUG)
  31
  32 # ==============================================================================
  33 class FDApproximation:
  34     """
  35     Cette classe sert d'interface pour définir les opérateurs approximés. A la
  36     création d'un objet, en fournissant une fonction "Function", on obtient un
  37     objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
  38     "AdjointOperator". On contrôle l'approximation DF avec l'incrément
  39     multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
  40     "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
  41     centrées si le booléen "centeredDF" est vrai.
  42     """
  43     def __init__(self, Function = None, centeredDF = False, increment = 0.01, dX = None):
  44         self.__userFunction = Function
  45         self.__centeredDF = bool(centeredDF)
  46         if float(increment) <> 0.:
  47             self.__increment  = float(increment)
  48         else:
  49             self.__increment  = 0.01
  50         if dX is None:
  51             self.__dX     = None
  52         else:
  53             self.__dX     = numpy.asmatrix(numpy.ravel( dX )).T
  54
  55     # ---------------------------------------------------------
  56     def DirectOperator(self, X ):
  57         """
  58         Calcul du direct à l'aide de la fonction fournie.
  59         """
  60         _X = numpy.asmatrix(numpy.ravel( X )).T
  61         _HX  = self.__userFunction( _X )
  62         return numpy.ravel( _HX )
  63
  64     # ---------------------------------------------------------
  65     def TangentMatrix(self, X ):
  66         """
  67         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
  68         c'est-à-dire le gradient de H en X. On utilise des différences finies
  69         directionnelles autour du point X. X est un numpy.matrix.
  70
  71         Différences finies centrées :
  72         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
  73            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
  74            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
  75            H( X_moins_dXi )
  76         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
  77            le pas 2*dXi
  78         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
  79
  80         Différences finies non centrées :
  81         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
  82            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
  83            HX_plus_dXi = H( X_plus_dXi )
  84         2/ On calcule la valeur centrale HX = H(X)
  85         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
  86            le pas dXi
  87         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
  88
  89         """
  90         logging.debug("  == Calcul de la Jacobienne")
  91         logging.debug("     Incrément de............: %s*X"%float(self.__increment))
  92         logging.debug("     Approximation centrée...: %s"%(self.__centeredDF))
  93         #
  94         _X = numpy.asmatrix(numpy.ravel( X )).T
  95         #
  96         if self.__dX is None:
  97             _dX  = self.__increment * _X
  98         else:
  99             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
 100         #
 101         if (_dX == 0.).any():
 102             moyenne = _dX.mean()
 103             if moyenne == 0.:
 104                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 105             else:
 106                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 107         #
 108         if self.__centeredDF:
 109             #
 110             # Boucle de calcul des colonnes de la Jacobienne
 111             # ----------------------------------------------
 112             _Jacobienne  = []
 113             for i in range( len(_dX) ):
 114                 _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 115                 _X_plus_dXi[i]  = _X[i] + _dX[i]
 116                 _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 117                 _X_moins_dXi[i] = _X[i] - _dX[i]
 118                 #
 119                 _HX_plus_dXi  = self.DirectOperator( _X_plus_dXi )
 120                 _HX_moins_dXi = self.DirectOperator( _X_moins_dXi )
 121                 #
 122                 _HX_Diff = numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dX[i])
 123                 #
 124                 _Jacobienne.append( _HX_Diff )
 125             #
 126         else:
 127             #
 128             # Boucle de calcul des colonnes de la Jacobienne
 129             # ----------------------------------------------
 130             _HX_plus_dX = []
 131             for i in range( len(_dX) ):
 132                 _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 133                 _X_plus_dXi[i] = _X[i] + _dX[i]
 134                 #
 135                 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 136                 #
 137                 _HX_plus_dX.append( _HX_plus_dXi )
 138             #
 139             # Calcul de la valeur centrale
 140             # ----------------------------
 141             _HX = self.DirectOperator( _X )
 142             #
 143             # Calcul effectif de la Jacobienne par différences finies
 144             # -------------------------------------------------------
 145             _Jacobienne = []
 146             for i in range( len(_dX) ):
 147                 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 148         #
 149         _Jacobienne = numpy.matrix( numpy.vstack( _Jacobienne ) ).T
 150         logging.debug("  == Fin du calcul de la Jacobienne")
 151         #
 152         return _Jacobienne
 153
 154     # ---------------------------------------------------------
 155     def TangentOperator(self, (X, dX) ):
 156         """
 157         Calcul du tangent à l'aide de la Jacobienne.
 158         """
 159         _Jacobienne = self.TangentMatrix( X )
 160         if dX is None or len(dX) == 0:
 161             #
 162             # Calcul de la forme matricielle si le second argument est None
 163             # -------------------------------------------------------------
 164             return _Jacobienne
 165         else:
 166             #
 167             # Calcul de la valeur linéarisée de H en X appliqué à dX
 168             # ------------------------------------------------------
 169             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 170             _HtX = numpy.dot(_Jacobienne, _dX)
 171             return _HtX.A1
 172
 173     # ---------------------------------------------------------
 174     def AdjointOperator(self, (X, Y) ):
 175         """
 176         Calcul de l'adjoint à l'aide de la Jacobienne.
 177         """
 178         _JacobienneT = self.TangentMatrix( X ).T
 179         if Y is None or len(Y) == 0:
 180             #
 181             # Calcul de la forme matricielle si le second argument est None
 182             # -------------------------------------------------------------
 183             return _JacobienneT
 184         else:
 185             #
 186             # Calcul de la valeur de l'adjoint en X appliqué à Y
 187             # --------------------------------------------------
 188             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 189             _HaY = numpy.dot(_JacobienneT, _Y)
 190             return _HaY.A1
 191
 192 # ==============================================================================
 193 #
 194 def test1( XX ):
 195     """ Direct non-linear simulation operator """
 196     #
 197     # NEED TO BE COMPLETED
 198     # NEED TO BE COMPLETED
 199     # NEED TO BE COMPLETED
 200     #
 201     # --------------------------------------> # EXAMPLE TO BE REMOVED
 202     # Example of Identity operator            # EXAMPLE TO BE REMOVED
 203     if type(XX) is type(numpy.matrix([])):    # EXAMPLE TO BE REMOVED
 204         HX = XX.A1.tolist()                   # EXAMPLE TO BE REMOVED
 205     elif type(XX) is type(numpy.array([])):   # EXAMPLE TO BE REMOVED
 206         HX = numpy.matrix(XX).A1.tolist()     # EXAMPLE TO BE REMOVED
 207     else:                                     # EXAMPLE TO BE REMOVED
 208         HX = XX                               # EXAMPLE TO BE REMOVED
 209     #                                         # EXAMPLE TO BE REMOVED
 210     HHX = []                                  # EXAMPLE TO BE REMOVED
 211     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 212     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 213     # --------------------------------------> # EXAMPLE TO BE REMOVED
 214     #
 215     return numpy.array( HHX )
 216
 217 # ==============================================================================
 218 if __name__ == "__main__":
 219
 220     print
 221     print "AUTODIAGNOSTIC"
 222     print "=============="
 223
 224     X0 = [1, 2, 3]
 225
 226     FDA = FDApproximation( test1 )
 227     print "H(X)       =",   FDA.DirectOperator( X0 )
 228     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 229     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 230     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
 231     print
 232     del FDA
 233     FDA = FDApproximation( test1, centeredDF=True )
 234     print "H(X)       =",   FDA.DirectOperator( X0 )
 235     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 236     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 237     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
 238     print
 239     del FDA
 240
 241     print "=============="
 242     print
 243     X0 = range(5)
 244
 245     FDA = FDApproximation( test1 )
 246     print "H(X)       =",   FDA.DirectOperator( X0 )
 247     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 248     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 249     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 250     print
 251     del FDA
 252     FDA = FDApproximation( test1, centeredDF=True )
 253     print "H(X)       =",   FDA.DirectOperator( X0 )
 254     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 255     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 256     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 257     print
 258     del FDA
 259
 260     print "=============="
 261     print
 262     X0 = numpy.arange(3)
 263
 264     FDA = FDApproximation( test1 )
 265     print "H(X)       =",   FDA.DirectOperator( X0 )
 266     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 267     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 268     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 269     print
 270     del FDA
 271     FDA = FDApproximation( test1, centeredDF=True )
 272     print "H(X)       =",   FDA.DirectOperator( X0 )
 273     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 274     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 275     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 276     print
 277     del FDA
 278
 279     print "=============="
 280     print
 281     X0 = numpy.asmatrix(numpy.arange(4)).T
 282
 283     FDA = FDApproximation( test1 )
 284     print "H(X)       =",   FDA.DirectOperator( X0 )
 285     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 286     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 287     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 288     print
 289     del FDA
 290     FDA = FDApproximation( test1, centeredDF=True )
 291     print "H(X)       =",   FDA.DirectOperator( X0 )
 292     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 293     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 294     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 295     print
 296     del FDA
 297