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, copy
  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-18,
  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.__lenghtRH = int(lenghtOfRedundancy)
  58             self.__lenghtRJ = int(lenghtOfRedundancy)
  59             self.__listDPCP = [] # Direct Operator Previous Calculated Points
  60             self.__listDPCR = [] # Direct Operator Previous Calculated Results
  61             self.__listDPCN = [] # Direct Operator Previous Calculated Point Norms
  62             self.__listJPCP = [] # Jacobian Previous Calculated Points
  63             self.__listJPCI = [] # Jacobian Previous Calculated Increment
  64             self.__listJPCR = [] # Jacobian Previous Calculated Results
  65             self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
  66             self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
  67         else:
  68             self.__avoidRC = False
  69         if float(increment) <> 0.:
  70             self.__increment  = float(increment)
  71         else:
  72             self.__increment  = 0.01
  73         if dX is None:
  74             self.__dX     = None
  75         else:
  76             self.__dX     = numpy.asmatrix(numpy.ravel( dX )).T
  77         logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
  78         if self.__avoidRC:
  79             logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
  80
  81     # ---------------------------------------------------------
  82     def __doublon__(self, e, l, n, v=""):
  83         __ac, __iac = False, -1
  84         for i in xrange(len(l)-1,-1,-1):
  85             if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
  86                 __ac, __iac = True, i
  87                 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
  88                 break
  89         return __ac, __iac
  90
  91     # ---------------------------------------------------------
  92     def DirectOperator(self, X ):
  93         """
  94         Calcul du direct à l'aide de la fonction fournie.
  95         """
  96         _X = numpy.asmatrix(numpy.ravel( X )).T
  97         #
  98         __alreadyCalculated = False
  99         if self.__avoidRC:
 100             __alreadyCalculated, __i = self.__doublon__(_X, self.__listDPCP, self.__listDPCN, " H")
 101         #
 102         if __alreadyCalculated:
 103             logging.debug("FDA Calcul DirectOperator (par récupération du doublon %i)"%__i)
 104             _HX = self.__listDPCR[__i]
 105         else:
 106             logging.debug("FDA Calcul DirectOperator (explicite)")
 107             _HX = numpy.ravel(self.__userFunction( _X ))
 108             if self.__avoidRC:
 109                 if self.__lenghtRH < 0: self.__lenghtRH = 2 * _X.size
 110                 if len(self.__listDPCP) > self.__lenghtRH:
 111                     logging.debug("FDA Réduction de la liste de H à %i éléments"%self.__lenghtRH)
 112                     self.__listDPCP.pop(0)
 113                     self.__listDPCR.pop(0)
 114                     self.__listDPCN.pop(0)
 115                 self.__listDPCP.append( copy.copy(_X) )
 116                 self.__listDPCR.append( copy.copy(_HX) )
 117                 self.__listDPCN.append( numpy.linalg.norm(_X) )
 118         #
 119         return _HX
 120
 121     # ---------------------------------------------------------
 122     def TangentMatrix(self, X ):
 123         """
 124         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
 125         c'est-à-dire le gradient de H en X. On utilise des différences finies
 126         directionnelles autour du point X. X est un numpy.matrix.
 127
 128         Différences finies centrées (approximation d'ordre 2):
 129         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
 130            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
 131            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
 132            H( X_moins_dXi )
 133         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
 134            le pas 2*dXi
 135         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 136
 137         Différences finies non centrées (approximation d'ordre 1):
 138         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
 139            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
 140            HX_plus_dXi = H( X_plus_dXi )
 141         2/ On calcule la valeur centrale HX = H(X)
 142         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
 143            le pas dXi
 144         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 145
 146         """
 147         logging.debug("FDA Calcul de la Jacobienne")
 148         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
 149         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
 150         #
 151         if X is None or len(X)==0:
 152             raise ValueError("Nominal point X for approximate derivatives can not be None or void.")
 153         #
 154         _X = numpy.asmatrix(numpy.ravel( X )).T
 155         #
 156         if self.__dX is None:
 157             _dX  = self.__increment * _X
 158         else:
 159             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
 160         #
 161         if (_dX == 0.).any():
 162             moyenne = _dX.mean()
 163             if moyenne == 0.:
 164                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 165             else:
 166                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 167         #
 168         __alreadyCalculated  = False
 169         if self.__avoidRC:
 170             __bidon, __alreadyCalculatedP = self.__doublon__(_X,  self.__listJPCP, self.__listJPPN, None)
 171             __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
 172             if __alreadyCalculatedP == __alreadyCalculatedI > -1:
 173                 __alreadyCalculated, __i = True, __alreadyCalculatedP
 174                 logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i)
 175         #
 176         if __alreadyCalculated:
 177             logging.debug("FDA   Calcul Jacobienne (par récupération du doublon %i)"%__i)
 178             _Jacobienne = self.__listJPCR[__i]
 179         else:
 180             logging.debug("FDA   Calcul Jacobienne (explicite)")
 181             if self.__centeredDF:
 182                 #
 183                 _Jacobienne  = []
 184                 for i in range( _dX.size ):
 185                     _dXi            = _dX[i]
 186                     _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 187                     _X_plus_dXi[i]  = _X[i] + _dXi
 188                     _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 189                     _X_moins_dXi[i] = _X[i] - _dXi
 190                     #
 191                     _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 192                     _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 193                     #
 194                     _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
 195                 #
 196             else:
 197                 #
 198                 _Jacobienne  = []
 199                 _HX = self.DirectOperator( _X )
 200                 for i in range( _dX.size ):
 201                     _dXi            = _dX[i]
 202                     _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 203                     _X_plus_dXi[i]  = _X[i] + _dXi
 204                     #
 205                     _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 206                     #
 207                     _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 208                 #
 209             #
 210             _Jacobienne = numpy.matrix( numpy.vstack( _Jacobienne ) ).T
 211             if self.__avoidRC:
 212                 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
 213                 if len(self.__listJPCP) > self.__lenghtRJ:
 214                     logging.debug("FDA Réduction de la liste de J à %i éléments"%self.__lenghtRJ)
 215                     self.__listJPCP.pop(0)
 216                     self.__listJPCI.pop(0)
 217                     self.__listJPCR.pop(0)
 218                     self.__listJPPN.pop(0)
 219                     self.__listJPIN.pop(0)
 220                 self.__listJPCP.append( copy.copy(_X) )
 221                 self.__listJPCI.append( copy.copy(_dX) )
 222                 self.__listJPCR.append( copy.copy(_Jacobienne) )
 223                 self.__listJPPN.append( numpy.linalg.norm(_X) )
 224                 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
 225         #
 226         logging.debug("FDA Fin du calcul de la Jacobienne")
 227         #
 228         return _Jacobienne
 229
 230     # ---------------------------------------------------------
 231     def TangentOperator(self, (X, dX) ):
 232         """
 233         Calcul du tangent à l'aide de la Jacobienne.
 234         """
 235         _Jacobienne = self.TangentMatrix( X )
 236         if dX is None or len(dX) == 0:
 237             #
 238             # Calcul de la forme matricielle si le second argument est None
 239             # -------------------------------------------------------------
 240             return _Jacobienne
 241         else:
 242             #
 243             # Calcul de la valeur linéarisée de H en X appliqué à dX
 244             # ------------------------------------------------------
 245             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 246             _HtX = numpy.dot(_Jacobienne, _dX)
 247             return _HtX.A1
 248
 249     # ---------------------------------------------------------
 250     def AdjointOperator(self, (X, Y) ):
 251         """
 252         Calcul de l'adjoint à l'aide de la Jacobienne.
 253         """
 254         _JacobienneT = self.TangentMatrix( X ).T
 255         if Y is None or len(Y) == 0:
 256             #
 257             # Calcul de la forme matricielle si le second argument est None
 258             # -------------------------------------------------------------
 259             return _JacobienneT
 260         else:
 261             #
 262             # Calcul de la valeur de l'adjoint en X appliqué à Y
 263             # --------------------------------------------------
 264             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 265             _HaY = numpy.dot(_JacobienneT, _Y)
 266             return _HaY.A1
 267
 268 # ==============================================================================
 269 if __name__ == "__main__":
 270     print '\n AUTODIAGNOSTIC \n'