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("  == 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         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             # Boucle de calcul des colonnes de la Jacobienne
 114             # ----------------------------------------------
 115             _Jacobienne  = []
 116             for i in range( len(_dX) ):
 117                 _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 118                 _X_plus_dXi[i]  = _X[i] + _dX[i]
 119                 _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 120                 _X_moins_dXi[i] = _X[i] - _dX[i]
 121                 #
 122                 _HX_plus_dXi  = self.DirectOperator( _X_plus_dXi )
 123                 _HX_moins_dXi = self.DirectOperator( _X_moins_dXi )
 124                 #
 125                 _HX_Diff = numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dX[i])
 126                 #
 127                 _Jacobienne.append( _HX_Diff )
 128             #
 129         else:
 130             #
 131             # Boucle de calcul des colonnes de la Jacobienne
 132             # ----------------------------------------------
 133             _HX_plus_dX = []
 134             for i in range( len(_dX) ):
 135                 _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 136                 _X_plus_dXi[i] = _X[i] + _dX[i]
 137                 #
 138                 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 139                 #
 140                 _HX_plus_dX.append( _HX_plus_dXi )
 141             #
 142             # Calcul de la valeur centrale
 143             # ----------------------------
 144             _HX = self.DirectOperator( _X )
 145             #
 146             # Calcul effectif de la Jacobienne par différences finies
 147             # -------------------------------------------------------
 148             _Jacobienne = []
 149             for i in range( len(_dX) ):
 150                 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 151         #
 152         _Jacobienne = numpy.matrix( numpy.vstack( _Jacobienne ) ).T
 153         logging.debug("  == Fin du calcul de la Jacobienne")
 154         #
 155         return _Jacobienne
 156
 157     # ---------------------------------------------------------
 158     def TangentOperator(self, (X, dX) ):
 159         """
 160         Calcul du tangent à l'aide de la Jacobienne.
 161         """
 162         _Jacobienne = self.TangentMatrix( X )
 163         if dX is None or len(dX) == 0:
 164             #
 165             # Calcul de la forme matricielle si le second argument est None
 166             # -------------------------------------------------------------
 167             return _Jacobienne
 168         else:
 169             #
 170             # Calcul de la valeur linéarisée de H en X appliqué à dX
 171             # ------------------------------------------------------
 172             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 173             _HtX = numpy.dot(_Jacobienne, _dX)
 174             return _HtX.A1
 175
 176     # ---------------------------------------------------------
 177     def AdjointOperator(self, (X, Y) ):
 178         """
 179         Calcul de l'adjoint à l'aide de la Jacobienne.
 180         """
 181         _JacobienneT = self.TangentMatrix( X ).T
 182         if Y is None or len(Y) == 0:
 183             #
 184             # Calcul de la forme matricielle si le second argument est None
 185             # -------------------------------------------------------------
 186             return _JacobienneT
 187         else:
 188             #
 189             # Calcul de la valeur de l'adjoint en X appliqué à Y
 190             # --------------------------------------------------
 191             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 192             _HaY = numpy.dot(_JacobienneT, _Y)
 193             return _HaY.A1
 194
 195 # ==============================================================================
 196 #
 197 def test1( XX ):
 198     """ Direct non-linear simulation operator """
 199     #
 200     # NEED TO BE COMPLETED
 201     # NEED TO BE COMPLETED
 202     # NEED TO BE COMPLETED
 203     #
 204     # --------------------------------------> # EXAMPLE TO BE REMOVED
 205     # Example of Identity operator            # EXAMPLE TO BE REMOVED
 206     if type(XX) is type(numpy.matrix([])):    # EXAMPLE TO BE REMOVED
 207         HX = XX.A1.tolist()                   # EXAMPLE TO BE REMOVED
 208     elif type(XX) is type(numpy.array([])):   # EXAMPLE TO BE REMOVED
 209         HX = numpy.matrix(XX).A1.tolist()     # EXAMPLE TO BE REMOVED
 210     else:                                     # EXAMPLE TO BE REMOVED
 211         HX = XX                               # EXAMPLE TO BE REMOVED
 212     #                                         # EXAMPLE TO BE REMOVED
 213     HHX = []                                  # EXAMPLE TO BE REMOVED
 214     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 215     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 216     # --------------------------------------> # EXAMPLE TO BE REMOVED
 217     #
 218     return numpy.array( HHX )
 219
 220 # ==============================================================================
 221 if __name__ == "__main__":
 222
 223     print
 224     print "AUTODIAGNOSTIC"
 225     print "=============="
 226
 227     X0 = [1, 2, 3]
 228
 229     FDA = FDApproximation( test1 )
 230     print "H(X)       =",   FDA.DirectOperator( X0 )
 231     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 232     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 233     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
 234     print
 235     del FDA
 236     FDA = FDApproximation( test1, centeredDF=True )
 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(3,3+2*len(X0))) )
 241     print
 242     del FDA
 243
 244     print "=============="
 245     print
 246     X0 = range(5)
 247
 248     FDA = FDApproximation( test1 )
 249     print "H(X)       =",   FDA.DirectOperator( X0 )
 250     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 251     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 252     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 253     print
 254     del FDA
 255     FDA = FDApproximation( test1, centeredDF=True )
 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
 263     print "=============="
 264     print
 265     X0 = numpy.arange(3)
 266
 267     FDA = FDApproximation( test1 )
 268     print "H(X)       =",   FDA.DirectOperator( X0 )
 269     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 270     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 271     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 272     print
 273     del FDA
 274     FDA = FDApproximation( test1, centeredDF=True )
 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
 282     print "=============="
 283     print
 284     X0 = numpy.asmatrix(numpy.arange(4)).T
 285
 286     FDA = FDApproximation( test1 )
 287     print "H(X)       =",   FDA.DirectOperator( X0 )
 288     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 289     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 290     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 291     print
 292     del FDA
 293     FDA = FDApproximation( test1, centeredDF=True )
 294     print "H(X)       =",   FDA.DirectOperator( X0 )
 295     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 296     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 297     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 298     print
 299     del FDA
 300