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,
  44             Function              = None,
  45             centeredDF            = False,
  46             increment             = 0.01,
  47             dX                    = None,
  48             avoidingRedundancy    = True,
  49             toleranceInRedundancy = 1.e-15,
  50             lenghtOfRedundancy    = -1,
  51             ):
  52         self.__userFunction = Function
  53         self.__centeredDF = bool(centeredDF)
  54         if avoidingRedundancy:
  55             self.__avoidRC = True
  56             self.__tolerBP = float(toleranceInRedundancy)
  57             self.__lenghtR = int(lenghtOfRedundancy)
  58             self.__listDPCP = [] # Direct Operator Previous Calculated Points
  59             self.__listDPCR = [] # Direct Operator Previous Calculated Results
  60             self.__listDPCN = [] # Direct Operator Previous Calculated Point Norms
  61             self.__listJPCP = [] # Jacobian Previous Calculated Points
  62             self.__listJPCI = [] # Jacobian Previous Calculated Increment
  63             self.__listJPCR = [] # Jacobian Previous Calculated Results
  64             self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
  65             self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
  66         else:
  67             self.__avoidRC = False
  68         if float(increment) <> 0.:
  69             self.__increment  = float(increment)
  70         else:
  71             self.__increment  = 0.01
  72         if dX is None:
  73             self.__dX     = None
  74         else:
  75             self.__dX     = numpy.asmatrix(numpy.ravel( dX )).T
  76         logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
  77         if self.__avoidRC:
  78             logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
  79
  80     # ---------------------------------------------------------
  81     def __doublon__(self, e, l, n, v=""):
  82         __ac, __i = False, -1
  83         for i in xrange(len(l)-1,-1,-1):
  84             if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
  85                 __ac, __i = True, i
  86                 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon"%v)
  87                 break
  88         return __ac, __i
  89
  90     # ---------------------------------------------------------
  91     def DirectOperator(self, X ):
  92         """
  93         Calcul du direct à l'aide de la fonction fournie.
  94         """
  95         _X = numpy.asmatrix(numpy.ravel( X )).T
  96         #
  97         __alreadyCalculated = False
  98         if self.__avoidRC:
  99             __alreadyCalculated, __i = self.__doublon__(_X, self.__listDPCP, self.__listDPCN, " H")
 100         #
 101         if __alreadyCalculated:
 102             logging.debug("FDA Calcul DirectOperator (par récupération)")
 103             _HX = self.__listDPCR[__i]
 104         else:
 105             logging.debug("FDA Calcul DirectOperator (explicite)")
 106             _HX  = self.__userFunction( _X )
 107             if self.__avoidRC:
 108                 if self.__lenghtR < 0: self.__lenghtR = 2 * _X.size
 109                 if len(self.__listDPCP) > self.__lenghtR:
 110                     self.__listDPCP.pop(0)
 111                     self.__listDPCR.pop(0)
 112                     self.__listDPCN.pop(0)
 113                 self.__listDPCP.append( _X )
 114                 self.__listDPCR.append( _HX )
 115                 self.__listDPCN.append( numpy.linalg.norm(_X) )
 116         #
 117         return numpy.ravel( _HX )
 118
 119     # ---------------------------------------------------------
 120     def TangentMatrix(self, X ):
 121         """
 122         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
 123         c'est-à-dire le gradient de H en X. On utilise des différences finies
 124         directionnelles autour du point X. X est un numpy.matrix.
 125
 126         Différences finies centrées (approximation d'ordre 2):
 127         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
 128            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
 129            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
 130            H( X_moins_dXi )
 131         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
 132            le pas 2*dXi
 133         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 134
 135         Différences finies non centrées (approximation d'ordre 1):
 136         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
 137            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
 138            HX_plus_dXi = H( X_plus_dXi )
 139         2/ On calcule la valeur centrale HX = H(X)
 140         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
 141            le pas dXi
 142         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 143
 144         """
 145         logging.debug("FDA Calcul de la Jacobienne")
 146         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
 147         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
 148         #
 149         if X is None or len(X)==0:
 150             raise ValueError("Nominal point X for approximate derivatives can not be None or void.")
 151         #
 152         _X = numpy.asmatrix(numpy.ravel( X )).T
 153         #
 154         if self.__dX is None:
 155             _dX  = self.__increment * _X
 156         else:
 157             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
 158         #
 159         if (_dX == 0.).any():
 160             moyenne = _dX.mean()
 161             if moyenne == 0.:
 162                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 163             else:
 164                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 165         #
 166         __alreadyCalculated  = False
 167         if self.__avoidRC:
 168             __bidon, __alreadyCalculatedP = self.__doublon__(_X,  self.__listJPCP, self.__listJPPN, None)
 169             __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
 170             if __alreadyCalculatedP == __alreadyCalculatedI > -1:
 171                 __alreadyCalculated, __i = True, __alreadyCalculatedP
 172         #
 173         if __alreadyCalculated:
 174             logging.debug("FDA   Calcul Jacobienne (par récupération)")
 175             _Jacobienne = self.__listJPCR[__i]
 176         else:
 177             logging.debug("FDA   Calcul Jacobienne (explicite)")
 178             if self.__centeredDF:
 179                 #
 180                 _Jacobienne  = []
 181                 for i in range( _dX.size ):
 182                     _dXi            = _dX[i]
 183                     _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 184                     _X_plus_dXi[i]  = _X[i] + _dXi
 185                     _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 186                     _X_moins_dXi[i] = _X[i] - _dXi
 187                     #
 188                     _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 189                     _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 190                     #
 191                     _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
 192                 #
 193             else:
 194                 #
 195                 _Jacobienne  = []
 196                 _HX = self.DirectOperator( _X )
 197                 for i in range( _dX.size ):
 198                     _dXi            = _dX[i]
 199                     _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 200                     _X_plus_dXi[i]  = _X[i] + _dXi
 201                     #
 202                     _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 203                     #
 204                     _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 205                 #
 206             #
 207             _Jacobienne = numpy.matrix( numpy.vstack( _Jacobienne ) ).T
 208             if self.__avoidRC:
 209                 if self.__lenghtR < 0: self.__lenghtR = 2 * _X.size
 210                 if len(self.__listJPCP) > self.__lenghtR:
 211                     self.__listJPCP.pop(0)
 212                     self.__listJPCI.pop(0)
 213                     self.__listJPCR.pop(0)
 214                     self.__listJPPN.pop(0)
 215                     self.__listJPIN.pop(0)
 216                 self.__listJPCP.append( _X )
 217                 self.__listJPCI.append( _dX )
 218                 self.__listJPCR.append( _Jacobienne )
 219                 self.__listJPPN.append( numpy.linalg.norm(_X) )
 220                 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
 221         #
 222         logging.debug("FDA Fin du calcul de la Jacobienne")
 223         #
 224         return _Jacobienne
 225
 226     # ---------------------------------------------------------
 227     def TangentOperator(self, (X, dX) ):
 228         """
 229         Calcul du tangent à l'aide de la Jacobienne.
 230         """
 231         _Jacobienne = self.TangentMatrix( X )
 232         if dX is None or len(dX) == 0:
 233             #
 234             # Calcul de la forme matricielle si le second argument est None
 235             # -------------------------------------------------------------
 236             return _Jacobienne
 237         else:
 238             #
 239             # Calcul de la valeur linéarisée de H en X appliqué à dX
 240             # ------------------------------------------------------
 241             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 242             _HtX = numpy.dot(_Jacobienne, _dX)
 243             return _HtX.A1
 244
 245     # ---------------------------------------------------------
 246     def AdjointOperator(self, (X, Y) ):
 247         """
 248         Calcul de l'adjoint à l'aide de la Jacobienne.
 249         """
 250         _JacobienneT = self.TangentMatrix( X ).T
 251         if Y is None or len(Y) == 0:
 252             #
 253             # Calcul de la forme matricielle si le second argument est None
 254             # -------------------------------------------------------------
 255             return _JacobienneT
 256         else:
 257             #
 258             # Calcul de la valeur de l'adjoint en X appliqué à Y
 259             # --------------------------------------------------
 260             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 261             _HaY = numpy.dot(_JacobienneT, _Y)
 262             return _HaY.A1
 263
 264 # ==============================================================================
 265 #
 266 def test1( XX ):
 267     """ Direct non-linear simulation operator """
 268     #
 269     # NEED TO BE COMPLETED
 270     # NEED TO BE COMPLETED
 271     # NEED TO BE COMPLETED
 272     #
 273     # --------------------------------------> # EXAMPLE TO BE REMOVED
 274     # Example of Identity operator            # EXAMPLE TO BE REMOVED
 275     if type(XX) is type(numpy.matrix([])):    # EXAMPLE TO BE REMOVED
 276         HX = XX.A1.tolist()                   # EXAMPLE TO BE REMOVED
 277     elif type(XX) is type(numpy.array([])):   # EXAMPLE TO BE REMOVED
 278         HX = numpy.matrix(XX).A1.tolist()     # EXAMPLE TO BE REMOVED
 279     else:                                     # EXAMPLE TO BE REMOVED
 280         HX = XX                               # EXAMPLE TO BE REMOVED
 281     #                                         # EXAMPLE TO BE REMOVED
 282     HHX = []                                  # EXAMPLE TO BE REMOVED
 283     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 284     HHX.extend( HX )                          # EXAMPLE TO BE REMOVED
 285     # --------------------------------------> # EXAMPLE TO BE REMOVED
 286     #
 287     return numpy.array( HHX )
 288
 289 # ==============================================================================
 290 if __name__ == "__main__":
 291
 292     print
 293     print "AUTODIAGNOSTIC"
 294     print "=============="
 295
 296     X0 = [1, 2, 3]
 297
 298     FDA = FDApproximation( test1 )
 299     print "H(X)       =",   FDA.DirectOperator( X0 )
 300     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 301     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 302     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
 303     print
 304     del FDA
 305     FDA = FDApproximation( test1, centeredDF=True )
 306     print "H(X)       =",   FDA.DirectOperator( X0 )
 307     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 308     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 309     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
 310     print
 311     del FDA
 312
 313     print "=============="
 314     print
 315     X0 = range(5)
 316
 317     FDA = FDApproximation( test1 )
 318     print "H(X)       =",   FDA.DirectOperator( X0 )
 319     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 320     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 321     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 322     print
 323     del FDA
 324     FDA = FDApproximation( test1, centeredDF=True )
 325     print "H(X)       =",   FDA.DirectOperator( X0 )
 326     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 327     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 328     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 329     print
 330     del FDA
 331
 332     print "=============="
 333     print
 334     X0 = numpy.arange(3)
 335
 336     FDA = FDApproximation( test1 )
 337     print "H(X)       =",   FDA.DirectOperator( X0 )
 338     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 339     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 340     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 341     print
 342     del FDA
 343     FDA = FDApproximation( test1, centeredDF=True )
 344     print "H(X)       =",   FDA.DirectOperator( X0 )
 345     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 346     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 347     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 348     print
 349     del FDA
 350
 351     print "=============="
 352     print
 353     X0 = numpy.asmatrix(numpy.arange(4)).T
 354
 355     FDA = FDApproximation( test1 )
 356     print "H(X)       =",   FDA.DirectOperator( X0 )
 357     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 358     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 359     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 360     print
 361     del FDA
 362     FDA = FDApproximation( test1, centeredDF=True )
 363     print "H(X)       =",   FDA.DirectOperator( X0 )
 364     print "Tg matrice =\n", FDA.TangentMatrix( X0 )
 365     print "Tg(X)      =",   FDA.TangentOperator( (X0, X0) )
 366     print "Ad((X,Y))  =",   FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
 367     print
 368     del FDA
 369