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.DirectOperator = 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 TangentMatrix(self, X ):
  57         """
  58         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
  59         c'est-à-dire le gradient de H en X. On utilise des différences finies
  60         directionnelles autour du point X. X est un numpy.matrix.
  61
  62         Différences finies centrées :
  63         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
  64            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
  65            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
  66            H( X_moins_dXi )
  67         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
  68            le pas 2*dXi
  69         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
  70
  71         Différences finies non centrées :
  72         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
  73            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
  74            HX_plus_dXi = H( X_plus_dXi )
  75         2/ On calcule la valeur centrale HX = H(X)
  76         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
  77            le pas dXi
  78         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
  79
  80         """
  81         logging.debug("  == Calcul de la Jacobienne")
  82         logging.debug("     Incrément de............: %s*X"%float(self.__increment))
  83         logging.debug("     Approximation centrée...: %s"%(self.__centeredDF))
  84         #
  85         _X = numpy.asmatrix(numpy.ravel( X )).T
  86         #
  87         if self.__dX is None:
  88             _dX  = self.__increment * _X
  89         else:
  90             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
  91         #
  92         if (_dX == 0.).any():
  93             moyenne = _dX.mean()
  94             if moyenne == 0.:
  95                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
  96             else:
  97                 _dX = numpy.where( _dX == 0., moyenne, _dX )
  98         #
  99         if self.__centeredDF:
 100             #
 101             # Boucle de calcul des colonnes de la Jacobienne
 102             # ----------------------------------------------
 103             Jacobienne  = []
 104             for i in range( len(_dX) ):
 105                 X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 106                 X_plus_dXi[i]  = _X[i] + _dX[i]
 107                 X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 108                 X_moins_dXi[i] = _X[i] - _dX[i]
 109                 #
 110                 HX_plus_dXi  = self.DirectOperator( X_plus_dXi )
 111                 HX_moins_dXi = self.DirectOperator( X_moins_dXi )
 112                 #
 113                 HX_Diff = ( HX_plus_dXi - HX_moins_dXi ) / (2.*_dX[i])
 114                 #
 115                 Jacobienne.append( HX_Diff )
 116             #
 117         else:
 118             #
 119             # Boucle de calcul des colonnes de la Jacobienne
 120             # ----------------------------------------------
 121             HX_plus_dX = []
 122             for i in range( len(_dX) ):
 123                 X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 124                 X_plus_dXi[i] = _X[i] + _dX[i]
 125                 #
 126                 HX_plus_dXi = self.DirectOperator( X_plus_dXi )
 127                 #
 128                 HX_plus_dX.append( HX_plus_dXi )
 129             #
 130             # Calcul de la valeur centrale
 131             # ----------------------------
 132             HX = self.DirectOperator( _X )
 133             #
 134             # Calcul effectif de la Jacobienne par différences finies
 135             # -------------------------------------------------------
 136             Jacobienne = []
 137             for i in range( len(_dX) ):
 138                 Jacobienne.append( numpy.ravel(( HX_plus_dX[i] - HX ) / _dX[i]) )
 139         #
 140         Jacobienne = numpy.matrix( numpy.vstack( Jacobienne ) ).T
 141         logging.debug("  == 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
 225     FDA = FDApproximation( test1, centeredDF=True )
 226     print "H(X)       =",   FDA.DirectOperator( X0 )
 227     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 228     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 229     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
 230     print
 231     del FDA
 232
 233     X0 = range(5)
 234
 235     FDA = FDApproximation( test1 )
 236     print "H(X)       =",   FDA.DirectOperator( X0 )
 237     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 238     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 239     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 240     print
 241     del FDA