src/daComposant/daNumerics/ApproximatedDerivatives.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2020 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, copy, types, sys
  29 import logging
  30 from daCore.BasicObjects import Operator
  31 # logging.getLogger().setLevel(logging.DEBUG)
  32
  33 # ==============================================================================
  34 def ExecuteFunction( paire ):
  35     assert len(paire) == 2, "Incorrect number of arguments"
  36     X, funcrepr = paire
  37     __X = numpy.asmatrix(numpy.ravel( X )).T
  38     __sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"])
  39     __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
  40     __fonction = getattr(__module,funcrepr["__userFunction__name"])
  41     sys.path = __sys_path_tmp ; del __sys_path_tmp
  42     __HX  = __fonction( __X )
  43     return numpy.ravel( __HX )
  44
  45 # ==============================================================================
  46 class FDApproximation(object):
  47     """
  48     Cette classe sert d'interface pour définir les opérateurs approximés. A la
  49     création d'un objet, en fournissant une fonction "Function", on obtient un
  50     objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
  51     "AdjointOperator". On contrôle l'approximation DF avec l'incrément
  52     multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
  53     "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
  54     centrées si le booléen "centeredDF" est vrai.
  55     """
  56     def __init__(self,
  57             Function              = None,
  58             centeredDF            = False,
  59             increment             = 0.01,
  60             dX                    = None,
  61             avoidingRedundancy    = True,
  62             toleranceInRedundancy = 1.e-18,
  63             lenghtOfRedundancy    = -1,
  64             mpEnabled             = False,
  65             mpWorkers             = None,
  66             mfEnabled             = False,
  67             ):
  68         if mpEnabled:
  69             try:
  70                 import multiprocessing
  71                 self.__mpEnabled = True
  72             except ImportError:
  73                 self.__mpEnabled = False
  74         else:
  75             self.__mpEnabled = False
  76         self.__mpWorkers = mpWorkers
  77         if self.__mpWorkers is not None and self.__mpWorkers < 1:
  78             self.__mpWorkers = None
  79         logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
  80         #
  81         if mfEnabled:
  82             self.__mfEnabled = True
  83         else:
  84             self.__mfEnabled = False
  85         logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
  86         #
  87         if avoidingRedundancy:
  88             self.__avoidRC = True
  89             self.__tolerBP = float(toleranceInRedundancy)
  90             self.__lenghtRJ = int(lenghtOfRedundancy)
  91             self.__listJPCP = [] # Jacobian Previous Calculated Points
  92             self.__listJPCI = [] # Jacobian Previous Calculated Increment
  93             self.__listJPCR = [] # Jacobian Previous Calculated Results
  94             self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
  95             self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
  96         else:
  97             self.__avoidRC = False
  98         #
  99         if self.__mpEnabled:
 100             if isinstance(Function,types.FunctionType):
 101                 logging.debug("FDA Calculs en multiprocessing : FunctionType")
 102                 self.__userFunction__name = Function.__name__
 103                 try:
 104                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 105                 except:
 106                     mod = os.path.abspath(Function.__globals__['__file__'])
 107                 if not os.path.isfile(mod):
 108                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 109                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 110                 self.__userFunction__path = os.path.dirname(mod)
 111                 del mod
 112                 self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 113                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 114             elif isinstance(Function,types.MethodType):
 115                 logging.debug("FDA Calculs en multiprocessing : MethodType")
 116                 self.__userFunction__name = Function.__name__
 117                 try:
 118                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 119                 except:
 120                     mod = os.path.abspath(Function.__func__.__globals__['__file__'])
 121                 if not os.path.isfile(mod):
 122                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 123                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 124                 self.__userFunction__path = os.path.dirname(mod)
 125                 del mod
 126                 self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 127                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 128             else:
 129                 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
 130         else:
 131             self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 132             self.__userFunction = self.__userOperator.appliedTo
 133         #
 134         self.__centeredDF = bool(centeredDF)
 135         if abs(float(increment)) > 1.e-15:
 136             self.__increment  = float(increment)
 137         else:
 138             self.__increment  = 0.01
 139         if dX is None:
 140             self.__dX     = None
 141         else:
 142             self.__dX     = numpy.asmatrix(numpy.ravel( dX )).T
 143         logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
 144         if self.__avoidRC:
 145             logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
 146
 147     # ---------------------------------------------------------
 148     def __doublon__(self, e, l, n, v=None):
 149         __ac, __iac = False, -1
 150         for i in range(len(l)-1,-1,-1):
 151             if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
 152                 __ac, __iac = True, i
 153                 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
 154                 break
 155         return __ac, __iac
 156
 157     # ---------------------------------------------------------
 158     def DirectOperator(self, X ):
 159         """
 160         Calcul du direct à l'aide de la fonction fournie.
 161         """
 162         logging.debug("FDA Calcul DirectOperator (explicite)")
 163         if self.__mfEnabled:
 164             _HX = self.__userFunction( X, argsAsSerie = True )
 165         else:
 166             _X = numpy.asmatrix(numpy.ravel( X )).T
 167             _HX = numpy.ravel(self.__userFunction( _X ))
 168         #
 169         return _HX
 170
 171     # ---------------------------------------------------------
 172     def TangentMatrix(self, X ):
 173         """
 174         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
 175         c'est-à-dire le gradient de H en X. On utilise des différences finies
 176         directionnelles autour du point X. X est un numpy.matrix.
 177
 178         Différences finies centrées (approximation d'ordre 2):
 179         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
 180            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
 181            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
 182            H( X_moins_dXi )
 183         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
 184            le pas 2*dXi
 185         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 186
 187         Différences finies non centrées (approximation d'ordre 1):
 188         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
 189            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
 190            HX_plus_dXi = H( X_plus_dXi )
 191         2/ On calcule la valeur centrale HX = H(X)
 192         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
 193            le pas dXi
 194         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 195
 196         """
 197         logging.debug("FDA Début du calcul de la Jacobienne")
 198         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
 199         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
 200         #
 201         if X is None or len(X)==0:
 202             raise ValueError("Nominal point X for approximate derivatives can not be None or void (X=%s)."%(str(X),))
 203         #
 204         _X = numpy.asmatrix(numpy.ravel( X )).T
 205         #
 206         if self.__dX is None:
 207             _dX  = self.__increment * _X
 208         else:
 209             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
 210         #
 211         if (_dX == 0.).any():
 212             moyenne = _dX.mean()
 213             if moyenne == 0.:
 214                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 215             else:
 216                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 217         #
 218         __alreadyCalculated  = False
 219         if self.__avoidRC:
 220             __bidon, __alreadyCalculatedP = self.__doublon__(_X,  self.__listJPCP, self.__listJPPN, None)
 221             __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
 222             if __alreadyCalculatedP == __alreadyCalculatedI > -1:
 223                 __alreadyCalculated, __i = True, __alreadyCalculatedP
 224                 logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i)
 225         #
 226         if __alreadyCalculated:
 227             logging.debug("FDA   Calcul Jacobienne (par récupération du doublon %i)"%__i)
 228             _Jacobienne = self.__listJPCR[__i]
 229         else:
 230             logging.debug("FDA   Calcul Jacobienne (explicite)")
 231             if self.__centeredDF:
 232                 #
 233                 if self.__mpEnabled and not self.__mfEnabled:
 234                     funcrepr = {
 235                         "__userFunction__path" : self.__userFunction__path,
 236                         "__userFunction__modl" : self.__userFunction__modl,
 237                         "__userFunction__name" : self.__userFunction__name,
 238                     }
 239                     _jobs = []
 240                     for i in range( len(_dX) ):
 241                         _dXi            = _dX[i]
 242                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 243                         _X_plus_dXi[i]  = _X[i] + _dXi
 244                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 245                         _X_moins_dXi[i] = _X[i] - _dXi
 246                         #
 247                         _jobs.append( (_X_plus_dXi,  funcrepr) )
 248                         _jobs.append( (_X_moins_dXi, funcrepr) )
 249                     #
 250                     import multiprocessing
 251                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 252                     _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
 253                     self.__pool.close()
 254                     self.__pool.join()
 255                     #
 256                     _Jacobienne  = []
 257                     for i in range( len(_dX) ):
 258                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 259                     #
 260                 elif self.__mfEnabled:
 261                     _xserie = []
 262                     for i in range( len(_dX) ):
 263                         _dXi            = _dX[i]
 264                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 265                         _X_plus_dXi[i]  = _X[i] + _dXi
 266                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 267                         _X_moins_dXi[i] = _X[i] - _dXi
 268                         #
 269                         _xserie.append( _X_plus_dXi )
 270                         _xserie.append( _X_moins_dXi )
 271                     #
 272                     _HX_plusmoins_dX = self.DirectOperator( _xserie )
 273                      #
 274                     _Jacobienne  = []
 275                     for i in range( len(_dX) ):
 276                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 277                     #
 278                 else:
 279                     _Jacobienne  = []
 280                     for i in range( _dX.size ):
 281                         _dXi            = _dX[i]
 282                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 283                         _X_plus_dXi[i]  = _X[i] + _dXi
 284                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 285                         _X_moins_dXi[i] = _X[i] - _dXi
 286                         #
 287                         _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 288                         _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 289                         #
 290                         _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
 291                 #
 292             else:
 293                 #
 294                 if self.__mpEnabled and not self.__mfEnabled:
 295                     funcrepr = {
 296                         "__userFunction__path" : self.__userFunction__path,
 297                         "__userFunction__modl" : self.__userFunction__modl,
 298                         "__userFunction__name" : self.__userFunction__name,
 299                     }
 300                     _jobs = []
 301                     _jobs.append( (_X.A1, funcrepr) )
 302                     for i in range( len(_dX) ):
 303                         _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 304                         _X_plus_dXi[i] = _X[i] + _dX[i]
 305                         #
 306                         _jobs.append( (_X_plus_dXi, funcrepr) )
 307                     #
 308                     import multiprocessing
 309                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 310                     _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
 311                     self.__pool.close()
 312                     self.__pool.join()
 313                     #
 314                     _HX = _HX_plus_dX.pop(0)
 315                     #
 316                     _Jacobienne = []
 317                     for i in range( len(_dX) ):
 318                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 319                     #
 320                 elif self.__mfEnabled:
 321                     _xserie = []
 322                     _xserie.append( _X.A1 )
 323                     for i in range( len(_dX) ):
 324                         _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 325                         _X_plus_dXi[i] = _X[i] + _dX[i]
 326                         #
 327                         _xserie.append( _X_plus_dXi )
 328                     #
 329                     _HX_plus_dX = self.DirectOperator( _xserie )
 330                     #
 331                     _HX = _HX_plus_dX.pop(0)
 332                     #
 333                     _Jacobienne = []
 334                     for i in range( len(_dX) ):
 335                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 336                    #
 337                 else:
 338                     _Jacobienne  = []
 339                     _HX = self.DirectOperator( _X )
 340                     for i in range( _dX.size ):
 341                         _dXi            = _dX[i]
 342                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 343                         _X_plus_dXi[i]  = _X[i] + _dXi
 344                         #
 345                         _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 346                         #
 347                         _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 348                 #
 349             #
 350             _Jacobienne = numpy.asmatrix( numpy.vstack( _Jacobienne ) ).T
 351             if self.__avoidRC:
 352                 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
 353                 while len(self.__listJPCP) > self.__lenghtRJ:
 354                     self.__listJPCP.pop(0)
 355                     self.__listJPCI.pop(0)
 356                     self.__listJPCR.pop(0)
 357                     self.__listJPPN.pop(0)
 358                     self.__listJPIN.pop(0)
 359                 self.__listJPCP.append( copy.copy(_X) )
 360                 self.__listJPCI.append( copy.copy(_dX) )
 361                 self.__listJPCR.append( copy.copy(_Jacobienne) )
 362                 self.__listJPPN.append( numpy.linalg.norm(_X) )
 363                 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
 364         #
 365         logging.debug("FDA Fin du calcul de la Jacobienne")
 366         #
 367         return _Jacobienne
 368
 369     # ---------------------------------------------------------
 370     def TangentOperator(self, paire ):
 371         """
 372         Calcul du tangent à l'aide de la Jacobienne.
 373         """
 374         if self.__mfEnabled:
 375             assert len(paire) == 1, "Incorrect lenght of arguments"
 376             _paire = paire[0]
 377             assert len(_paire) == 2, "Incorrect number of arguments"
 378         else:
 379             assert len(paire) == 2, "Incorrect number of arguments"
 380             _paire = paire
 381         X, dX = _paire
 382         _Jacobienne = self.TangentMatrix( X )
 383         if dX is None or len(dX) == 0:
 384             #
 385             # Calcul de la forme matricielle si le second argument est None
 386             # -------------------------------------------------------------
 387             if self.__mfEnabled: return [_Jacobienne,]
 388             else:                return _Jacobienne
 389         else:
 390             #
 391             # Calcul de la valeur linéarisée de H en X appliqué à dX
 392             # ------------------------------------------------------
 393             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 394             _HtX = numpy.dot(_Jacobienne, _dX)
 395             if self.__mfEnabled: return [_HtX.A1,]
 396             else:                return _HtX.A1
 397
 398     # ---------------------------------------------------------
 399     def AdjointOperator(self, paire ):
 400         """
 401         Calcul de l'adjoint à l'aide de la Jacobienne.
 402         """
 403         if self.__mfEnabled:
 404             assert len(paire) == 1, "Incorrect lenght of arguments"
 405             _paire = paire[0]
 406             assert len(_paire) == 2, "Incorrect number of arguments"
 407         else:
 408             assert len(paire) == 2, "Incorrect number of arguments"
 409             _paire = paire
 410         X, Y = _paire
 411         _JacobienneT = self.TangentMatrix( X ).T
 412         if Y is None or len(Y) == 0:
 413             #
 414             # Calcul de la forme matricielle si le second argument est None
 415             # -------------------------------------------------------------
 416             if self.__mfEnabled: return [_JacobienneT,]
 417             else:                return _JacobienneT
 418         else:
 419             #
 420             # Calcul de la valeur de l'adjoint en X appliqué à Y
 421             # --------------------------------------------------
 422             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 423             _HaY = numpy.dot(_JacobienneT, _Y)
 424             if self.__mfEnabled: return [_HaY.A1,]
 425             else:                return _HaY.A1
 426
 427 # ==============================================================================
 428 if __name__ == "__main__":
 429     print('\n AUTODIAGNOSTIC\n')