1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2020 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, time, copy, types, sys, logging
29 import math, numpy, scipy
30 from daCore.BasicObjects import Operator
31 from daCore.PlatformInfo import PlatformInfo
32 mpr = PlatformInfo().MachinePrecision()
33 mfp = PlatformInfo().MaximumPrecision()
34 # logging.getLogger().setLevel(logging.DEBUG)
36 # ==============================================================================
37 def ExecuteFunction( paire ):
38 assert len(paire) == 2, "Incorrect number of arguments"
40 __X = numpy.asmatrix(numpy.ravel( X )).T
41 __sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"])
42 __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
43 __fonction = getattr(__module,funcrepr["__userFunction__name"])
44 sys.path = __sys_path_tmp ; del __sys_path_tmp
45 __HX = __fonction( __X )
46 return numpy.ravel( __HX )
48 # ==============================================================================
49 class FDApproximation(object):
51 Cette classe sert d'interface pour définir les opérateurs approximés. A la
52 création d'un objet, en fournissant une fonction "Function", on obtient un
53 objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
54 "AdjointOperator". On contrôle l'approximation DF avec l'incrément
55 multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
56 "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
57 centrées si le booléen "centeredDF" est vrai.
64 avoidingRedundancy = True,
65 toleranceInRedundancy = 1.e-18,
66 lenghtOfRedundancy = -1,
73 import multiprocessing
74 self.__mpEnabled = True
76 self.__mpEnabled = False
78 self.__mpEnabled = False
79 self.__mpWorkers = mpWorkers
80 if self.__mpWorkers is not None and self.__mpWorkers < 1:
81 self.__mpWorkers = None
82 logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
85 self.__mfEnabled = True
87 self.__mfEnabled = False
88 logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
90 if avoidingRedundancy:
92 self.__tolerBP = float(toleranceInRedundancy)
93 self.__lenghtRJ = int(lenghtOfRedundancy)
94 self.__listJPCP = [] # Jacobian Previous Calculated Points
95 self.__listJPCI = [] # Jacobian Previous Calculated Increment
96 self.__listJPCR = [] # Jacobian Previous Calculated Results
97 self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
98 self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
100 self.__avoidRC = False
103 if isinstance(Function,types.FunctionType):
104 logging.debug("FDA Calculs en multiprocessing : FunctionType")
105 self.__userFunction__name = Function.__name__
107 mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
109 mod = os.path.abspath(Function.__globals__['__file__'])
110 if not os.path.isfile(mod):
111 raise ImportError("No user defined function or method found with the name %s"%(mod,))
112 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
113 self.__userFunction__path = os.path.dirname(mod)
115 self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
116 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
117 elif isinstance(Function,types.MethodType):
118 logging.debug("FDA Calculs en multiprocessing : MethodType")
119 self.__userFunction__name = Function.__name__
121 mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
123 mod = os.path.abspath(Function.__func__.__globals__['__file__'])
124 if not os.path.isfile(mod):
125 raise ImportError("No user defined function or method found with the name %s"%(mod,))
126 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
127 self.__userFunction__path = os.path.dirname(mod)
129 self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
130 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
132 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
134 self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
135 self.__userFunction = self.__userOperator.appliedTo
137 self.__centeredDF = bool(centeredDF)
138 if abs(float(increment)) > 1.e-15:
139 self.__increment = float(increment)
141 self.__increment = 0.01
145 self.__dX = numpy.asmatrix(numpy.ravel( dX )).T
146 logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
148 logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
150 # ---------------------------------------------------------
151 def __doublon__(self, e, l, n, v=None):
152 __ac, __iac = False, -1
153 for i in range(len(l)-1,-1,-1):
154 if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
155 __ac, __iac = True, i
156 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
160 # ---------------------------------------------------------
161 def DirectOperator(self, X ):
163 Calcul du direct à l'aide de la fonction fournie.
165 logging.debug("FDA Calcul DirectOperator (explicite)")
167 _HX = self.__userFunction( X, argsAsSerie = True )
169 _X = numpy.asmatrix(numpy.ravel( X )).T
170 _HX = numpy.ravel(self.__userFunction( _X ))
174 # ---------------------------------------------------------
175 def TangentMatrix(self, X ):
177 Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
178 c'est-à-dire le gradient de H en X. On utilise des différences finies
179 directionnelles autour du point X. X est un numpy.matrix.
181 Différences finies centrées (approximation d'ordre 2):
182 1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
183 dX[i] à la composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
184 on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
186 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
188 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
190 Différences finies non centrées (approximation d'ordre 1):
191 1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
192 composante X[i] pour composer X_plus_dXi, et on calcule la réponse
193 HX_plus_dXi = H( X_plus_dXi )
194 2/ On calcule la valeur centrale HX = H(X)
195 3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
197 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
200 logging.debug("FDA Début du calcul de la Jacobienne")
201 logging.debug("FDA Incrément de............: %s*X"%float(self.__increment))
202 logging.debug("FDA Approximation centrée...: %s"%(self.__centeredDF))
204 if X is None or len(X)==0:
205 raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
207 _X = numpy.asmatrix(numpy.ravel( X )).T
209 if self.__dX is None:
210 _dX = self.__increment * _X
212 _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
214 if (_dX == 0.).any():
217 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
219 _dX = numpy.where( _dX == 0., moyenne, _dX )
221 __alreadyCalculated = False
223 __bidon, __alreadyCalculatedP = self.__doublon__(_X, self.__listJPCP, self.__listJPPN, None)
224 __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
225 if __alreadyCalculatedP == __alreadyCalculatedI > -1:
226 __alreadyCalculated, __i = True, __alreadyCalculatedP
227 logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i)
229 if __alreadyCalculated:
230 logging.debug("FDA Calcul Jacobienne (par récupération du doublon %i)"%__i)
231 _Jacobienne = self.__listJPCR[__i]
233 logging.debug("FDA Calcul Jacobienne (explicite)")
234 if self.__centeredDF:
236 if self.__mpEnabled and not self.__mfEnabled:
238 "__userFunction__path" : self.__userFunction__path,
239 "__userFunction__modl" : self.__userFunction__modl,
240 "__userFunction__name" : self.__userFunction__name,
243 for i in range( len(_dX) ):
245 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
246 _X_plus_dXi[i] = _X[i] + _dXi
247 _X_moins_dXi = numpy.array( _X.A1, dtype=float )
248 _X_moins_dXi[i] = _X[i] - _dXi
250 _jobs.append( (_X_plus_dXi, funcrepr) )
251 _jobs.append( (_X_moins_dXi, funcrepr) )
253 import multiprocessing
254 self.__pool = multiprocessing.Pool(self.__mpWorkers)
255 _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
260 for i in range( len(_dX) ):
261 _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
263 elif self.__mfEnabled:
265 for i in range( len(_dX) ):
267 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
268 _X_plus_dXi[i] = _X[i] + _dXi
269 _X_moins_dXi = numpy.array( _X.A1, dtype=float )
270 _X_moins_dXi[i] = _X[i] - _dXi
272 _xserie.append( _X_plus_dXi )
273 _xserie.append( _X_moins_dXi )
275 _HX_plusmoins_dX = self.DirectOperator( _xserie )
278 for i in range( len(_dX) ):
279 _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
283 for i in range( _dX.size ):
285 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
286 _X_plus_dXi[i] = _X[i] + _dXi
287 _X_moins_dXi = numpy.array( _X.A1, dtype=float )
288 _X_moins_dXi[i] = _X[i] - _dXi
290 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
291 _HX_moins_dXi = self.DirectOperator( _X_moins_dXi )
293 _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
297 if self.__mpEnabled and not self.__mfEnabled:
299 "__userFunction__path" : self.__userFunction__path,
300 "__userFunction__modl" : self.__userFunction__modl,
301 "__userFunction__name" : self.__userFunction__name,
304 _jobs.append( (_X.A1, funcrepr) )
305 for i in range( len(_dX) ):
306 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
307 _X_plus_dXi[i] = _X[i] + _dX[i]
309 _jobs.append( (_X_plus_dXi, funcrepr) )
311 import multiprocessing
312 self.__pool = multiprocessing.Pool(self.__mpWorkers)
313 _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
317 _HX = _HX_plus_dX.pop(0)
320 for i in range( len(_dX) ):
321 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
323 elif self.__mfEnabled:
325 _xserie.append( _X.A1 )
326 for i in range( len(_dX) ):
327 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
328 _X_plus_dXi[i] = _X[i] + _dX[i]
330 _xserie.append( _X_plus_dXi )
332 _HX_plus_dX = self.DirectOperator( _xserie )
334 _HX = _HX_plus_dX.pop(0)
337 for i in range( len(_dX) ):
338 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
342 _HX = self.DirectOperator( _X )
343 for i in range( _dX.size ):
345 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
346 _X_plus_dXi[i] = _X[i] + _dXi
348 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
350 _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
353 _Jacobienne = numpy.asmatrix( numpy.vstack( _Jacobienne ) ).T
355 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
356 while len(self.__listJPCP) > self.__lenghtRJ:
357 self.__listJPCP.pop(0)
358 self.__listJPCI.pop(0)
359 self.__listJPCR.pop(0)
360 self.__listJPPN.pop(0)
361 self.__listJPIN.pop(0)
362 self.__listJPCP.append( copy.copy(_X) )
363 self.__listJPCI.append( copy.copy(_dX) )
364 self.__listJPCR.append( copy.copy(_Jacobienne) )
365 self.__listJPPN.append( numpy.linalg.norm(_X) )
366 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
368 logging.debug("FDA Fin du calcul de la Jacobienne")
372 # ---------------------------------------------------------
373 def TangentOperator(self, paire ):
375 Calcul du tangent à l'aide de la Jacobienne.
378 assert len(paire) == 1, "Incorrect lenght of arguments"
380 assert len(_paire) == 2, "Incorrect number of arguments"
382 assert len(paire) == 2, "Incorrect number of arguments"
385 _Jacobienne = self.TangentMatrix( X )
386 if dX is None or len(dX) == 0:
388 # Calcul de la forme matricielle si le second argument est None
389 # -------------------------------------------------------------
390 if self.__mfEnabled: return [_Jacobienne,]
391 else: return _Jacobienne
394 # Calcul de la valeur linéarisée de H en X appliqué à dX
395 # ------------------------------------------------------
396 _dX = numpy.asmatrix(numpy.ravel( dX )).T
397 _HtX = numpy.dot(_Jacobienne, _dX)
398 if self.__mfEnabled: return [_HtX.A1,]
401 # ---------------------------------------------------------
402 def AdjointOperator(self, paire ):
404 Calcul de l'adjoint à l'aide de la Jacobienne.
407 assert len(paire) == 1, "Incorrect lenght of arguments"
409 assert len(_paire) == 2, "Incorrect number of arguments"
411 assert len(paire) == 2, "Incorrect number of arguments"
414 _JacobienneT = self.TangentMatrix( X ).T
415 if Y is None or len(Y) == 0:
417 # Calcul de la forme matricielle si le second argument est None
418 # -------------------------------------------------------------
419 if self.__mfEnabled: return [_JacobienneT,]
420 else: return _JacobienneT
423 # Calcul de la valeur de l'adjoint en X appliqué à Y
424 # --------------------------------------------------
425 _Y = numpy.asmatrix(numpy.ravel( Y )).T
426 _HaY = numpy.dot(_JacobienneT, _Y)
427 if self.__mfEnabled: return [_HaY.A1,]
430 # ==============================================================================
442 Implémentation informatique de l'algorithme MMQR, basée sur la publication :
443 David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
444 Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
447 # Recuperation des donnees et informations initiales
448 # --------------------------------------------------
449 variables = numpy.ravel( x0 )
450 mesures = numpy.ravel( y )
451 increment = sys.float_info[0]
454 quantile = float(quantile)
456 # Calcul des parametres du MM
457 # ---------------------------
458 tn = float(toler) / n
459 e0 = -tn / math.log(tn)
460 epsilon = (e0-tn)/(1+math.log(e0))
462 # Calculs d'initialisation
463 # ------------------------
464 residus = mesures - numpy.ravel( func( variables ) )
465 poids = 1./(epsilon+numpy.abs(residus))
466 veps = 1. - 2. * quantile - residus * poids
467 lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus)
470 # Recherche iterative
471 # -------------------
472 while (increment > toler) and (iteration < maxfun) :
475 Derivees = numpy.array(fprime(variables))
476 Derivees = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
477 DeriveesT = Derivees.transpose()
478 M = numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) )
479 SM = numpy.transpose(numpy.dot( DeriveesT , veps ))
480 step = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0]
482 variables = variables + step
483 if bounds is not None:
484 while( (variables < numpy.ravel(numpy.asmatrix(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asmatrix(bounds)[:,1])).any() ):
486 variables = variables - step
487 residus = mesures - numpy.ravel( func(variables) )
488 surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
490 while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) :
492 variables = variables - step
493 residus = mesures - numpy.ravel( func(variables) )
494 surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
496 increment = lastsurrogate-surrogate
497 poids = 1./(epsilon+numpy.abs(residus))
498 veps = 1. - 2. * quantile - residus * poids
499 lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus)
503 Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] )
505 return variables, Ecart, [n,p,iteration,increment,0]
507 # ==============================================================================
509 def _BackgroundEnsembleGeneration( _bgcenter, _bgcovariance, _nbmembers, _withSVD = True):
510 "Génération d'un ensemble d'ébauche de taille _nbmembers-1"
511 # ~ numpy.random.seed(1234567)
513 raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
515 U, s, V = numpy.linalg.svd(_bgcovariance, full_matrices=False)
516 _nbctl = len(_bgcenter)
517 if _nbmembers > _nbctl:
518 _Z = numpy.concatenate((numpy.dot(
519 numpy.diag(numpy.sqrt(s[:_nbctl])), V[:_nbctl]),
520 numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1-_nbctl)), axis = 0)
522 _Z = numpy.dot(numpy.diag(numpy.sqrt(s[:_nbmembers-1])), V[:_nbmembers-1])
523 _Zca = _CenteredAnomalies(_Z, _nbmembers)
524 BackgroundEnsemble = (_bgcenter + _Zca.T).T
526 if max(abs(_bgcovariance.flatten())) > 0:
527 _nbctl = len(_bgcenter)
528 _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1)
529 _Zca = _CenteredAnomalies(_Z, _nbmembers)
530 BackgroundEnsemble = (_bgcenter + _Zca.T).T
532 BackgroundEnsemble = numpy.tile([_bgcenter],(_nbmembers,1)).T
533 return BackgroundEnsemble
535 def _CenteredAnomalies(Zr, N):
537 Génère une matrice d'anomalies centrées selon les notes manuscrites de MB
538 et conforme au code de PS avec eps = -1
541 Q = numpy.eye(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
542 Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
543 R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
548 def _IEnKF_cycle_Lag_1_SDA_GN(
554 variant = "IEnKF", # IEnKF or IEKF
559 nbPS = 0, # nbPreviousSteps
562 if logging.getLogger().level < logging.WARNING:
563 assert len(E0.shape) == 2, "Ensemble E0 is not well formed: not of shape 2!"
564 assert len(RIdemi.shape) == 2, "R^{-1/2} is not well formed: not of shape 2!"
565 assert variant in ("IEnKF", "IEKF"), "Variant has to be IEnKF or IEKF"
567 nbCtl, nbMbr = E0.shape
570 if logging.getLogger().level < logging.WARNING:
571 assert RIdemi.shape[0] == RIdemi.shape[1] == nbObs, "R^{-1} not of good size: not of size nbObs!"
573 yo = yObs.reshape((nbObs,1))
574 IN = numpy.identity(nbMbr)
575 if variant == "IEnKF":
576 T = numpy.identity(nbMbr)
577 Tinv = numpy.identity(nbMbr)
578 x00 = numpy.mean(E0, axis = 1)
580 Ap0 = numpy.linalg.pinv( Ah0.T.dot(Ah0) )
581 if logging.getLogger().level < logging.WARNING:
582 assert len(Ah0.shape) == 2, "Ensemble A0 is not well formed, of shape 2!"
583 assert Ah0.shape[0] == nbCtl and Ah0.shape[1] == nbMbr, "Ensemble A0 is not well shaped!"
584 assert abs(max(numpy.mean(Ah0, axis = 1))) < nbMbr*mpr, "Ensemble A0 seems not to be centered!"
586 def _convergence_condition(j, dx, JCurr, JPrev):
588 logging.debug("Convergence on maximum number of iterations per cycle, that reach the limit of %i."%iMaximum)
594 _deltaOnJ = abs(JCurr - JPrev) / JPrev
595 if _deltaOnJ <= jTolerance:
596 logging.debug("Convergence on cost decrement tolerance, that is below the threshold of %.1e."%jTolerance)
599 _deltaOnX = numpy.linalg.norm(dx)
600 if _deltaOnX <= sTolerance:
601 logging.debug("Convergence on norm of state correction, that is below the threshold of %.1e."%sTolerance)
602 return True # En correction de l'état
606 St = dict([(k,[]) for k in [
607 "CurrentState", "CurrentEnsemble",
608 "CostFunctionJb", "CostFunctionJo", "CostFunctionJ",
611 j, convergence, JPrev = 1, False, numpy.nan
613 while not convergence:
614 logging.debug("Internal IEnKS step number %i"%j)
615 St["CurrentState"].append( x1.squeeze() )
616 if variant == "IEnKF": # Transform
619 E1 = x1 + epsilonE * Ah0
620 St["CurrentEnsemble"].append( E1 )
621 E2 = numpy.array([Mnnpu(_x) for _x in E1.T]).reshape((nbCtl, nbMbr)) # Evolution 1->2
622 HEL = numpy.array([Hn(_x) for _x in E2.T]).T # Observation à 2
623 yLm = numpy.mean( HEL, axis = 1).reshape((nbObs,1))
625 if variant == "IEnKF":
629 RIdemidy = RIdemi.dot(yo - yLm)
630 xs = RIdemidy / math.sqrt(nbMbr-1)
631 ES = RIdemi.dot(HA2) / math.sqrt(nbMbr-1)
632 G = numpy.linalg.inv(IN + ES.T.dot(ES))
633 xb = G.dot(ES.T.dot(xs))
634 dx = Ah0.dot(xb) + Ah0.dot(G.dot(Ap0.dot(Ah0.T.dot(x00 - x1))))
636 Jb = float(dx.T.dot(dx))
637 Jo = float(RIdemidy.T.dot(RIdemidy))
639 logging.debug("Values for cost functions are: J = %.5e Jo = %.5e Jb = %.5e"%(J,Jo,Jb))
640 St["CostFunctionJb"].append( Jb )
641 St["CostFunctionJo"].append( Jo )
642 St["CostFunctionJ"].append( J )
646 convergence = _convergence_condition(j, dx, J, JPrev)
649 if variant == "IEnKF":
650 T = numpy.real_if_close(scipy.linalg.sqrtm(G))
651 Tinv = numpy.linalg.inv(T)
653 # Stocke le dernier pas
654 x2 = numpy.mean( E2, axis = 1)
655 if variant == "IEKF":
657 A2 = A2.dot(numpy.linalg.cholesky(G)) / epsilonE
659 St["CurrentState"].append( x2.squeeze() )
660 St["CurrentEnsemble"].append( E2 )
662 IndexMin = numpy.argmin( St["CostFunctionJ"][nbPS:] ) + nbPS
663 xa = St["CurrentState"][IndexMin]
664 Ea = St["CurrentEnsemble"][IndexMin]
669 xb = None, # Background (None si E0)
670 E0 = None, # Background ensemble (None si xb)
671 yObs = None, # Observation (série)
673 RIdemi = None, # R^(-1/2)
674 Mnnpu = None, # Evolution operator
675 Hn = None, # Observation operator
676 variant = "IEnKF", # IEnKF or IEKF
677 nMembers = 5, # Number of members
678 sMaximum = 0, # Number of spinup steps
679 cMaximum = 15000, # Number of steps or cycles
680 iMaximum = 15000, # Number of iterations per cycle
681 sTolerance = mfp, # State correction tolerance
682 jTolerance = mfp, # Cost decrement tolerance
685 nbPS = 0, # Number of previous steps
690 if setSeed is not None: numpy.random.seed(setSeed)
691 if E0 is None: E0 = _BackgroundEnsembleGeneration( xb, B, nMembers)
698 xa, Ea, Sa = [xb,], [E0,], [{}]
699 for step in range(cMaximum):
700 if hasattr(yObs,"store"): Ynpu = numpy.ravel( yObs[step+1] )
701 elif type(yObs) in [list, tuple]: Ynpu = numpy.ravel( yObs[step+1] )
702 else: Ynpu = numpy.ravel( yObs )
704 (xa_c, Ea_c, Sa_c) = _IEnKF_cycle_Lag_1_SDA_GN(
721 # Inflation for next cycle
722 E0 = xa_c + inflation * (Ea_c - xa_c)
726 def _IEnKS_cycle_Lag_L_SDA_GN(
732 method = "Transform",
741 if logging.getLogger().level < logging.WARNING:
742 assert len(E0.shape) == 2, "Ensemble E0 is not well formed: not of shape 2!"
743 assert len(RIdemi.shape) == 2, "R^{-1/2} is not well formed: not of shape 2!"
744 assert method in ("Transform", "Bundle"), "Method has to be Transform or Bundle"
746 nbCtl, nbMbr = E0.shape
749 if logging.getLogger().level < logging.WARNING:
750 assert RIdemi.shape[0] == RIdemi.shape[1] == nbObs, "R^{-1} not of good size: not of size nbObs!"
752 yo = yObs.reshape((nbObs,1))
753 IN = numpy.identity(nbMbr)
754 if method == "Transform":
755 T = numpy.identity(nbMbr)
756 Tinv = numpy.identity(nbMbr)
757 x00 = numpy.mean(E0, axis = 1)
759 Am0 = (1/math.sqrt(nbMbr - 1)) * Ah0
760 w = numpy.zeros((nbMbr,1))
761 if logging.getLogger().level < logging.WARNING:
762 assert len(Ah0.shape) == 2, "Ensemble A0 is not well formed, of shape 2!"
763 assert Ah0.shape[0] == nbCtl and Ah0.shape[1] == nbMbr, "Ensemble A0 is not well shaped!"
764 assert abs(max(numpy.mean(Ah0, axis = 1))) < nbMbr*mpr, "Ensemble A0 seems not to be centered!"
766 def _convergence_condition(j, dw, JCurr, JPrev):
768 logging.debug("Convergence on maximum number of iterations per cycle, that reach the limit of %i."%iMaximum)
774 _deltaOnJ = abs(JCurr - JPrev) / JPrev
775 if _deltaOnJ <= jTolerance:
776 logging.debug("Convergence on cost decrement tolerance, that is below the threshold of %.1e."%jTolerance)
779 _deltaOnW = numpy.sqrt(numpy.mean(dw.squeeze()**2))
780 if _deltaOnW <= sTolerance:
781 logging.debug("Convergence on norm of weights correction, that is below the threshold of %.1e."%sTolerance)
782 return True # En correction des poids
786 St = dict([(k,[]) for k in [
787 "CurrentState", "CurrentEnsemble", "CurrentWeights",
788 "CostFunctionJb", "CostFunctionJo", "CostFunctionJ",
791 j, convergence, JPrev = 1, False, numpy.nan
792 while not convergence:
793 logging.debug("Internal IEnKS step number %i"%j)
794 x0 = x00 + Am0.dot( w )
795 St["CurrentState"].append( x0.squeeze() )
796 if method == "Transform":
799 E0 = x0 + epsilon * Am0
800 St["CurrentEnsemble"].append( E0 )
802 yHmean = numpy.mean(E0, axis = 1)
803 for k in range(1, Lag+1):
804 Ek = numpy.array([Mnnpu(_x) for _x in Ek.T]).reshape((nbCtl, nbMbr)) # Evolution 0->L
805 if method == "Transform":
806 yHmean = Mnnpu(yHmean)
807 HEL = numpy.array([Hn(_x) for _x in Ek.T]).T # Observation à L
809 if method == "Transform":
810 yLm = Hn( yHmean ).reshape((nbObs,1))
811 YL = RIdemi.dot( (HEL - numpy.mean( HEL, axis = 1).reshape((nbObs,1))).dot(Tinv) ) / math.sqrt(nbMbr-1)
813 yLm = numpy.mean( HEL, axis = 1).reshape((nbObs,1))
814 YL = RIdemi.dot(HEL - yLm) / epsilon
815 dy = RIdemi.dot(yo - yLm)
817 Jb = float(w.T.dot(w))
818 Jo = float(dy.T.dot(dy))
820 logging.debug("Values for cost functions are: J = %.5e Jo = %.5e Jb = %.5e"%(J,Jo,Jb))
821 St["CurrentWeights"].append( w.squeeze() )
822 St["CostFunctionJb"].append( Jb )
823 St["CostFunctionJo"].append( Jo )
824 St["CostFunctionJ"].append( J )
825 if method == "Transform":
826 GradJ = w - YL.T.dot(dy)
827 HTild = IN + YL.T.dot(YL)
829 GradJ = (nbMbr - 1)*w - YL.T.dot(RIdemi.dot(dy))
830 HTild = (nbMbr - 1)*IN + YL.T.dot(RIdemi.dot(YL))
831 HTild = numpy.array(HTild, dtype=float)
832 dw = numpy.linalg.solve( HTild, numpy.array(GradJ, dtype=float) )
835 convergence = _convergence_condition(j, dw, J, JPrev)
838 if method == "Transform":
839 (U, s, _) = numpy.linalg.svd(HTild, full_matrices=False) # Hess = U s V
840 T = U.dot(numpy.diag(numpy.sqrt(1./s)).dot(U.T)) # T = Hess^(-1/2)
841 Tinv = U.dot(numpy.diag(numpy.sqrt(s)).dot(U.T)) # Tinv = T^(-1)
843 # Stocke le dernier pas
844 St["CurrentState"].append( numpy.mean( Ek, axis = 1).squeeze() )
845 St["CurrentEnsemble"].append( Ek )
847 IndexMin = numpy.argmin( St["CostFunctionJ"][nbPS:] ) + nbPS
848 xa = St["CurrentState"][IndexMin]
849 Ea = St["CurrentEnsemble"][IndexMin]
854 xb = None, # Background
855 yObs = None, # Observation (série)
856 E0 = None, # Background ensemble
858 RIdemi = None, # R^(-1/2)
859 Mnnpu = None, # Evolution operator
860 Hn = None, # Observation operator
861 method = "Transform", # Bundle ou Transform
862 nMembers = 5, # Number of members
863 cMaximum = 15000, # Number of steps or cycles
864 iMaximum = 15000, # Number of iterations per cycle
865 sTolerance = mfp, # Weights correction tolerance
866 jTolerance = mfp, # Cost decrement tolerance
867 Lag = 1, # Lenght of smoothing window
870 nbPS = 0, # Number of previous steps
875 if setSeed is not None: numpy.random.seed(setSeed)
876 if E0 is None: E0 = _BackgroundEnsembleGeneration( xb, B, nMembers)
883 xa, Ea, Sa = [], [], []
884 for i in range(Lag): # Lag void results
888 for i in range(Lag,cMaximum):
889 (xa_c, Ea_c, Sa_c) = _IEnKS_cycle_Lag_L_SDA_GN(
907 # Inflation for next cycle
908 E0 = xa_c + inflation * (Ea_c - xa_c)
912 # ==============================================================================
913 if __name__ == "__main__":
914 print('\n AUTODIAGNOSTIC\n')