src/daComposant/daNumerics/ApproximatedDerivatives.py

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