1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2014 EDF R&D
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.
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.
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
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 Définit les versions approximées des opérateurs tangents et adjoints.
26 __author__ = "Jean-Philippe ARGAUD"
28 import os, numpy, time, copy, types, sys
30 from daCore.BasicObjects import Operator
31 # logging.getLogger().setLevel(logging.DEBUG)
33 # ==============================================================================
34 class FDApproximation:
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.
49 avoidingRedundancy = True,
50 toleranceInRedundancy = 1.e-18,
51 lenghtOfRedundancy = -1,
53 self.__userOperator = Operator( fromMethod = Function )
54 self.__userFunction = self.__userOperator.appliedTo
55 self.__centeredDF = bool(centeredDF)
56 if avoidingRedundancy:
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
70 self.__avoidRC = False
71 if float(increment) <> 0.:
72 self.__increment = float(increment)
74 self.__increment = 0.01
78 self.__dX = numpy.asmatrix(numpy.ravel( dX )).T
79 logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
81 logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
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]:
89 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
93 # ---------------------------------------------------------
94 def DirectOperator(self, X ):
96 Calcul du direct à l'aide de la fonction fournie.
98 _X = numpy.asmatrix(numpy.ravel( X )).T
100 __alreadyCalculated = False
102 __alreadyCalculated, __i = self.__doublon__(_X, self.__listDPCP, self.__listDPCN, " H")
104 if __alreadyCalculated:
105 logging.debug("FDA Calcul DirectOperator (par récupération du doublon %i)"%__i)
106 _HX = self.__listDPCR[__i]
108 logging.debug("FDA Calcul DirectOperator (explicite)")
109 _HX = numpy.ravel(self.__userFunction( _X ))
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) )
123 # ---------------------------------------------------------
124 def TangentMatrix(self, X ):
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.
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 =
135 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
137 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
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
146 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
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))
153 if X is None or len(X)==0:
154 raise ValueError("Nominal point X for approximate derivatives can not be None or void.")
156 _X = numpy.asmatrix(numpy.ravel( X )).T
158 if self.__dX is None:
159 _dX = self.__increment * _X
161 _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
163 if (_dX == 0.).any():
166 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
168 _dX = numpy.where( _dX == 0., moyenne, _dX )
170 __alreadyCalculated = False
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)
178 if __alreadyCalculated:
179 logging.debug("FDA Calcul Jacobienne (par récupération du doublon %i)"%__i)
180 _Jacobienne = self.__listJPCR[__i]
182 logging.debug("FDA Calcul Jacobienne (explicite)")
183 if self.__centeredDF:
186 for i in range( _dX.size ):
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
193 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
194 _HX_moins_dXi = self.DirectOperator( _X_moins_dXi )
196 _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
201 _HX = self.DirectOperator( _X )
202 for i in range( _dX.size ):
204 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
205 _X_plus_dXi[i] = _X[i] + _dXi
207 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
209 _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
212 _Jacobienne = numpy.matrix( numpy.vstack( _Jacobienne ) ).T
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) )
228 logging.debug("FDA Fin du calcul de la Jacobienne")
232 # ---------------------------------------------------------
233 def TangentOperator(self, (X, dX) ):
235 Calcul du tangent à l'aide de la Jacobienne.
237 _Jacobienne = self.TangentMatrix( X )
238 if dX is None or len(dX) == 0:
240 # Calcul de la forme matricielle si le second argument est None
241 # -------------------------------------------------------------
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)
251 # ---------------------------------------------------------
252 def AdjointOperator(self, (X, Y) ):
254 Calcul de l'adjoint à l'aide de la Jacobienne.
256 _JacobienneT = self.TangentMatrix( X ).T
257 if Y is None or len(Y) == 0:
259 # Calcul de la forme matricielle si le second argument est None
260 # -------------------------------------------------------------
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)
270 # ==============================================================================
271 if __name__ == "__main__":
272 print '\n AUTODIAGNOSTIC \n'