src/daComposant/daNumerics/ApproximatedDerivatives.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2013 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 (approximation d'ordre 2):
  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 (approximation d'ordre 1):
  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("FDA Calcul de la Jacobienne")
  91         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
  92         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
  93         #
  94         if X is None or len(X)==0:
  95             raise ValueError("Nominal point X for approximate derivatives can not be None or void.")
  96         #
  97         _X = numpy.asmatrix(numpy.ravel( X )).T
  98         #
  99         if self.__dX is None:
 100             _dX  = self.__increment * _X
 101         else:
 102             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
 103         #
 104         if (_dX == 0.).any():
 105             moyenne = _dX.mean()
 106             if moyenne == 0.:
 107                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 108             else:
 109                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 110         #
 111         if self.__centeredDF:
 112             #
 113             _Jacobienne  = []
 114             for i in range( len(_dX) ):
 115                 _dXi            = _dX[i]
 116                 _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 117                 _X_plus_dXi[i]  = _X[i] + _dXi
 118                 _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 119                 _X_moins_dXi[i] = _X[i] - _dXi
 120                 #
 121                 _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 122                 _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 123                 #
 124                 _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
 125             #
 126         else:
 127             #
 128             _Jacobienne  = []
 129             _HX = self.DirectOperator( _X )
 130             for i in range( len(_dX) ):
 131                 _dXi            = _dX[i]
 132                 _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 133                 _X_plus_dXi[i]  = _X[i] + _dXi
 134                 #
 135                 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 136                 #
 137                 _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 138             #
 139         #
 140         _Jacobienne = numpy.matrix( numpy.vstack( _Jacobienne ) ).T
 141         logging.debug("FDA Fin du calcul de la Jacobienne")
 142         #
 143         return _Jacobienne
 144
 145     # ---------------------------------------------------------
 146     def TangentOperator(self, (X, dX) ):
 147         """
 148         Calcul du tangent à l'aide de la Jacobienne.
 149         """
 150         _Jacobienne = self.TangentMatrix( X )
 151         if dX is None or len(dX) == 0:
 152             #
 153             # Calcul de la forme matricielle si le second argument est None
 154             # -------------------------------------------------------------
 155             return _Jacobienne
 156         else:
 157             #
 158             # Calcul de la valeur linéarisée de H en X appliqué à dX
 159             # ------------------------------------------------------
 160             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 161             _HtX = numpy.dot(_Jacobienne, _dX)
 162             return _HtX.A1
 163
 164     # ---------------------------------------------------------
 165     def AdjointOperator(self, (X, Y) ):
 166         """
 167         Calcul de l'adjoint à l'aide de la Jacobienne.
 168         """
 169         _JacobienneT = self.TangentMatrix( X ).T
 170         if Y is None or len(Y) == 0:
 171             #
 172             # Calcul de la forme matricielle si le second argument est None
 173             # -------------------------------------------------------------
 174             return _JacobienneT
 175         else:
 176             #
 177             # Calcul de la valeur de l'adjoint en X appliqué à Y
 178             # --------------------------------------------------
 179             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 180             _HaY = numpy.dot(_JacobienneT, _Y)
 181             return _HaY.A1
 182
 183 # ==============================================================================
 184 #
 185 def test1( XX ):
 186     """ Direct non-linear simulation operator """
 187     #
 188     # NEED TO BE COMPLETED
 189     # NEED TO BE COMPLETED
 190     # NEED TO BE COMPLETED
 191     #
 192     # --------------------------------------> # EXAMPLE TO BE REMOVED
 193     # Example of Identity operator            # EXAMPLE TO BE REMOVED
 194     if type(XX) is type(numpy.matrix([])):    # EXAMPLE TO BE REMOVED
 195         HX = XX.A1.tolist()                   # EXAMPLE TO BE REMOVED
 196     elif type(XX) is type(numpy.array([])):   # EXAMPLE TO BE REMOVED
 197         HX = numpy.matrix(XX).A1.tolist()     # EXAMPLE TO BE REMOVED
 198     else:                                     # EXAMPLE TO BE REMOVED
 199         HX = XX                               # EXAMPLE TO BE REMOVED
 200     #                                         # EXAMPLE TO BE REMOVED
 201     HHX = []                                  # EXAMPLE TO BE REMOVED
 202     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 203     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 204     # --------------------------------------> # EXAMPLE TO BE REMOVED
 205     #
 206     return numpy.array( HHX )
 207
 208 # ==============================================================================
 209 if __name__ == "__main__":
 210
 211     print
 212     print "AUTODIAGNOSTIC"
 213     print "=============="
 214
 215     X0 = [1, 2, 3]
 216
 217     FDA = FDApproximation( test1 )
 218     print "H(X)       =",   FDA.DirectOperator( X0 )
 219     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 220     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 221     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
 222     print
 223     del FDA
 224     FDA = FDApproximation( test1, centeredDF=True )
 225     print "H(X)       =",   FDA.DirectOperator( X0 )
 226     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 227     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 228     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
 229     print
 230     del FDA
 231
 232     print "=============="
 233     print
 234     X0 = range(5)
 235
 236     FDA = FDApproximation( test1 )
 237     print "H(X)       =",   FDA.DirectOperator( X0 )
 238     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 239     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 240     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 241     print
 242     del FDA
 243     FDA = FDApproximation( test1, centeredDF=True )
 244     print "H(X)       =",   FDA.DirectOperator( X0 )
 245     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 246     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 247     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 248     print
 249     del FDA
 250
 251     print "=============="
 252     print
 253     X0 = numpy.arange(3)
 254
 255     FDA = FDApproximation( test1 )
 256     print "H(X)       =",   FDA.DirectOperator( X0 )
 257     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 258     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 259     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 260     print
 261     del FDA
 262     FDA = FDApproximation( test1, centeredDF=True )
 263     print "H(X)       =",   FDA.DirectOperator( X0 )
 264     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 265     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 266     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 267     print
 268     del FDA
 269
 270     print "=============="
 271     print
 272     X0 = numpy.asmatrix(numpy.arange(4)).T
 273
 274     FDA = FDApproximation( test1 )
 275     print "H(X)       =",   FDA.DirectOperator( X0 )
 276     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 277     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 278     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 279     print
 280     del FDA
 281     FDA = FDApproximation( test1, centeredDF=True )
 282     print "H(X)       =",   FDA.DirectOperator( X0 )
 283     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 284     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 285     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 286     print
 287     del FDA
 288