1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2013 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
30 # logging.getLogger().setLevel(logging.DEBUG)
32 # ==============================================================================
33 class FDApproximation:
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.
43 def __init__(self, Function = None, centeredDF = False, increment = 0.01, dX = None):
44 self.__userFunction = Function
45 self.__centeredDF = bool(centeredDF)
46 if float(increment) <> 0.:
47 self.__increment = float(increment)
49 self.__increment = 0.01
53 self.__dX = numpy.asmatrix(numpy.ravel( dX )).T
55 # ---------------------------------------------------------
56 def DirectOperator(self, X ):
58 Calcul du direct à l'aide de la fonction fournie.
60 _X = numpy.asmatrix(numpy.ravel( X )).T
61 _HX = self.__userFunction( _X )
62 return numpy.ravel( _HX )
64 # ---------------------------------------------------------
65 def TangentMatrix(self, X ):
67 Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
68 c'est-à-dire le gradient de H en X. On utilise des différences finies
69 directionnelles autour du point X. X est un numpy.matrix.
71 Différences finies centrées :
72 1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
73 dX[i] à la composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
74 on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
76 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
78 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
80 Différences finies non centrées :
81 1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
82 composante X[i] pour composer X_plus_dXi, et on calcule la réponse
83 HX_plus_dXi = H( X_plus_dXi )
84 2/ On calcule la valeur centrale HX = H(X)
85 3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
87 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
90 logging.debug(" == Calcul de la Jacobienne")
91 logging.debug(" Incrément de............: %s*X"%float(self.__increment))
92 logging.debug(" Approximation centrée...: %s"%(self.__centeredDF))
94 if X is None or len(X)==0:
95 raise ValueError("Nominal point X for approximate derivatives can not be None or void.")
97 _X = numpy.asmatrix(numpy.ravel( X )).T
100 _dX = self.__increment * _X
102 _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
104 if (_dX == 0.).any():
107 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
109 _dX = numpy.where( _dX == 0., moyenne, _dX )
111 if self.__centeredDF:
113 # Boucle de calcul des colonnes de la Jacobienne
114 # ----------------------------------------------
116 for i in range( len(_dX) ):
117 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
118 _X_plus_dXi[i] = _X[i] + _dX[i]
119 _X_moins_dXi = numpy.array( _X.A1, dtype=float )
120 _X_moins_dXi[i] = _X[i] - _dX[i]
122 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
123 _HX_moins_dXi = self.DirectOperator( _X_moins_dXi )
125 _HX_Diff = numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dX[i])
127 _Jacobienne.append( _HX_Diff )
131 # Boucle de calcul des colonnes de la Jacobienne
132 # ----------------------------------------------
134 for i in range( len(_dX) ):
135 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
136 _X_plus_dXi[i] = _X[i] + _dX[i]
138 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
140 _HX_plus_dX.append( _HX_plus_dXi )
142 # Calcul de la valeur centrale
143 # ----------------------------
144 _HX = self.DirectOperator( _X )
146 # Calcul effectif de la Jacobienne par différences finies
147 # -------------------------------------------------------
149 for i in range( len(_dX) ):
150 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
152 _Jacobienne = numpy.matrix( numpy.vstack( _Jacobienne ) ).T
153 logging.debug(" == Fin du calcul de la Jacobienne")
157 # ---------------------------------------------------------
158 def TangentOperator(self, (X, dX) ):
160 Calcul du tangent à l'aide de la Jacobienne.
162 _Jacobienne = self.TangentMatrix( X )
163 if dX is None or len(dX) == 0:
165 # Calcul de la forme matricielle si le second argument est None
166 # -------------------------------------------------------------
170 # Calcul de la valeur linéarisée de H en X appliqué à dX
171 # ------------------------------------------------------
172 _dX = numpy.asmatrix(numpy.ravel( dX )).T
173 _HtX = numpy.dot(_Jacobienne, _dX)
176 # ---------------------------------------------------------
177 def AdjointOperator(self, (X, Y) ):
179 Calcul de l'adjoint à l'aide de la Jacobienne.
181 _JacobienneT = self.TangentMatrix( X ).T
182 if Y is None or len(Y) == 0:
184 # Calcul de la forme matricielle si le second argument est None
185 # -------------------------------------------------------------
189 # Calcul de la valeur de l'adjoint en X appliqué à Y
190 # --------------------------------------------------
191 _Y = numpy.asmatrix(numpy.ravel( Y )).T
192 _HaY = numpy.dot(_JacobienneT, _Y)
195 # ==============================================================================
198 """ Direct non-linear simulation operator """
200 # NEED TO BE COMPLETED
201 # NEED TO BE COMPLETED
202 # NEED TO BE COMPLETED
204 # --------------------------------------> # EXAMPLE TO BE REMOVED
205 # Example of Identity operator # EXAMPLE TO BE REMOVED
206 if type(XX) is type(numpy.matrix([])): # EXAMPLE TO BE REMOVED
207 HX = XX.A1.tolist() # EXAMPLE TO BE REMOVED
208 elif type(XX) is type(numpy.array([])): # EXAMPLE TO BE REMOVED
209 HX = numpy.matrix(XX).A1.tolist() # EXAMPLE TO BE REMOVED
210 else: # EXAMPLE TO BE REMOVED
211 HX = XX # EXAMPLE TO BE REMOVED
212 # # EXAMPLE TO BE REMOVED
213 HHX = [] # EXAMPLE TO BE REMOVED
214 HHX.extend( HX ) # EXAMPLE TO BE REMOVED
215 HHX.extend( HX ) # EXAMPLE TO BE REMOVED
216 # --------------------------------------> # EXAMPLE TO BE REMOVED
218 return numpy.array( HHX )
220 # ==============================================================================
221 if __name__ == "__main__":
224 print "AUTODIAGNOSTIC"
225 print "=============="
229 FDA = FDApproximation( test1 )
230 print "H(X) =", FDA.DirectOperator( X0 )
231 print "Tg matrice =\n", FDA.TangentMatrix( X0 )
232 print "Tg(X) =", FDA.TangentOperator( (X0, X0) )
233 print "Ad((X,Y)) =", FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
236 FDA = FDApproximation( test1, centeredDF=True )
237 print "H(X) =", FDA.DirectOperator( X0 )
238 print "Tg matrice =\n", FDA.TangentMatrix( X0 )
239 print "Tg(X) =", FDA.TangentOperator( (X0, X0) )
240 print "Ad((X,Y)) =", FDA.AdjointOperator( (X0,range(3,3+2*len(X0))) )
244 print "=============="
248 FDA = FDApproximation( test1 )
249 print "H(X) =", FDA.DirectOperator( X0 )
250 print "Tg matrice =\n", FDA.TangentMatrix( X0 )
251 print "Tg(X) =", FDA.TangentOperator( (X0, X0) )
252 print "Ad((X,Y)) =", FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
255 FDA = FDApproximation( test1, centeredDF=True )
256 print "H(X) =", FDA.DirectOperator( X0 )
257 print "Tg matrice =\n", FDA.TangentMatrix( X0 )
258 print "Tg(X) =", FDA.TangentOperator( (X0, X0) )
259 print "Ad((X,Y)) =", FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
263 print "=============="
267 FDA = FDApproximation( test1 )
268 print "H(X) =", FDA.DirectOperator( X0 )
269 print "Tg matrice =\n", FDA.TangentMatrix( X0 )
270 print "Tg(X) =", FDA.TangentOperator( (X0, X0) )
271 print "Ad((X,Y)) =", FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
274 FDA = FDApproximation( test1, centeredDF=True )
275 print "H(X) =", FDA.DirectOperator( X0 )
276 print "Tg matrice =\n", FDA.TangentMatrix( X0 )
277 print "Tg(X) =", FDA.TangentOperator( (X0, X0) )
278 print "Ad((X,Y)) =", FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
282 print "=============="
284 X0 = numpy.asmatrix(numpy.arange(4)).T
286 FDA = FDApproximation( test1 )
287 print "H(X) =", FDA.DirectOperator( X0 )
288 print "Tg matrice =\n", FDA.TangentMatrix( X0 )
289 print "Tg(X) =", FDA.TangentOperator( (X0, X0) )
290 print "Ad((X,Y)) =", FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )
293 FDA = FDApproximation( test1, centeredDF=True )
294 print "H(X) =", FDA.DirectOperator( X0 )
295 print "Tg matrice =\n", FDA.TangentMatrix( X0 )
296 print "Tg(X) =", FDA.TangentOperator( (X0, X0) )
297 print "Ad((X,Y)) =", FDA.AdjointOperator( (X0,range(7,7+2*len(X0))) )