1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2022 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 objets numériques génériques.
26 __author__ = "Jean-Philippe ARGAUD"
28 import os, copy, types, sys, logging, numpy
29 from daCore.BasicObjects import Operator, Covariance, PartialAlgorithm
30 from daCore.PlatformInfo import PlatformInfo
31 mpr = PlatformInfo().MachinePrecision()
32 mfp = PlatformInfo().MaximumPrecision()
33 # logging.getLogger().setLevel(logging.DEBUG)
35 # ==============================================================================
36 def ExecuteFunction( triplet ):
37 assert len(triplet) == 3, "Incorrect number of arguments"
38 X, xArgs, funcrepr = triplet
39 __X = numpy.ravel( X ).reshape((-1,1))
40 __sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"])
41 __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
42 __fonction = getattr(__module,funcrepr["__userFunction__name"])
43 sys.path = __sys_path_tmp ; del __sys_path_tmp
44 if isinstance(xArgs, dict):
45 __HX = __fonction( __X, **xArgs )
47 __HX = __fonction( __X )
48 return numpy.ravel( __HX )
50 # ==============================================================================
51 class FDApproximation(object):
53 Cette classe sert d'interface pour définir les opérateurs approximés. A la
54 création d'un objet, en fournissant une fonction "Function", on obtient un
55 objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
56 "AdjointOperator". On contrôle l'approximation DF avec l'incrément
57 multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
58 "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
59 centrées si le booléen "centeredDF" est vrai.
62 name = "FDApproximation",
67 extraArguments = None,
68 reducingMemoryUse = False,
69 avoidingRedundancy = True,
70 toleranceInRedundancy = 1.e-18,
71 lenghtOfRedundancy = -1,
76 self.__name = str(name)
77 self.__extraArgs = extraArguments
81 import multiprocessing
82 self.__mpEnabled = True
84 self.__mpEnabled = False
86 self.__mpEnabled = False
87 self.__mpWorkers = mpWorkers
88 if self.__mpWorkers is not None and self.__mpWorkers < 1:
89 self.__mpWorkers = None
90 logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
92 self.__mfEnabled = bool(mfEnabled)
93 logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
95 self.__rmEnabled = bool(reducingMemoryUse)
96 logging.debug("FDA Calculs avec réduction mémoire : %s"%(self.__rmEnabled,))
98 if avoidingRedundancy:
100 self.__tolerBP = float(toleranceInRedundancy)
101 self.__lenghtRJ = int(lenghtOfRedundancy)
102 self.__listJPCP = [] # Jacobian Previous Calculated Points
103 self.__listJPCI = [] # Jacobian Previous Calculated Increment
104 self.__listJPCR = [] # Jacobian Previous Calculated Results
105 self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
106 self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
108 self.__avoidRC = False
109 logging.debug("FDA Calculs avec réduction des doublons : %s"%self.__avoidRC)
111 logging.debug("FDA Tolérance de détermination des doublons : %.2e"%self.__tolerBP)
114 if isinstance(Function,types.FunctionType):
115 logging.debug("FDA Calculs en multiprocessing : FunctionType")
116 self.__userFunction__name = Function.__name__
118 mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
120 mod = os.path.abspath(Function.__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)
126 self.__userOperator = Operator(
128 fromMethod = Function,
129 avoidingRedundancy = self.__avoidRC,
130 inputAsMultiFunction = self.__mfEnabled,
131 extraArguments = self.__extraArgs )
132 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
133 elif isinstance(Function,types.MethodType):
134 logging.debug("FDA Calculs en multiprocessing : MethodType")
135 self.__userFunction__name = Function.__name__
137 mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
139 mod = os.path.abspath(Function.__func__.__globals__['__file__'])
140 if not os.path.isfile(mod):
141 raise ImportError("No user defined function or method found with the name %s"%(mod,))
142 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
143 self.__userFunction__path = os.path.dirname(mod)
145 self.__userOperator = Operator(
147 fromMethod = Function,
148 avoidingRedundancy = self.__avoidRC,
149 inputAsMultiFunction = self.__mfEnabled,
150 extraArguments = self.__extraArgs )
151 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
153 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
155 self.__userOperator = Operator(
157 fromMethod = Function,
158 avoidingRedundancy = self.__avoidRC,
159 inputAsMultiFunction = self.__mfEnabled,
160 extraArguments = self.__extraArgs )
161 self.__userFunction = self.__userOperator.appliedTo
163 self.__centeredDF = bool(centeredDF)
164 if abs(float(increment)) > 1.e-15:
165 self.__increment = float(increment)
167 self.__increment = 0.01
171 self.__dX = numpy.ravel( dX )
173 # ---------------------------------------------------------
174 def __doublon__(self, __e, __l, __n, __v=None):
175 __ac, __iac = False, -1
176 for i in range(len(__l)-1,-1,-1):
177 if numpy.linalg.norm(__e - __l[i]) < self.__tolerBP * __n[i]:
178 __ac, __iac = True, i
179 if __v is not None: logging.debug("FDA Cas%s déjà calculé, récupération du doublon %i"%(__v,__iac))
183 # ---------------------------------------------------------
184 def __listdotwith__(self, __LMatrix, __dotWith = None, __dotTWith = None):
185 "Produit incrémental d'une matrice liste de colonnes avec un vecteur"
186 if not isinstance(__LMatrix, (list,tuple)):
187 raise TypeError("Columnwise list matrix has not the proper type: %s"%type(__LMatrix))
188 if __dotWith is not None:
189 __Idwx = numpy.ravel( __dotWith )
190 assert len(__LMatrix) == __Idwx.size, "Incorrect size of elements"
191 __Produit = numpy.zeros(__LMatrix[0].size)
192 for i, col in enumerate(__LMatrix):
193 __Produit += float(__Idwx[i]) * col
195 elif __dotTWith is not None:
196 _Idwy = numpy.ravel( __dotTWith ).T
197 assert __LMatrix[0].size == _Idwy.size, "Incorrect size of elements"
198 __Produit = numpy.zeros(len(__LMatrix))
199 for i, col in enumerate(__LMatrix):
200 __Produit[i] = float( _Idwy @ col)
206 # ---------------------------------------------------------
207 def DirectOperator(self, X, **extraArgs ):
209 Calcul du direct à l'aide de la fonction fournie.
211 NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
212 ne doivent pas être données ici à la fonction utilisateur.
214 logging.debug("FDA Calcul DirectOperator (explicite)")
216 _HX = self.__userFunction( X, argsAsSerie = True )
218 _HX = numpy.ravel(self.__userFunction( numpy.ravel(X) ))
222 # ---------------------------------------------------------
223 def TangentMatrix(self, X, dotWith = None, dotTWith = None ):
225 Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
226 c'est-à-dire le gradient de H en X. On utilise des différences finies
227 directionnelles autour du point X. X est un numpy.ndarray.
229 Différences finies centrées (approximation d'ordre 2):
230 1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
231 dX[i] à la composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
232 on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
234 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
236 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
238 Différences finies non centrées (approximation d'ordre 1):
239 1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
240 composante X[i] pour composer X_plus_dXi, et on calcule la réponse
241 HX_plus_dXi = H( X_plus_dXi )
242 2/ On calcule la valeur centrale HX = H(X)
243 3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
245 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
248 logging.debug("FDA Début du calcul de la Jacobienne")
249 logging.debug("FDA Incrément de............: %s*X"%float(self.__increment))
250 logging.debug("FDA Approximation centrée...: %s"%(self.__centeredDF))
252 if X is None or len(X)==0:
253 raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
255 _X = numpy.ravel( X )
257 if self.__dX is None:
258 _dX = self.__increment * _X
260 _dX = numpy.ravel( self.__dX )
261 assert len(_X) == len(_dX), "Inconsistent dX increment length with respect to the X one"
262 assert _X.size == _dX.size, "Inconsistent dX increment size with respect to the X one"
264 if (_dX == 0.).any():
267 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
269 _dX = numpy.where( _dX == 0., moyenne, _dX )
271 __alreadyCalculated = False
273 __bidon, __alreadyCalculatedP = self.__doublon__(_X, self.__listJPCP, self.__listJPPN, None)
274 __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
275 if __alreadyCalculatedP == __alreadyCalculatedI > -1:
276 __alreadyCalculated, __i = True, __alreadyCalculatedP
277 logging.debug("FDA Cas J déjà calculé, récupération du doublon %i"%__i)
279 if __alreadyCalculated:
280 logging.debug("FDA Calcul Jacobienne (par récupération du doublon %i)"%__i)
281 _Jacobienne = self.__listJPCR[__i]
282 logging.debug("FDA Fin du calcul de la Jacobienne")
283 if dotWith is not None:
284 return numpy.dot(_Jacobienne, numpy.ravel( dotWith ))
285 elif dotTWith is not None:
286 return numpy.dot(_Jacobienne.T, numpy.ravel( dotTWith ))
288 logging.debug("FDA Calcul Jacobienne (explicite)")
289 if self.__centeredDF:
291 if self.__mpEnabled and not self.__mfEnabled:
293 "__userFunction__path" : self.__userFunction__path,
294 "__userFunction__modl" : self.__userFunction__modl,
295 "__userFunction__name" : self.__userFunction__name,
298 for i in range( len(_dX) ):
300 _X_plus_dXi = numpy.array( _X, dtype=float )
301 _X_plus_dXi[i] = _X[i] + _dXi
302 _X_moins_dXi = numpy.array( _X, dtype=float )
303 _X_moins_dXi[i] = _X[i] - _dXi
305 _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
306 _jobs.append( (_X_moins_dXi, self.__extraArgs, funcrepr) )
308 import multiprocessing
309 self.__pool = multiprocessing.Pool(self.__mpWorkers)
310 _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
315 for i in range( len(_dX) ):
316 _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
318 elif self.__mfEnabled:
320 for i in range( len(_dX) ):
322 _X_plus_dXi = numpy.array( _X, dtype=float )
323 _X_plus_dXi[i] = _X[i] + _dXi
324 _X_moins_dXi = numpy.array( _X, dtype=float )
325 _X_moins_dXi[i] = _X[i] - _dXi
327 _xserie.append( _X_plus_dXi )
328 _xserie.append( _X_moins_dXi )
330 _HX_plusmoins_dX = self.DirectOperator( _xserie )
333 for i in range( len(_dX) ):
334 _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
338 for i in range( _dX.size ):
340 _X_plus_dXi = numpy.array( _X, dtype=float )
341 _X_plus_dXi[i] = _X[i] + _dXi
342 _X_moins_dXi = numpy.array( _X, dtype=float )
343 _X_moins_dXi[i] = _X[i] - _dXi
345 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
346 _HX_moins_dXi = self.DirectOperator( _X_moins_dXi )
348 _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
352 if self.__mpEnabled and not self.__mfEnabled:
354 "__userFunction__path" : self.__userFunction__path,
355 "__userFunction__modl" : self.__userFunction__modl,
356 "__userFunction__name" : self.__userFunction__name,
359 _jobs.append( (_X, self.__extraArgs, funcrepr) )
360 for i in range( len(_dX) ):
361 _X_plus_dXi = numpy.array( _X, dtype=float )
362 _X_plus_dXi[i] = _X[i] + _dX[i]
364 _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
366 import multiprocessing
367 self.__pool = multiprocessing.Pool(self.__mpWorkers)
368 _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
372 _HX = _HX_plus_dX.pop(0)
375 for i in range( len(_dX) ):
376 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
378 elif self.__mfEnabled:
381 for i in range( len(_dX) ):
382 _X_plus_dXi = numpy.array( _X, dtype=float )
383 _X_plus_dXi[i] = _X[i] + _dX[i]
385 _xserie.append( _X_plus_dXi )
387 _HX_plus_dX = self.DirectOperator( _xserie )
389 _HX = _HX_plus_dX.pop(0)
392 for i in range( len(_dX) ):
393 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
397 _HX = self.DirectOperator( _X )
398 for i in range( _dX.size ):
400 _X_plus_dXi = numpy.array( _X, dtype=float )
401 _X_plus_dXi[i] = _X[i] + _dXi
403 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
405 _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
407 if (dotWith is not None) or (dotTWith is not None):
408 __Produit = self.__listdotwith__(_Jacobienne, dotWith, dotTWith)
411 if __Produit is None or self.__avoidRC:
412 _Jacobienne = numpy.transpose( numpy.vstack( _Jacobienne ) )
414 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
415 while len(self.__listJPCP) > self.__lenghtRJ:
416 self.__listJPCP.pop(0)
417 self.__listJPCI.pop(0)
418 self.__listJPCR.pop(0)
419 self.__listJPPN.pop(0)
420 self.__listJPIN.pop(0)
421 self.__listJPCP.append( copy.copy(_X) )
422 self.__listJPCI.append( copy.copy(_dX) )
423 self.__listJPCR.append( copy.copy(_Jacobienne) )
424 self.__listJPPN.append( numpy.linalg.norm(_X) )
425 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
426 logging.debug("FDA Fin du calcul de la Jacobienne")
427 if __Produit is not None:
432 # ---------------------------------------------------------
433 def TangentOperator(self, paire, **extraArgs ):
435 Calcul du tangent à l'aide de la Jacobienne.
437 NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
438 ne doivent pas être données ici à la fonction utilisateur.
441 assert len(paire) == 1, "Incorrect length of arguments"
443 assert len(_paire) == 2, "Incorrect number of arguments"
445 assert len(paire) == 2, "Incorrect number of arguments"
448 if dX is None or len(dX) == 0:
450 # Calcul de la forme matricielle si le second argument est None
451 # -------------------------------------------------------------
452 _Jacobienne = self.TangentMatrix( X )
453 if self.__mfEnabled: return [_Jacobienne,]
454 else: return _Jacobienne
457 # Calcul de la valeur linéarisée de H en X appliqué à dX
458 # ------------------------------------------------------
459 _HtX = self.TangentMatrix( X, dotWith = dX )
460 if self.__mfEnabled: return [_HtX,]
463 # ---------------------------------------------------------
464 def AdjointOperator(self, paire, **extraArgs ):
466 Calcul de l'adjoint à l'aide de la Jacobienne.
468 NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
469 ne doivent pas être données ici à la fonction utilisateur.
472 assert len(paire) == 1, "Incorrect length of arguments"
474 assert len(_paire) == 2, "Incorrect number of arguments"
476 assert len(paire) == 2, "Incorrect number of arguments"
479 if Y is None or len(Y) == 0:
481 # Calcul de la forme matricielle si le second argument est None
482 # -------------------------------------------------------------
483 _JacobienneT = self.TangentMatrix( X ).T
484 if self.__mfEnabled: return [_JacobienneT,]
485 else: return _JacobienneT
488 # Calcul de la valeur de l'adjoint en X appliqué à Y
489 # --------------------------------------------------
490 _HaY = self.TangentMatrix( X, dotTWith = Y )
491 if self.__mfEnabled: return [_HaY,]
494 # ==============================================================================
495 def SetInitialDirection( __Direction = [], __Amplitude = 1., __Position = None ):
496 "Établit ou élabore une direction avec une amplitude"
498 if len(__Direction) == 0 and __Position is None:
499 raise ValueError("If initial direction is void, current position has to be given")
500 if abs(float(__Amplitude)) < mpr:
501 raise ValueError("Amplitude of perturbation can not be zero")
503 if len(__Direction) > 0:
504 __dX0 = numpy.asarray(__Direction)
507 __X0 = numpy.ravel(numpy.asarray(__Position))
508 __mX0 = numpy.mean( __X0, dtype=mfp )
509 if abs(__mX0) < 2*mpr: __mX0 = 1. # Évite le problème de position nulle
512 __dX0.append( numpy.random.normal(0.,abs(v)) )
514 __dX0.append( numpy.random.normal(0.,__mX0) )
516 __dX0 = numpy.asarray(__dX0,float) # Évite le problème d'array de taille 1
517 __dX0 = numpy.ravel( __dX0 ) # Redresse les vecteurs
518 __dX0 = float(__Amplitude) * __dX0
522 # ==============================================================================
523 def EnsembleOfCenteredPerturbations( __bgCenter, __bgCovariance, __nbMembers ):
524 "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
526 __bgCenter = numpy.ravel(__bgCenter)[:,None]
528 raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
530 if __bgCovariance is None:
531 _Perturbations = numpy.tile( __bgCenter, __nbMembers)
533 _Z = numpy.random.multivariate_normal(numpy.zeros(__bgCenter.size), __bgCovariance, size=__nbMembers).T
534 _Perturbations = numpy.tile( __bgCenter, __nbMembers) + _Z
536 return _Perturbations
538 # ==============================================================================
539 def EnsembleOfBackgroundPerturbations( __bgCenter, __bgCovariance, __nbMembers, __withSVD = True):
540 "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
541 def __CenteredRandomAnomalies(Zr, N):
543 Génère une matrice de N anomalies aléatoires centrées sur Zr selon les
544 notes manuscrites de MB et conforme au code de PS avec eps = -1
547 Q = numpy.identity(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
548 Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
549 R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
554 __bgCenter = numpy.ravel(__bgCenter).reshape((-1,1))
556 raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
557 if __bgCovariance is None:
558 _Perturbations = numpy.tile( __bgCenter, __nbMembers)
561 _U, _s, _V = numpy.linalg.svd(__bgCovariance, full_matrices=False)
562 _nbctl = __bgCenter.size
563 if __nbMembers > _nbctl:
564 _Z = numpy.concatenate((numpy.dot(
565 numpy.diag(numpy.sqrt(_s[:_nbctl])), _V[:_nbctl]),
566 numpy.random.multivariate_normal(numpy.zeros(_nbctl),__bgCovariance,__nbMembers-1-_nbctl)), axis = 0)
568 _Z = numpy.dot(numpy.diag(numpy.sqrt(_s[:__nbMembers-1])), _V[:__nbMembers-1])
569 _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
570 _Perturbations = __bgCenter + _Zca
572 if max(abs(__bgCovariance.flatten())) > 0:
573 _nbctl = __bgCenter.size
574 _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),__bgCovariance,__nbMembers-1)
575 _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
576 _Perturbations = __bgCenter + _Zca
578 _Perturbations = numpy.tile( __bgCenter, __nbMembers)
580 return _Perturbations
582 # ==============================================================================
583 def EnsembleMean( __Ensemble ):
584 "Renvoie la moyenne empirique d'un ensemble"
585 return numpy.asarray(__Ensemble).mean(axis=1, dtype=mfp).astype('float').reshape((-1,1))
587 # ==============================================================================
588 def EnsembleOfAnomalies( __Ensemble, __OptMean = None, __Normalisation = 1.):
589 "Renvoie les anomalies centrées à partir d'un ensemble"
590 if __OptMean is None:
591 __Em = EnsembleMean( __Ensemble )
593 __Em = numpy.ravel( __OptMean ).reshape((-1,1))
595 return __Normalisation * (numpy.asarray( __Ensemble ) - __Em)
597 # ==============================================================================
598 def EnsembleErrorCovariance( __Ensemble, __Quick = False ):
599 "Renvoie l'estimation empirique de la covariance d'ensemble"
601 # Covariance rapide mais rarement définie positive
602 __Covariance = numpy.cov( __Ensemble )
604 # Résultat souvent identique à numpy.cov, mais plus robuste
605 __n, __m = numpy.asarray( __Ensemble ).shape
606 __Anomalies = EnsembleOfAnomalies( __Ensemble )
607 # Estimation empirique
608 __Covariance = ( __Anomalies @ __Anomalies.T ) / (__m-1)
610 __Covariance = ( __Covariance + __Covariance.T ) * 0.5
611 # Assure la positivité
612 __epsilon = mpr*numpy.trace( __Covariance )
613 __Covariance = __Covariance + __epsilon * numpy.identity(__n)
617 # ==============================================================================
618 def EnsemblePerturbationWithGivenCovariance(
623 "Ajout d'une perturbation à chaque membre d'un ensemble selon une covariance prescrite"
624 if hasattr(__Covariance,"assparsematrix"):
625 if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance.assparsematrix())/abs(__Ensemble).mean() < mpr).all():
626 # Traitement d'une covariance nulle ou presque
628 if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance.assparsematrix()) < mpr).all():
629 # Traitement d'une covariance nulle ou presque
632 if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance)/abs(__Ensemble).mean() < mpr).all():
633 # Traitement d'une covariance nulle ou presque
635 if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance) < mpr).all():
636 # Traitement d'une covariance nulle ou presque
639 __n, __m = __Ensemble.shape
640 if __Seed is not None: numpy.random.seed(__Seed)
642 if hasattr(__Covariance,"isscalar") and __Covariance.isscalar():
643 # Traitement d'une covariance multiple de l'identité
645 __std = numpy.sqrt(__Covariance.assparsematrix())
646 __Ensemble += numpy.random.normal(__zero, __std, size=(__m,__n)).T
648 elif hasattr(__Covariance,"isvector") and __Covariance.isvector():
649 # Traitement d'une covariance diagonale avec variances non identiques
650 __zero = numpy.zeros(__n)
651 __std = numpy.sqrt(__Covariance.assparsematrix())
652 __Ensemble += numpy.asarray([numpy.random.normal(__zero, __std) for i in range(__m)]).T
654 elif hasattr(__Covariance,"ismatrix") and __Covariance.ismatrix():
655 # Traitement d'une covariance pleine
656 __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance.asfullmatrix(__n), size=__m).T
658 elif isinstance(__Covariance, numpy.ndarray):
659 # Traitement d'une covariance numpy pleine, sachant qu'on arrive ici en dernier
660 __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance, size=__m).T
663 raise ValueError("Error in ensemble perturbation with inadequate covariance specification")
667 # ==============================================================================
668 def CovarianceInflation(
670 __InflationType = None,
671 __InflationFactor = None,
672 __BackgroundCov = None,
675 Inflation applicable soit sur Pb ou Pa, soit sur les ensembles EXb ou EXa
677 Synthèse : Hunt 2007, section 2.3.5
679 if __InflationFactor is None:
680 return __InputCovOrEns
682 __InflationFactor = float(__InflationFactor)
684 __InputCovOrEns = numpy.asarray(__InputCovOrEns)
685 if __InputCovOrEns.size == 0: return __InputCovOrEns
687 if __InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
688 if __InflationFactor < 1.:
689 raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
690 if __InflationFactor < 1.+mpr: # No inflation = 1
691 return __InputCovOrEns
692 __OutputCovOrEns = __InflationFactor**2 * __InputCovOrEns
694 elif __InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
695 if __InflationFactor < 1.:
696 raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
697 if __InflationFactor < 1.+mpr: # No inflation = 1
698 return __InputCovOrEns
699 __InputCovOrEnsMean = __InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
700 __OutputCovOrEns = __InputCovOrEnsMean[:,numpy.newaxis] \
701 + __InflationFactor * (__InputCovOrEns - __InputCovOrEnsMean[:,numpy.newaxis])
703 elif __InflationType in ["AdditiveOnAnalysisCovariance", "AdditiveOnBackgroundCovariance"]:
704 if __InflationFactor < 0.:
705 raise ValueError("Inflation factor for additive inflation has to be greater or equal than 0.")
706 if __InflationFactor < mpr: # No inflation = 0
707 return __InputCovOrEns
708 __n, __m = __InputCovOrEns.shape
710 raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
711 __tr = __InputCovOrEns.trace()/__n
712 if __InflationFactor > __tr:
713 raise ValueError("Inflation factor for additive inflation has to be small over %.0e."%__tr)
714 __OutputCovOrEns = (1. - __InflationFactor)*__InputCovOrEns + __InflationFactor * numpy.identity(__n)
716 elif __InflationType == "HybridOnBackgroundCovariance":
717 if __InflationFactor < 0.:
718 raise ValueError("Inflation factor for hybrid inflation has to be greater or equal than 0.")
719 if __InflationFactor < mpr: # No inflation = 0
720 return __InputCovOrEns
721 __n, __m = __InputCovOrEns.shape
723 raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
724 if __BackgroundCov is None:
725 raise ValueError("Background covariance matrix B has to be given for hybrid inflation.")
726 if __InputCovOrEns.shape != __BackgroundCov.shape:
727 raise ValueError("Ensemble covariance matrix has to be of same size than background covariance matrix B.")
728 __OutputCovOrEns = (1. - __InflationFactor) * __InputCovOrEns + __InflationFactor * __BackgroundCov
730 elif __InflationType == "Relaxation":
731 raise NotImplementedError("Relaxation inflation type not implemented")
734 raise ValueError("Error in inflation type, '%s' is not a valid keyword."%__InflationType)
736 return __OutputCovOrEns
738 # ==============================================================================
739 def HessienneEstimation(__selfA, __nb, __HaM, __HtM, __BI, __RI):
740 "Estimation de la Hessienne"
743 for i in range(int(__nb)):
744 __ee = numpy.zeros((__nb,1))
746 __HtEE = numpy.dot(__HtM,__ee).reshape((-1,1))
747 __HessienneI.append( numpy.ravel( __BI * __ee + __HaM * (__RI * __HtEE) ) )
749 __A = numpy.linalg.inv(numpy.array( __HessienneI ))
750 __A = (__A + __A.T) * 0.5 # Symétrie
751 __A = __A + mpr*numpy.trace( __A ) * numpy.identity(__nb) # Positivité
753 if min(__A.shape) != max(__A.shape):
755 "The %s a posteriori covariance matrix A"%(__selfA._name,)+\
756 " is of shape %s, despites it has to be a"%(str(__A.shape),)+\
757 " squared matrix. There is an error in the observation operator,"+\
759 if (numpy.diag(__A) < 0).any():
761 "The %s a posteriori covariance matrix A"%(__selfA._name,)+\
762 " has at least one negative value on its diagonal. There is an"+\
763 " error in the observation operator, please check it.")
764 if logging.getLogger().level < logging.WARNING: # La vérification n'a lieu qu'en debug
766 numpy.linalg.cholesky( __A )
769 "The %s a posteriori covariance matrix A"%(__selfA._name,)+\
770 " is not symmetric positive-definite. Please check your a"+\
771 " priori covariances and your observation operator.")
775 # ==============================================================================
776 def QuantilesEstimations(selfA, A, Xa, HXa = None, Hm = None, HtM = None):
777 "Estimation des quantiles a posteriori à partir de A>0 (selfA est modifié)"
778 nbsamples = selfA._parameters["NumberOfSamplesForQuantiles"]
780 # Traitement des bornes
781 if "StateBoundsForQuantiles" in selfA._parameters:
782 LBounds = selfA._parameters["StateBoundsForQuantiles"] # Prioritaire
783 elif "Bounds" in selfA._parameters:
784 LBounds = selfA._parameters["Bounds"] # Défaut raisonnable
787 if LBounds is not None:
788 LBounds = ForceNumericBounds( LBounds )
789 __Xa = numpy.ravel(Xa)
791 # Échantillonnage des états
794 for i in range(nbsamples):
795 if selfA._parameters["SimulationForQuantiles"] == "Linear" and HtM is not None and HXa is not None:
796 dXr = (numpy.random.multivariate_normal(__Xa,A) - __Xa).reshape((-1,1))
797 if LBounds is not None: # "EstimateProjection" par défaut
798 dXr = numpy.max(numpy.hstack((dXr,LBounds[:,0].reshape((-1,1))) - __Xa.reshape((-1,1))),axis=1)
799 dXr = numpy.min(numpy.hstack((dXr,LBounds[:,1].reshape((-1,1))) - __Xa.reshape((-1,1))),axis=1)
801 Yr = HXa.reshape((-1,1)) + dYr
802 if selfA._toStore("SampledStateForQuantiles"): Xr = __Xa + numpy.ravel(dXr)
803 elif selfA._parameters["SimulationForQuantiles"] == "NonLinear" and Hm is not None:
804 Xr = numpy.random.multivariate_normal(__Xa,A)
805 if LBounds is not None: # "EstimateProjection" par défaut
806 Xr = numpy.max(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,0].reshape((-1,1)))),axis=1)
807 Xr = numpy.min(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,1].reshape((-1,1)))),axis=1)
808 Yr = numpy.asarray(Hm( Xr ))
810 raise ValueError("Quantile simulations has only to be Linear or NonLinear.")
813 YfQ = Yr.reshape((-1,1))
814 if selfA._toStore("SampledStateForQuantiles"): EXr = Xr.reshape((-1,1))
816 YfQ = numpy.hstack((YfQ,Yr.reshape((-1,1))))
817 if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.hstack((EXr,Xr.reshape((-1,1))))
819 # Extraction des quantiles
822 for quantile in selfA._parameters["Quantiles"]:
823 if not (0. <= float(quantile) <= 1.): continue
824 indice = int(nbsamples * float(quantile) - 1./nbsamples)
825 if YQ is None: YQ = YfQ[:,indice].reshape((-1,1))
826 else: YQ = numpy.hstack((YQ,YfQ[:,indice].reshape((-1,1))))
827 if YQ is not None: # Liste non vide de quantiles
828 selfA.StoredVariables["SimulationQuantiles"].store( YQ )
829 if selfA._toStore("SampledStateForQuantiles"):
830 selfA.StoredVariables["SampledStateForQuantiles"].store( EXr )
834 # ==============================================================================
835 def ForceNumericBounds( __Bounds, __infNumbers = True ):
836 "Force les bornes à être des valeurs numériques, sauf si globalement None"
837 # Conserve une valeur par défaut à None s'il n'y a pas de bornes
838 if __Bounds is None: return None
839 # Converti toutes les bornes individuelles None à +/- l'infini chiffré
840 __Bounds = numpy.asarray( __Bounds, dtype=float )
841 if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
842 raise ValueError("Incorrectly shaped bounds data")
844 __Bounds[numpy.isnan(__Bounds[:,0]),0] = -float('inf')
845 __Bounds[numpy.isnan(__Bounds[:,1]),1] = float('inf')
847 __Bounds[numpy.isnan(__Bounds[:,0]),0] = -sys.float_info.max
848 __Bounds[numpy.isnan(__Bounds[:,1]),1] = sys.float_info.max
851 # ==============================================================================
852 def RecentredBounds( __Bounds, __Center, __Scale = None):
853 "Recentre les bornes autour de 0, sauf si globalement None"
854 # Conserve une valeur par défaut à None s'il n'y a pas de bornes
855 if __Bounds is None: return None
857 # Recentre les valeurs numériques de bornes
858 return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1,1))
860 # Recentre les valeurs numériques de bornes et change l'échelle par une matrice
861 return __Scale @ (ForceNumericBounds( __Bounds, False ) - numpy.ravel( __Center ).reshape((-1,1)))
863 # ==============================================================================
864 def ApplyBounds( __Vector, __Bounds, __newClip = True):
865 "Applique des bornes numériques à un point"
866 # Conserve une valeur par défaut s'il n'y a pas de bornes
867 if __Bounds is None: return __Vector
869 if not isinstance(__Vector, numpy.ndarray): # Is an array
870 raise ValueError("Incorrect array definition of vector data")
871 if not isinstance(__Bounds, numpy.ndarray): # Is an array
872 raise ValueError("Incorrect array definition of bounds data")
873 if 2*__Vector.size != __Bounds.size: # Is a 2 column array of vector lenght
874 raise ValueError("Incorrect bounds number (%i) to be applied for this vector (of size %i)"%(__Bounds.size,__Vector.size))
875 if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
876 raise ValueError("Incorrectly shaped bounds data")
879 __Vector = __Vector.clip(
880 __Bounds[:,0].reshape(__Vector.shape),
881 __Bounds[:,1].reshape(__Vector.shape),
884 __Vector = numpy.max(numpy.hstack((__Vector.reshape((-1,1)),numpy.asmatrix(__Bounds)[:,0])),axis=1)
885 __Vector = numpy.min(numpy.hstack((__Vector.reshape((-1,1)),numpy.asmatrix(__Bounds)[:,1])),axis=1)
886 __Vector = numpy.asarray(__Vector)
890 # ==============================================================================
891 def Apply3DVarRecentringOnEnsemble(__EnXn, __EnXf, __Ynpu, __HO, __R, __B, __SuppPars):
892 "Recentre l'ensemble Xn autour de l'analyse 3DVAR"
893 __Betaf = __SuppPars["HybridCovarianceEquilibrium"]
895 Xf = EnsembleMean( __EnXf )
896 Pf = Covariance( asCovariance=EnsembleErrorCovariance(__EnXf) )
897 Pf = (1 - __Betaf) * __B.asfullmatrix(Xf.size) + __Betaf * Pf
899 selfB = PartialAlgorithm("3DVAR")
900 selfB._parameters["Minimizer"] = "LBFGSB"
901 selfB._parameters["MaximumNumberOfIterations"] = __SuppPars["HybridMaximumNumberOfIterations"]
902 selfB._parameters["CostDecrementTolerance"] = __SuppPars["HybridCostDecrementTolerance"]
903 selfB._parameters["ProjectedGradientTolerance"] = -1
904 selfB._parameters["GradientNormTolerance"] = 1.e-05
905 selfB._parameters["StoreInternalVariables"] = False
906 selfB._parameters["optiprint"] = -1
907 selfB._parameters["optdisp"] = 0
908 selfB._parameters["Bounds"] = None
909 selfB._parameters["InitializationPoint"] = Xf
910 from daAlgorithms.Atoms import std3dvar
911 std3dvar.std3dvar(selfB, Xf, __Ynpu, None, __HO, None, __R, Pf)
912 Xa = selfB.get("Analysis")[-1].reshape((-1,1))
915 return Xa + EnsembleOfAnomalies( __EnXn )
917 # ==============================================================================
918 def multiXOsteps(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle,
919 __CovForecast = False, __LinEvolution = False,
922 Prévision multi-pas avec une correction par pas (multi-méthodes)
927 if selfA._parameters["EstimationOf"] == "State":
928 if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
929 Xn = numpy.asarray(Xb)
930 if __CovForecast: Pn = B
931 selfA.StoredVariables["Analysis"].store( Xn )
932 if selfA._toStore("APosterioriCovariance"):
933 if hasattr(B,"asfullmatrix"):
934 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(Xn.size) )
936 selfA.StoredVariables["APosterioriCovariance"].store( B )
937 selfA._setInternalState("seed", numpy.random.get_state())
938 elif selfA._parameters["nextStep"]:
939 Xn = selfA._getInternalState("Xn")
940 if __CovForecast: Pn = selfA._getInternalState("Pn")
942 Xn = numpy.asarray(Xb)
943 if __CovForecast: Pn = B
945 if hasattr(Y,"stepnumber"):
946 duration = Y.stepnumber()
952 for step in range(duration-1):
953 selfA.StoredVariables["CurrentStepNumber"].store( len(selfA.StoredVariables["Analysis"]) )
955 if hasattr(Y,"store"):
956 Ynpu = numpy.asarray( Y[step+1] ).reshape((-1,1))
958 Ynpu = numpy.asarray( Y ).reshape((-1,1))
961 if hasattr(U,"store") and len(U)>1:
962 Un = numpy.asarray( U[step] ).reshape((-1,1))
963 elif hasattr(U,"store") and len(U)==1:
964 Un = numpy.asarray( U[0] ).reshape((-1,1))
966 Un = numpy.asarray( U ).reshape((-1,1))
970 # Predict (Time Update)
971 # ---------------------
972 if selfA._parameters["EstimationOf"] == "State":
973 if __CovForecast or __LinEvolution:
974 Mt = EM["Tangent"].asMatrix(Xn)
975 Mt = Mt.reshape(Xn.size,Xn.size) # ADAO & check shape
977 Ma = EM["Adjoint"].asMatrix(Xn)
978 Ma = Ma.reshape(Xn.size,Xn.size) # ADAO & check shape
979 Pn_predicted = Q + Mt @ (Pn @ Ma)
981 Xn_predicted = Mt @ Xn
983 M = EM["Direct"].appliedControledFormTo
984 Xn_predicted = M( (Xn, Un) )
985 if CM is not None and "Tangent" in CM and Un is not None: # Attention : si Cm est aussi dans M, doublon !
986 Cm = CM["Tangent"].asMatrix(Xn_predicted)
987 Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape
988 Xn_predicted = Xn_predicted + (Cm @ Un).reshape((-1,1))
989 elif selfA._parameters["EstimationOf"] == "Parameters": # No forecast
990 # --- > Par principe, M = Id, Q = 0
992 if __CovForecast: Pn_predicted = Pn
993 Xn_predicted = numpy.asarray(Xn_predicted).reshape((-1,1))
994 if selfA._toStore("ForecastState"):
995 selfA.StoredVariables["ForecastState"].store( Xn_predicted )
997 if hasattr(Pn_predicted,"asfullmatrix"):
998 Pn_predicted = Pn_predicted.asfullmatrix(Xn.size)
1000 Pn_predicted = numpy.asarray(Pn_predicted).reshape((Xn.size,Xn.size))
1001 if selfA._toStore("ForecastCovariance"):
1002 selfA.StoredVariables["ForecastCovariance"].store( Pn_predicted )
1004 # Correct (Measurement Update)
1005 # ----------------------------
1007 oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, Pn_predicted, True)
1009 oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, B, True)
1011 #--------------------------
1012 Xn = selfA._getInternalState("Xn")
1013 if __CovForecast: Pn = selfA._getInternalState("Pn")
1017 # ==============================================================================
1018 if __name__ == "__main__":
1019 print('\n AUTODIAGNOSTIC\n')