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( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
127 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
128 elif isinstance(Function,types.MethodType):
129 logging.debug("FDA Calculs en multiprocessing : MethodType")
130 self.__userFunction__name = Function.__name__
132 mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
134 mod = os.path.abspath(Function.__func__.__globals__['__file__'])
135 if not os.path.isfile(mod):
136 raise ImportError("No user defined function or method found with the name %s"%(mod,))
137 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
138 self.__userFunction__path = os.path.dirname(mod)
140 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
141 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
143 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
145 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
146 self.__userFunction = self.__userOperator.appliedTo
148 self.__centeredDF = bool(centeredDF)
149 if abs(float(increment)) > 1.e-15:
150 self.__increment = float(increment)
152 self.__increment = 0.01
156 self.__dX = numpy.ravel( dX )
158 # ---------------------------------------------------------
159 def __doublon__(self, e, l, n, v=None):
160 __ac, __iac = False, -1
161 for i in range(len(l)-1,-1,-1):
162 if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
163 __ac, __iac = True, i
164 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
168 # ---------------------------------------------------------
169 def __listdotwith__(self, __LMatrix, __dotWith = None, __dotTWith = None):
170 "Produit incrémental d'une matrice liste de colonnes avec un vecteur"
171 if not isinstance(__LMatrix, (list,tuple)):
172 raise TypeError("Columnwise list matrix has not the proper type: %s"%type(__LMatrix))
173 if __dotWith is not None:
174 __Idwx = numpy.ravel( __dotWith )
175 assert len(__LMatrix) == __Idwx.size, "Incorrect size of elements"
176 __Produit = numpy.zeros(__LMatrix[0].size)
177 for i, col in enumerate(__LMatrix):
178 __Produit += float(__Idwx[i]) * col
180 elif __dotTWith is not None:
181 _Idwy = numpy.ravel( __dotTWith ).T
182 assert __LMatrix[0].size == _Idwy.size, "Incorrect size of elements"
183 __Produit = numpy.zeros(len(__LMatrix))
184 for i, col in enumerate(__LMatrix):
185 __Produit[i] = float( _Idwy @ col)
191 # ---------------------------------------------------------
192 def DirectOperator(self, X, **extraArgs ):
194 Calcul du direct à l'aide de la fonction fournie.
196 NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
197 ne doivent pas être données ici à la fonction utilisateur.
199 logging.debug("FDA Calcul DirectOperator (explicite)")
201 _HX = self.__userFunction( X, argsAsSerie = True )
203 _HX = numpy.ravel(self.__userFunction( numpy.ravel(X) ))
207 # ---------------------------------------------------------
208 def TangentMatrix(self, X, dotWith = None, dotTWith = None ):
210 Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
211 c'est-à-dire le gradient de H en X. On utilise des différences finies
212 directionnelles autour du point X. X est un numpy.ndarray.
214 Différences finies centrées (approximation d'ordre 2):
215 1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
216 dX[i] à la composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
217 on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
219 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
221 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
223 Différences finies non centrées (approximation d'ordre 1):
224 1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
225 composante X[i] pour composer X_plus_dXi, et on calcule la réponse
226 HX_plus_dXi = H( X_plus_dXi )
227 2/ On calcule la valeur centrale HX = H(X)
228 3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
230 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
233 logging.debug("FDA Début du calcul de la Jacobienne")
234 logging.debug("FDA Incrément de............: %s*X"%float(self.__increment))
235 logging.debug("FDA Approximation centrée...: %s"%(self.__centeredDF))
237 if X is None or len(X)==0:
238 raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
240 _X = numpy.ravel( X )
242 if self.__dX is None:
243 _dX = self.__increment * _X
245 _dX = numpy.ravel( self.__dX )
246 assert len(_X) == len(_dX), "Inconsistent dX increment length with respect to the X one"
247 assert _X.size == _dX.size, "Inconsistent dX increment size with respect to the X one"
249 if (_dX == 0.).any():
252 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
254 _dX = numpy.where( _dX == 0., moyenne, _dX )
256 __alreadyCalculated = False
258 __bidon, __alreadyCalculatedP = self.__doublon__(_X, self.__listJPCP, self.__listJPPN, None)
259 __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
260 if __alreadyCalculatedP == __alreadyCalculatedI > -1:
261 __alreadyCalculated, __i = True, __alreadyCalculatedP
262 logging.debug("FDA Cas J déjà calculé, récupération du doublon %i"%__i)
264 if __alreadyCalculated:
265 logging.debug("FDA Calcul Jacobienne (par récupération du doublon %i)"%__i)
266 _Jacobienne = self.__listJPCR[__i]
267 logging.debug("FDA Fin du calcul de la Jacobienne")
268 if dotWith is not None:
269 return numpy.dot(_Jacobienne, numpy.ravel( dotWith ))
270 elif dotTWith is not None:
271 return numpy.dot(_Jacobienne.T, numpy.ravel( dotTWith ))
273 logging.debug("FDA Calcul Jacobienne (explicite)")
274 if self.__centeredDF:
276 if self.__mpEnabled and not self.__mfEnabled:
278 "__userFunction__path" : self.__userFunction__path,
279 "__userFunction__modl" : self.__userFunction__modl,
280 "__userFunction__name" : self.__userFunction__name,
283 for i in range( len(_dX) ):
285 _X_plus_dXi = numpy.array( _X, dtype=float )
286 _X_plus_dXi[i] = _X[i] + _dXi
287 _X_moins_dXi = numpy.array( _X, dtype=float )
288 _X_moins_dXi[i] = _X[i] - _dXi
290 _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
291 _jobs.append( (_X_moins_dXi, self.__extraArgs, funcrepr) )
293 import multiprocessing
294 self.__pool = multiprocessing.Pool(self.__mpWorkers)
295 _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
300 for i in range( len(_dX) ):
301 _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
303 elif self.__mfEnabled:
305 for i in range( len(_dX) ):
307 _X_plus_dXi = numpy.array( _X, dtype=float )
308 _X_plus_dXi[i] = _X[i] + _dXi
309 _X_moins_dXi = numpy.array( _X, dtype=float )
310 _X_moins_dXi[i] = _X[i] - _dXi
312 _xserie.append( _X_plus_dXi )
313 _xserie.append( _X_moins_dXi )
315 _HX_plusmoins_dX = self.DirectOperator( _xserie )
318 for i in range( len(_dX) ):
319 _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
323 for i in range( _dX.size ):
325 _X_plus_dXi = numpy.array( _X, dtype=float )
326 _X_plus_dXi[i] = _X[i] + _dXi
327 _X_moins_dXi = numpy.array( _X, dtype=float )
328 _X_moins_dXi[i] = _X[i] - _dXi
330 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
331 _HX_moins_dXi = self.DirectOperator( _X_moins_dXi )
333 _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
337 if self.__mpEnabled and not self.__mfEnabled:
339 "__userFunction__path" : self.__userFunction__path,
340 "__userFunction__modl" : self.__userFunction__modl,
341 "__userFunction__name" : self.__userFunction__name,
344 _jobs.append( (_X, self.__extraArgs, funcrepr) )
345 for i in range( len(_dX) ):
346 _X_plus_dXi = numpy.array( _X, dtype=float )
347 _X_plus_dXi[i] = _X[i] + _dX[i]
349 _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
351 import multiprocessing
352 self.__pool = multiprocessing.Pool(self.__mpWorkers)
353 _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
357 _HX = _HX_plus_dX.pop(0)
360 for i in range( len(_dX) ):
361 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
363 elif self.__mfEnabled:
366 for i in range( len(_dX) ):
367 _X_plus_dXi = numpy.array( _X, dtype=float )
368 _X_plus_dXi[i] = _X[i] + _dX[i]
370 _xserie.append( _X_plus_dXi )
372 _HX_plus_dX = self.DirectOperator( _xserie )
374 _HX = _HX_plus_dX.pop(0)
377 for i in range( len(_dX) ):
378 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
382 _HX = self.DirectOperator( _X )
383 for i in range( _dX.size ):
385 _X_plus_dXi = numpy.array( _X, dtype=float )
386 _X_plus_dXi[i] = _X[i] + _dXi
388 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
390 _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
392 if (dotWith is not None) or (dotTWith is not None):
393 __Produit = self.__listdotwith__(_Jacobienne, dotWith, dotTWith)
396 if __Produit is None or self.__avoidRC:
397 _Jacobienne = numpy.transpose( numpy.vstack( _Jacobienne ) )
399 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
400 while len(self.__listJPCP) > self.__lenghtRJ:
401 self.__listJPCP.pop(0)
402 self.__listJPCI.pop(0)
403 self.__listJPCR.pop(0)
404 self.__listJPPN.pop(0)
405 self.__listJPIN.pop(0)
406 self.__listJPCP.append( copy.copy(_X) )
407 self.__listJPCI.append( copy.copy(_dX) )
408 self.__listJPCR.append( copy.copy(_Jacobienne) )
409 self.__listJPPN.append( numpy.linalg.norm(_X) )
410 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
411 logging.debug("FDA Fin du calcul de la Jacobienne")
412 if __Produit is not None:
417 # ---------------------------------------------------------
418 def TangentOperator(self, paire, **extraArgs ):
420 Calcul du tangent à l'aide de la Jacobienne.
422 NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
423 ne doivent pas être données ici à la fonction utilisateur.
426 assert len(paire) == 1, "Incorrect length of arguments"
428 assert len(_paire) == 2, "Incorrect number of arguments"
430 assert len(paire) == 2, "Incorrect number of arguments"
433 if dX is None or len(dX) == 0:
435 # Calcul de la forme matricielle si le second argument est None
436 # -------------------------------------------------------------
437 _Jacobienne = self.TangentMatrix( X )
438 if self.__mfEnabled: return [_Jacobienne,]
439 else: return _Jacobienne
442 # Calcul de la valeur linéarisée de H en X appliqué à dX
443 # ------------------------------------------------------
444 _HtX = self.TangentMatrix( X, dotWith = dX )
445 if self.__mfEnabled: return [_HtX,]
448 # ---------------------------------------------------------
449 def AdjointOperator(self, paire, **extraArgs ):
451 Calcul de l'adjoint à l'aide de la Jacobienne.
453 NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
454 ne doivent pas être données ici à la fonction utilisateur.
457 assert len(paire) == 1, "Incorrect length of arguments"
459 assert len(_paire) == 2, "Incorrect number of arguments"
461 assert len(paire) == 2, "Incorrect number of arguments"
464 if Y is None or len(Y) == 0:
466 # Calcul de la forme matricielle si le second argument est None
467 # -------------------------------------------------------------
468 _JacobienneT = self.TangentMatrix( X ).T
469 if self.__mfEnabled: return [_JacobienneT,]
470 else: return _JacobienneT
473 # Calcul de la valeur de l'adjoint en X appliqué à Y
474 # --------------------------------------------------
475 _HaY = self.TangentMatrix( X, dotTWith = Y )
476 if self.__mfEnabled: return [_HaY,]
479 # ==============================================================================
480 def EnsembleOfCenteredPerturbations( __bgCenter, __bgCovariance, __nbMembers ):
481 "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
483 __bgCenter = numpy.ravel(__bgCenter)[:,None]
485 raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
487 if __bgCovariance is None:
488 _Perturbations = numpy.tile( __bgCenter, __nbMembers)
490 _Z = numpy.random.multivariate_normal(numpy.zeros(__bgCenter.size), __bgCovariance, size=__nbMembers).T
491 _Perturbations = numpy.tile( __bgCenter, __nbMembers) + _Z
493 return _Perturbations
495 # ==============================================================================
496 def EnsembleOfBackgroundPerturbations( __bgCenter, __bgCovariance, __nbMembers, __withSVD = True):
497 "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
498 def __CenteredRandomAnomalies(Zr, N):
500 Génère une matrice de N anomalies aléatoires centrées sur Zr selon les
501 notes manuscrites de MB et conforme au code de PS avec eps = -1
504 Q = numpy.identity(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
505 Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
506 R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
511 __bgCenter = numpy.ravel(__bgCenter).reshape((-1,1))
513 raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
514 if __bgCovariance is None:
515 _Perturbations = numpy.tile( __bgCenter, __nbMembers)
518 _U, _s, _V = numpy.linalg.svd(__bgCovariance, full_matrices=False)
519 _nbctl = __bgCenter.size
520 if __nbMembers > _nbctl:
521 _Z = numpy.concatenate((numpy.dot(
522 numpy.diag(numpy.sqrt(_s[:_nbctl])), _V[:_nbctl]),
523 numpy.random.multivariate_normal(numpy.zeros(_nbctl),__bgCovariance,__nbMembers-1-_nbctl)), axis = 0)
525 _Z = numpy.dot(numpy.diag(numpy.sqrt(_s[:__nbMembers-1])), _V[:__nbMembers-1])
526 _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
527 _Perturbations = __bgCenter + _Zca
529 if max(abs(__bgCovariance.flatten())) > 0:
530 _nbctl = __bgCenter.size
531 _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),__bgCovariance,__nbMembers-1)
532 _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
533 _Perturbations = __bgCenter + _Zca
535 _Perturbations = numpy.tile( __bgCenter, __nbMembers)
537 return _Perturbations
539 # ==============================================================================
540 def EnsembleMean( __Ensemble ):
541 "Renvoie la moyenne empirique d'un ensemble"
542 return numpy.asarray(__Ensemble).mean(axis=1, dtype=mfp).astype('float').reshape((-1,1))
544 # ==============================================================================
545 def EnsembleOfAnomalies( __Ensemble, __OptMean = None, __Normalisation = 1.):
546 "Renvoie les anomalies centrées à partir d'un ensemble"
547 if __OptMean is None:
548 __Em = EnsembleMean( __Ensemble )
550 __Em = numpy.ravel( __OptMean ).reshape((-1,1))
552 return __Normalisation * (numpy.asarray( __Ensemble ) - __Em)
554 # ==============================================================================
555 def EnsembleErrorCovariance( __Ensemble, __Quick = False ):
556 "Renvoie l'estimation empirique de la covariance d'ensemble"
558 # Covariance rapide mais rarement définie positive
559 __Covariance = numpy.cov( __Ensemble )
561 # Résultat souvent identique à numpy.cov, mais plus robuste
562 __n, __m = numpy.asarray( __Ensemble ).shape
563 __Anomalies = EnsembleOfAnomalies( __Ensemble )
564 # Estimation empirique
565 __Covariance = ( __Anomalies @ __Anomalies.T ) / (__m-1)
567 __Covariance = ( __Covariance + __Covariance.T ) * 0.5
568 # Assure la positivité
569 __epsilon = mpr*numpy.trace( __Covariance )
570 __Covariance = __Covariance + __epsilon * numpy.identity(__n)
574 # ==============================================================================
575 def EnsemblePerturbationWithGivenCovariance(
580 "Ajout d'une perturbation à chaque membre d'un ensemble selon une covariance prescrite"
581 if hasattr(__Covariance,"assparsematrix"):
582 if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance.assparsematrix())/abs(__Ensemble).mean() < mpr).all():
583 # Traitement d'une covariance nulle ou presque
585 if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance.assparsematrix()) < mpr).all():
586 # Traitement d'une covariance nulle ou presque
589 if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance)/abs(__Ensemble).mean() < mpr).all():
590 # Traitement d'une covariance nulle ou presque
592 if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance) < mpr).all():
593 # Traitement d'une covariance nulle ou presque
596 __n, __m = __Ensemble.shape
597 if __Seed is not None: numpy.random.seed(__Seed)
599 if hasattr(__Covariance,"isscalar") and __Covariance.isscalar():
600 # Traitement d'une covariance multiple de l'identité
602 __std = numpy.sqrt(__Covariance.assparsematrix())
603 __Ensemble += numpy.random.normal(__zero, __std, size=(__m,__n)).T
605 elif hasattr(__Covariance,"isvector") and __Covariance.isvector():
606 # Traitement d'une covariance diagonale avec variances non identiques
607 __zero = numpy.zeros(__n)
608 __std = numpy.sqrt(__Covariance.assparsematrix())
609 __Ensemble += numpy.asarray([numpy.random.normal(__zero, __std) for i in range(__m)]).T
611 elif hasattr(__Covariance,"ismatrix") and __Covariance.ismatrix():
612 # Traitement d'une covariance pleine
613 __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance.asfullmatrix(__n), size=__m).T
615 elif isinstance(__Covariance, numpy.ndarray):
616 # Traitement d'une covariance numpy pleine, sachant qu'on arrive ici en dernier
617 __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance, size=__m).T
620 raise ValueError("Error in ensemble perturbation with inadequate covariance specification")
624 # ==============================================================================
625 def CovarianceInflation(
627 __InflationType = None,
628 __InflationFactor = None,
629 __BackgroundCov = None,
632 Inflation applicable soit sur Pb ou Pa, soit sur les ensembles EXb ou EXa
634 Synthèse : Hunt 2007, section 2.3.5
636 if __InflationFactor is None:
637 return __InputCovOrEns
639 __InflationFactor = float(__InflationFactor)
641 __InputCovOrEns = numpy.asarray(__InputCovOrEns)
642 if __InputCovOrEns.size == 0: return __InputCovOrEns
644 if __InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
645 if __InflationFactor < 1.:
646 raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
647 if __InflationFactor < 1.+mpr: # No inflation = 1
648 return __InputCovOrEns
649 __OutputCovOrEns = __InflationFactor**2 * __InputCovOrEns
651 elif __InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
652 if __InflationFactor < 1.:
653 raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
654 if __InflationFactor < 1.+mpr: # No inflation = 1
655 return __InputCovOrEns
656 __InputCovOrEnsMean = __InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
657 __OutputCovOrEns = __InputCovOrEnsMean[:,numpy.newaxis] \
658 + __InflationFactor * (__InputCovOrEns - __InputCovOrEnsMean[:,numpy.newaxis])
660 elif __InflationType in ["AdditiveOnAnalysisCovariance", "AdditiveOnBackgroundCovariance"]:
661 if __InflationFactor < 0.:
662 raise ValueError("Inflation factor for additive inflation has to be greater or equal than 0.")
663 if __InflationFactor < mpr: # No inflation = 0
664 return __InputCovOrEns
665 __n, __m = __InputCovOrEns.shape
667 raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
668 __tr = __InputCovOrEns.trace()/__n
669 if __InflationFactor > __tr:
670 raise ValueError("Inflation factor for additive inflation has to be small over %.0e."%__tr)
671 __OutputCovOrEns = (1. - __InflationFactor)*__InputCovOrEns + __InflationFactor * numpy.identity(__n)
673 elif __InflationType == "HybridOnBackgroundCovariance":
674 if __InflationFactor < 0.:
675 raise ValueError("Inflation factor for hybrid inflation has to be greater or equal than 0.")
676 if __InflationFactor < mpr: # No inflation = 0
677 return __InputCovOrEns
678 __n, __m = __InputCovOrEns.shape
680 raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
681 if __BackgroundCov is None:
682 raise ValueError("Background covariance matrix B has to be given for hybrid inflation.")
683 if __InputCovOrEns.shape != __BackgroundCov.shape:
684 raise ValueError("Ensemble covariance matrix has to be of same size than background covariance matrix B.")
685 __OutputCovOrEns = (1. - __InflationFactor) * __InputCovOrEns + __InflationFactor * __BackgroundCov
687 elif __InflationType == "Relaxation":
688 raise NotImplementedError("Relaxation inflation type not implemented")
691 raise ValueError("Error in inflation type, '%s' is not a valid keyword."%__InflationType)
693 return __OutputCovOrEns
695 # ==============================================================================
696 def HessienneEstimation(__selfA, __nb, __HaM, __HtM, __BI, __RI):
697 "Estimation de la Hessienne"
700 for i in range(int(__nb)):
701 __ee = numpy.zeros((__nb,1))
703 __HtEE = numpy.dot(__HtM,__ee).reshape((-1,1))
704 __HessienneI.append( numpy.ravel( __BI * __ee + __HaM * (__RI * __HtEE) ) )
706 __A = numpy.linalg.inv(numpy.array( __HessienneI ))
707 __A = (__A + __A.T) * 0.5 # Symétrie
708 __A = __A + mpr*numpy.trace( __A ) * numpy.identity(__nb) # Positivité
710 if min(__A.shape) != max(__A.shape):
711 raise ValueError("The %s a posteriori covariance matrix A is of shape %s, despites it has to be a squared matrix. There is an error in the observation operator, please check it."%(__selfA._name,str(__A.shape)))
712 if (numpy.diag(__A) < 0).any():
713 raise ValueError("The %s a posteriori covariance matrix A has at least one negative value on its diagonal. There is an error in the observation operator, please check it."%(__selfA._name,))
714 if logging.getLogger().level < logging.WARNING: # La vérification n'a lieu qu'en debug
716 numpy.linalg.cholesky( __A )
718 raise ValueError("The %s a posteriori covariance matrix A is not symmetric positive-definite. Please check your a priori covariances and your observation operator."%(__selfA._name,))
722 # ==============================================================================
723 def QuantilesEstimations(selfA, A, Xa, HXa = None, Hm = None, HtM = None):
724 "Estimation des quantiles a posteriori à partir de A>0 (selfA est modifié)"
725 nbsamples = selfA._parameters["NumberOfSamplesForQuantiles"]
727 # Traitement des bornes
728 if "StateBoundsForQuantiles" in selfA._parameters:
729 LBounds = selfA._parameters["StateBoundsForQuantiles"] # Prioritaire
730 elif "Bounds" in selfA._parameters:
731 LBounds = selfA._parameters["Bounds"] # Défaut raisonnable
734 if LBounds is not None:
735 LBounds = ForceNumericBounds( LBounds )
736 __Xa = numpy.ravel(Xa)
738 # Échantillonnage des états
741 for i in range(nbsamples):
742 if selfA._parameters["SimulationForQuantiles"] == "Linear" and HtM is not None and HXa is not None:
743 dXr = (numpy.random.multivariate_normal(__Xa,A) - __Xa).reshape((-1,1))
744 if LBounds is not None: # "EstimateProjection" par défaut
745 dXr = numpy.max(numpy.hstack((dXr,LBounds[:,0].reshape((-1,1))) - __Xa.reshape((-1,1))),axis=1)
746 dXr = numpy.min(numpy.hstack((dXr,LBounds[:,1].reshape((-1,1))) - __Xa.reshape((-1,1))),axis=1)
748 Yr = HXa.reshape((-1,1)) + dYr
749 if selfA._toStore("SampledStateForQuantiles"): Xr = __Xa + numpy.ravel(dXr)
750 elif selfA._parameters["SimulationForQuantiles"] == "NonLinear" and Hm is not None:
751 Xr = numpy.random.multivariate_normal(__Xa,A)
752 if LBounds is not None: # "EstimateProjection" par défaut
753 Xr = numpy.max(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,0].reshape((-1,1)))),axis=1)
754 Xr = numpy.min(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,1].reshape((-1,1)))),axis=1)
755 Yr = numpy.asarray(Hm( Xr ))
757 raise ValueError("Quantile simulations has only to be Linear or NonLinear.")
760 YfQ = Yr.reshape((-1,1))
761 if selfA._toStore("SampledStateForQuantiles"): EXr = Xr.reshape((-1,1))
763 YfQ = numpy.hstack((YfQ,Yr.reshape((-1,1))))
764 if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.hstack((EXr,Xr.reshape((-1,1))))
766 # Extraction des quantiles
769 for quantile in selfA._parameters["Quantiles"]:
770 if not (0. <= float(quantile) <= 1.): continue
771 indice = int(nbsamples * float(quantile) - 1./nbsamples)
772 if YQ is None: YQ = YfQ[:,indice].reshape((-1,1))
773 else: YQ = numpy.hstack((YQ,YfQ[:,indice].reshape((-1,1))))
774 if YQ is not None: # Liste non vide de quantiles
775 selfA.StoredVariables["SimulationQuantiles"].store( YQ )
776 if selfA._toStore("SampledStateForQuantiles"):
777 selfA.StoredVariables["SampledStateForQuantiles"].store( EXr )
781 # ==============================================================================
782 def ForceNumericBounds( __Bounds ):
783 "Force les bornes à être des valeurs numériques, sauf si globalement None"
784 # Conserve une valeur par défaut à None s'il n'y a pas de bornes
785 if __Bounds is None: return None
786 # Converti toutes les bornes individuelles None à +/- l'infini
787 __Bounds = numpy.asarray( __Bounds, dtype=float )
788 if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
789 raise ValueError("Incorrectly shaped bounds data")
790 __Bounds[numpy.isnan(__Bounds[:,0]),0] = -sys.float_info.max
791 __Bounds[numpy.isnan(__Bounds[:,1]),1] = sys.float_info.max
794 # ==============================================================================
795 def RecentredBounds( __Bounds, __Center, __Scale = None):
796 "Recentre les bornes autour de 0, sauf si globalement None"
797 # Conserve une valeur par défaut à None s'il n'y a pas de bornes
798 if __Bounds is None: return None
800 # Recentre les valeurs numériques de bornes
801 return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1,1))
803 # Recentre les valeurs numériques de bornes et change l'échelle par une matrice
804 return __Scale @ (ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1,1)))
806 # ==============================================================================
807 def ApplyBounds( __Vector, __Bounds, __newClip = True):
808 "Applique des bornes numériques à un point"
809 # Conserve une valeur par défaut s'il n'y a pas de bornes
810 if __Bounds is None: return __Vector
812 if not isinstance(__Vector, numpy.ndarray): # Is an array
813 raise ValueError("Incorrect array definition of vector data")
814 if not isinstance(__Bounds, numpy.ndarray): # Is an array
815 raise ValueError("Incorrect array definition of bounds data")
816 if 2*__Vector.size != __Bounds.size: # Is a 2 column array of vector lenght
817 raise ValueError("Incorrect bounds number (%i) to be applied for this vector (of size %i)"%(__Bounds.size,__Vector.size))
818 if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
819 raise ValueError("Incorrectly shaped bounds data")
822 __Vector = __Vector.clip(
823 __Bounds[:,0].reshape(__Vector.shape),
824 __Bounds[:,1].reshape(__Vector.shape),
827 __Vector = numpy.max(numpy.hstack((__Vector.reshape((-1,1)),numpy.asmatrix(__Bounds)[:,0])),axis=1)
828 __Vector = numpy.min(numpy.hstack((__Vector.reshape((-1,1)),numpy.asmatrix(__Bounds)[:,1])),axis=1)
829 __Vector = numpy.asarray(__Vector)
833 # ==============================================================================
834 def Apply3DVarRecentringOnEnsemble(__EnXn, __EnXf, __Ynpu, __HO, __R, __B, __Betaf):
835 "Recentre l'ensemble Xn autour de l'analyse 3DVAR"
837 Xf = EnsembleMean( __EnXf )
838 Pf = Covariance( asCovariance=EnsembleErrorCovariance(__EnXf) )
839 Pf = (1 - __Betaf) * __B.asfullmatrix(Xf.size) + __Betaf * Pf
841 selfB = PartialAlgorithm("3DVAR")
842 selfB._parameters["Minimizer"] = "LBFGSB"
843 selfB._parameters["MaximumNumberOfSteps"] = 15000
844 selfB._parameters["CostDecrementTolerance"] = 1.e-7
845 selfB._parameters["ProjectedGradientTolerance"] = -1
846 selfB._parameters["GradientNormTolerance"] = 1.e-05
847 selfB._parameters["StoreInternalVariables"] = False
848 selfB._parameters["optiprint"] = -1
849 selfB._parameters["optdisp"] = 0
850 selfB._parameters["Bounds"] = None
851 selfB._parameters["InitializationPoint"] = Xf
852 from daAlgorithms.Atoms import std3dvar
853 std3dvar.std3dvar(selfB, Xf, __Ynpu, None, __HO, None, __R, Pf)
854 Xa = selfB.get("Analysis")[-1].reshape((-1,1))
857 return Xa + EnsembleOfAnomalies( __EnXn )
859 # ==============================================================================
860 def multiXOsteps(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle,
861 __CovForecast = False, __LinEvolution = False,
864 Prévision multi-pas avec une correction par pas (multi-méthodes)
869 if selfA._parameters["EstimationOf"] == "State":
870 if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
871 Xn = numpy.asarray(Xb)
872 if __CovForecast: Pn = B
873 selfA.StoredVariables["Analysis"].store( Xn )
874 if selfA._toStore("APosterioriCovariance"):
875 if hasattr(B,"asfullmatrix"):
876 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(Xn.size) )
878 selfA.StoredVariables["APosterioriCovariance"].store( B )
879 selfA._setInternalState("seed", numpy.random.get_state())
880 elif selfA._parameters["nextStep"]:
881 Xn = selfA._getInternalState("Xn")
882 if __CovForecast: Pn = selfA._getInternalState("Pn")
884 Xn = numpy.asarray(Xb)
885 if __CovForecast: Pn = B
887 if hasattr(Y,"stepnumber"):
888 duration = Y.stepnumber()
894 for step in range(duration-1):
895 selfA.StoredVariables["CurrentStepNumber"].store( len(selfA.StoredVariables["Analysis"]) )
897 if hasattr(Y,"store"):
898 Ynpu = numpy.asarray( Y[step+1] ).reshape((-1,1))
900 Ynpu = numpy.asarray( Y ).reshape((-1,1))
903 if hasattr(U,"store") and len(U)>1:
904 Un = numpy.asarray( U[step] ).reshape((-1,1))
905 elif hasattr(U,"store") and len(U)==1:
906 Un = numpy.asarray( U[0] ).reshape((-1,1))
908 Un = numpy.asarray( U ).reshape((-1,1))
912 # Predict (Time Update)
913 # ---------------------
914 if selfA._parameters["EstimationOf"] == "State":
915 if __CovForecast or __LinEvolution:
916 Mt = EM["Tangent"].asMatrix(Xn)
917 Mt = Mt.reshape(Xn.size,Xn.size) # ADAO & check shape
919 Ma = EM["Adjoint"].asMatrix(Xn)
920 Ma = Ma.reshape(Xn.size,Xn.size) # ADAO & check shape
921 Pn_predicted = Q + Mt @ (Pn @ Ma)
923 Xn_predicted = Mt @ Xn
925 M = EM["Direct"].appliedControledFormTo
926 Xn_predicted = M( (Xn, Un) )
927 if CM is not None and "Tangent" in CM and Un is not None: # Attention : si Cm est aussi dans M, doublon !
928 Cm = CM["Tangent"].asMatrix(Xn_predicted)
929 Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape
930 Xn_predicted = Xn_predicted + (Cm @ Un).reshape((-1,1))
931 elif selfA._parameters["EstimationOf"] == "Parameters": # No forecast
932 # --- > Par principe, M = Id, Q = 0
934 if __CovForecast: Pn_predicted = Pn
935 Xn_predicted = numpy.asarray(Xn_predicted).reshape((-1,1))
936 if selfA._toStore("ForecastState"):
937 selfA.StoredVariables["ForecastState"].store( Xn_predicted )
939 if hasattr(Pn_predicted,"asfullmatrix"):
940 Pn_predicted = Pn_predicted.asfullmatrix(Xn.size)
942 Pn_predicted = numpy.asarray(Pn_predicted).reshape((Xn.size,Xn.size))
943 if selfA._toStore("ForecastCovariance"):
944 selfA.StoredVariables["ForecastCovariance"].store( Pn_predicted )
946 # Correct (Measurement Update)
947 # ----------------------------
949 oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, Pn_predicted, True)
951 oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, B, True)
953 #--------------------------
954 Xn = selfA._getInternalState("Xn")
955 if __CovForecast: Pn = selfA._getInternalState("Pn")
959 # ==============================================================================
960 if __name__ == "__main__":
961 print('\n AUTODIAGNOSTIC\n')