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 "FunctionH", on obtient un
  37     objet qui dispose de 3 méthodes "FunctionH", "TangentH" et "AdjointH". On
  38     contrôle l'approximation DF avec l'incrément multiplicatif "increment"
  39     valant par défaut 1%, ou avec l'incrément fixe "dX" qui sera multiplié par
  40     "increment" (donc en %), et on effectue de DF centrées si le booléen
  41     "centeredDF" est vrai.
  42     """
  43     def __init__(self, FunctionH = None, centeredDF = False, increment = 0.01, dX = None):
  44         self.FunctionH    = FunctionH
  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(dX).flatten().T
  54
  55     # ---------------------------------------------------------
  56     def TangentHMatrix(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(X).flatten().T
  86         #
  87         if self.__dX is None:
  88             _dX  = self.__increment * _X
  89         else:
  90             _dX = numpy.asmatrix(self.__dX).flatten().T
  91         logging.debug("     Incrément non corrigé...: %s"%(_dX))
  92         #
  93         if (_dX == 0.).any():
  94             moyenne = _dX.mean()
  95             if moyenne == 0.:
  96                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
  97             else:
  98                 _dX = numpy.where( _dX == 0., moyenne, _dX )
  99         logging.debug("     Incrément dX............: %s"%(_dX))
 100         #
 101         if self.__centeredDF:
 102             #
 103             # Boucle de calcul des colonnes de la Jacobienne
 104             # ----------------------------------------------
 105             Jacobienne  = []
 106             for i in range( len(_dX) ):
 107                 X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 108                 X_plus_dXi[i]  = _X[i] + _dX[i]
 109                 X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 110                 X_moins_dXi[i] = _X[i] - _dX[i]
 111                 #
 112                 HX_plus_dXi  = self.FunctionH( X_plus_dXi )
 113                 HX_moins_dXi = self.FunctionH( X_moins_dXi )
 114                 #
 115                 HX_Diff = ( HX_plus_dXi - HX_moins_dXi ) / (2.*_dX[i])
 116                 #
 117                 Jacobienne.append( HX_Diff )
 118             #
 119         else:
 120             #
 121             # Boucle de calcul des colonnes de la Jacobienne
 122             # ----------------------------------------------
 123             HX_plus_dX = []
 124             for i in range( len(_dX) ):
 125                 X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 126                 X_plus_dXi[i] = _X[i] + _dX[i]
 127                 #
 128                 HX_plus_dXi = self.FunctionH( X_plus_dXi )
 129                 #
 130                 HX_plus_dX.append( HX_plus_dXi )
 131             #
 132             # Calcul de la valeur centrale
 133             # ----------------------------
 134             HX = self.FunctionH( _X )
 135             #
 136             # Calcul effectif de la Jacobienne par différences finies
 137             # -------------------------------------------------------
 138             Jacobienne = []
 139             for i in range( len(_dX) ):
 140                 Jacobienne.append( numpy.asmatrix(( HX_plus_dX[i] - HX ) / _dX[i]).flatten().A1 )
 141         #
 142         Jacobienne = numpy.matrix( numpy.vstack( Jacobienne ) ).T
 143         logging.debug("  == Fin du calcul de la Jacobienne")
 144         #
 145         return Jacobienne
 146
 147     # ---------------------------------------------------------
 148     def TangentH(self, (X, dX) ):
 149         """
 150         Calcul du tangent à l'aide de la Jacobienne.
 151         """
 152         Jacobienne = self.TangentHMatrix( X )
 153         if dX is None or len(dX) == 0:
 154             #
 155             # Calcul de la forme matricielle si le second argument est None
 156             # -------------------------------------------------------------
 157             return Jacobienne
 158         else:
 159             #
 160             # Calcul de la valeur linéarisée de H en X appliqué à dX
 161             # ------------------------------------------------------
 162             dX = numpy.asmatrix(dX).flatten().T
 163             HtX = numpy.dot(Jacobienne, dX)
 164             return HtX.A1
 165
 166     # ---------------------------------------------------------
 167     def AdjointH(self, (X, Y) ):
 168         """
 169         Calcul de l'adjoint à l'aide de la Jacobienne.
 170         """
 171         JacobienneT = self.TangentHMatrix( X ).T
 172         if Y is None or len(Y) == 0:
 173             #
 174             # Calcul de la forme matricielle si le second argument est None
 175             # -------------------------------------------------------------
 176             return JacobienneT
 177         else:
 178             #
 179             # Calcul de la valeur de l'adjoint en X appliqué à Y
 180             # --------------------------------------------------
 181             Y = numpy.asmatrix(Y).flatten().T
 182             HaY = numpy.dot(JacobienneT, Y)
 183             return HaY.A1
 184
 185 # ==============================================================================
 186 #
 187 def test1FunctionH( XX ):
 188     """ Direct non-linear simulation operator """
 189     #
 190     # NEED TO BE COMPLETED
 191     # NEED TO BE COMPLETED
 192     # NEED TO BE COMPLETED
 193     #
 194     # --------------------------------------> # EXAMPLE TO BE REMOVED
 195     # Example of Identity operator            # EXAMPLE TO BE REMOVED
 196     if type(XX) is type(numpy.matrix([])):    # EXAMPLE TO BE REMOVED
 197         HX = XX.A1.tolist()                   # EXAMPLE TO BE REMOVED
 198     elif type(XX) is type(numpy.array([])):   # EXAMPLE TO BE REMOVED
 199         HX = numpy.matrix(XX).A1.tolist()     # EXAMPLE TO BE REMOVED
 200     else:                                     # EXAMPLE TO BE REMOVED
 201         HX = XX                               # EXAMPLE TO BE REMOVED
 202     #                                         # EXAMPLE TO BE REMOVED
 203     HHX = []                                  # EXAMPLE TO BE REMOVED
 204     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 205     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 206     # --------------------------------------> # EXAMPLE TO BE REMOVED
 207     #
 208     return numpy.array( HHX )
 209
 210 # ==============================================================================
 211 if __name__ == "__main__":
 212
 213     print
 214     print "AUTODIAGNOSTIC"
 215     print "=============="
 216
 217     X0 = [1, 2, 3]
 218
 219     FDA = FDApproximation( test1FunctionH )
 220     print "H(X)       =",   FDA.FunctionH( X0 )
 221     print "Tg matrice =\n", FDA.TangentHMatrix( X0 )
 222     print "Tg(X)      =",   FDA.TangentH( (X0, X0) )
 223     print "Ad((X,Y))  =",   FDA.AdjointH( (X0,range(3,3+2*len(X0))) )
 224     print
 225     del FDA
 226
 227     FDA = FDApproximation( test1FunctionH, centeredDF=True )
 228     print "H(X)       =",   FDA.FunctionH( X0 )
 229     print "Tg matrice =\n", FDA.TangentHMatrix( X0 )
 230     print "Tg(X)      =",   FDA.TangentH( (X0, X0) )
 231     print "Ad((X,Y))  =",   FDA.AdjointH( (X0,range(3,3+2*len(X0))) )
 232     print
 233     del FDA
 234
 235     X0 = range(5)
 236
 237     FDA = FDApproximation( test1FunctionH )
 238     print "H(X)       =",   FDA.FunctionH( X0 )
 239     print "Tg matrice =\n", FDA.TangentHMatrix( X0 )
 240     print "Tg(X)      =",   FDA.TangentH( (X0, X0) )
 241     print "Ad((X,Y))  =",   FDA.AdjointH( (X0,range(7,7+2*len(X0))) )
 242     print
 243     del FDA