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.
60 name = "FDApproximation",
65 avoidingRedundancy = True,
66 toleranceInRedundancy = 1.e-18,
67 lenghtOfRedundancy = -1,
72 self.__name = str(name)
75 import multiprocessing
76 self.__mpEnabled = True
78 self.__mpEnabled = False
80 self.__mpEnabled = False
81 self.__mpWorkers = mpWorkers
82 if self.__mpWorkers is not None and self.__mpWorkers < 1:
83 self.__mpWorkers = None
84 logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
87 self.__mfEnabled = True
89 self.__mfEnabled = False
90 logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
92 if avoidingRedundancy:
94 self.__tolerBP = float(toleranceInRedundancy)
95 self.__lenghtRJ = int(lenghtOfRedundancy)
96 self.__listJPCP = [] # Jacobian Previous Calculated Points
97 self.__listJPCI = [] # Jacobian Previous Calculated Increment
98 self.__listJPCR = [] # Jacobian Previous Calculated Results
99 self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
100 self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
102 self.__avoidRC = False
105 if isinstance(Function,types.FunctionType):
106 logging.debug("FDA Calculs en multiprocessing : FunctionType")
107 self.__userFunction__name = Function.__name__
109 mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
111 mod = os.path.abspath(Function.__globals__['__file__'])
112 if not os.path.isfile(mod):
113 raise ImportError("No user defined function or method found with the name %s"%(mod,))
114 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
115 self.__userFunction__path = os.path.dirname(mod)
117 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
118 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
119 elif isinstance(Function,types.MethodType):
120 logging.debug("FDA Calculs en multiprocessing : MethodType")
121 self.__userFunction__name = Function.__name__
123 mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
125 mod = os.path.abspath(Function.__func__.__globals__['__file__'])
126 if not os.path.isfile(mod):
127 raise ImportError("No user defined function or method found with the name %s"%(mod,))
128 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
129 self.__userFunction__path = os.path.dirname(mod)
131 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
132 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
134 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
136 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
137 self.__userFunction = self.__userOperator.appliedTo
139 self.__centeredDF = bool(centeredDF)
140 if abs(float(increment)) > 1.e-15:
141 self.__increment = float(increment)
143 self.__increment = 0.01
147 self.__dX = numpy.asmatrix(numpy.ravel( dX )).T
148 logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
150 logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
152 # ---------------------------------------------------------
153 def __doublon__(self, e, l, n, v=None):
154 __ac, __iac = False, -1
155 for i in range(len(l)-1,-1,-1):
156 if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
157 __ac, __iac = True, i
158 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
162 # ---------------------------------------------------------
163 def DirectOperator(self, X ):
165 Calcul du direct à l'aide de la fonction fournie.
167 logging.debug("FDA Calcul DirectOperator (explicite)")
169 _HX = self.__userFunction( X, argsAsSerie = True )
171 _X = numpy.asmatrix(numpy.ravel( X )).T
172 _HX = numpy.ravel(self.__userFunction( _X ))
176 # ---------------------------------------------------------
177 def TangentMatrix(self, X ):
179 Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
180 c'est-à-dire le gradient de H en X. On utilise des différences finies
181 directionnelles autour du point X. X est un numpy.matrix.
183 Différences finies centrées (approximation d'ordre 2):
184 1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
185 dX[i] à la composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
186 on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
188 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
190 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
192 Différences finies non centrées (approximation d'ordre 1):
193 1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
194 composante X[i] pour composer X_plus_dXi, et on calcule la réponse
195 HX_plus_dXi = H( X_plus_dXi )
196 2/ On calcule la valeur centrale HX = H(X)
197 3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
199 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
202 logging.debug("FDA Début du calcul de la Jacobienne")
203 logging.debug("FDA Incrément de............: %s*X"%float(self.__increment))
204 logging.debug("FDA Approximation centrée...: %s"%(self.__centeredDF))
206 if X is None or len(X)==0:
207 raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
209 _X = numpy.asmatrix(numpy.ravel( X )).T
211 if self.__dX is None:
212 _dX = self.__increment * _X
214 _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
216 if (_dX == 0.).any():
219 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
221 _dX = numpy.where( _dX == 0., moyenne, _dX )
223 __alreadyCalculated = False
225 __bidon, __alreadyCalculatedP = self.__doublon__(_X, self.__listJPCP, self.__listJPPN, None)
226 __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
227 if __alreadyCalculatedP == __alreadyCalculatedI > -1:
228 __alreadyCalculated, __i = True, __alreadyCalculatedP
229 logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i)
231 if __alreadyCalculated:
232 logging.debug("FDA Calcul Jacobienne (par récupération du doublon %i)"%__i)
233 _Jacobienne = self.__listJPCR[__i]
235 logging.debug("FDA Calcul Jacobienne (explicite)")
236 if self.__centeredDF:
238 if self.__mpEnabled and not self.__mfEnabled:
240 "__userFunction__path" : self.__userFunction__path,
241 "__userFunction__modl" : self.__userFunction__modl,
242 "__userFunction__name" : self.__userFunction__name,
245 for i in range( len(_dX) ):
247 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
248 _X_plus_dXi[i] = _X[i] + _dXi
249 _X_moins_dXi = numpy.array( _X.A1, dtype=float )
250 _X_moins_dXi[i] = _X[i] - _dXi
252 _jobs.append( (_X_plus_dXi, funcrepr) )
253 _jobs.append( (_X_moins_dXi, funcrepr) )
255 import multiprocessing
256 self.__pool = multiprocessing.Pool(self.__mpWorkers)
257 _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
262 for i in range( len(_dX) ):
263 _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
265 elif self.__mfEnabled:
267 for i in range( len(_dX) ):
269 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
270 _X_plus_dXi[i] = _X[i] + _dXi
271 _X_moins_dXi = numpy.array( _X.A1, dtype=float )
272 _X_moins_dXi[i] = _X[i] - _dXi
274 _xserie.append( _X_plus_dXi )
275 _xserie.append( _X_moins_dXi )
277 _HX_plusmoins_dX = self.DirectOperator( _xserie )
280 for i in range( len(_dX) ):
281 _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
285 for i in range( _dX.size ):
287 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
288 _X_plus_dXi[i] = _X[i] + _dXi
289 _X_moins_dXi = numpy.array( _X.A1, dtype=float )
290 _X_moins_dXi[i] = _X[i] - _dXi
292 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
293 _HX_moins_dXi = self.DirectOperator( _X_moins_dXi )
295 _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
299 if self.__mpEnabled and not self.__mfEnabled:
301 "__userFunction__path" : self.__userFunction__path,
302 "__userFunction__modl" : self.__userFunction__modl,
303 "__userFunction__name" : self.__userFunction__name,
306 _jobs.append( (_X.A1, funcrepr) )
307 for i in range( len(_dX) ):
308 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
309 _X_plus_dXi[i] = _X[i] + _dX[i]
311 _jobs.append( (_X_plus_dXi, funcrepr) )
313 import multiprocessing
314 self.__pool = multiprocessing.Pool(self.__mpWorkers)
315 _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
319 _HX = _HX_plus_dX.pop(0)
322 for i in range( len(_dX) ):
323 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
325 elif self.__mfEnabled:
327 _xserie.append( _X.A1 )
328 for i in range( len(_dX) ):
329 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
330 _X_plus_dXi[i] = _X[i] + _dX[i]
332 _xserie.append( _X_plus_dXi )
334 _HX_plus_dX = self.DirectOperator( _xserie )
336 _HX = _HX_plus_dX.pop(0)
339 for i in range( len(_dX) ):
340 _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
344 _HX = self.DirectOperator( _X )
345 for i in range( _dX.size ):
347 _X_plus_dXi = numpy.array( _X.A1, dtype=float )
348 _X_plus_dXi[i] = _X[i] + _dXi
350 _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
352 _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
355 _Jacobienne = numpy.asmatrix( numpy.vstack( _Jacobienne ) ).T
357 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
358 while len(self.__listJPCP) > self.__lenghtRJ:
359 self.__listJPCP.pop(0)
360 self.__listJPCI.pop(0)
361 self.__listJPCR.pop(0)
362 self.__listJPPN.pop(0)
363 self.__listJPIN.pop(0)
364 self.__listJPCP.append( copy.copy(_X) )
365 self.__listJPCI.append( copy.copy(_dX) )
366 self.__listJPCR.append( copy.copy(_Jacobienne) )
367 self.__listJPPN.append( numpy.linalg.norm(_X) )
368 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
370 logging.debug("FDA Fin du calcul de la Jacobienne")
374 # ---------------------------------------------------------
375 def TangentOperator(self, paire ):
377 Calcul du tangent à l'aide de la Jacobienne.
380 assert len(paire) == 1, "Incorrect lenght of arguments"
382 assert len(_paire) == 2, "Incorrect number of arguments"
384 assert len(paire) == 2, "Incorrect number of arguments"
387 _Jacobienne = self.TangentMatrix( X )
388 if dX is None or len(dX) == 0:
390 # Calcul de la forme matricielle si le second argument est None
391 # -------------------------------------------------------------
392 if self.__mfEnabled: return [_Jacobienne,]
393 else: return _Jacobienne
396 # Calcul de la valeur linéarisée de H en X appliqué à dX
397 # ------------------------------------------------------
398 _dX = numpy.asmatrix(numpy.ravel( dX )).T
399 _HtX = numpy.dot(_Jacobienne, _dX)
400 if self.__mfEnabled: return [_HtX.A1,]
403 # ---------------------------------------------------------
404 def AdjointOperator(self, paire ):
406 Calcul de l'adjoint à l'aide de la Jacobienne.
409 assert len(paire) == 1, "Incorrect lenght of arguments"
411 assert len(_paire) == 2, "Incorrect number of arguments"
413 assert len(paire) == 2, "Incorrect number of arguments"
416 _JacobienneT = self.TangentMatrix( X ).T
417 if Y is None or len(Y) == 0:
419 # Calcul de la forme matricielle si le second argument est None
420 # -------------------------------------------------------------
421 if self.__mfEnabled: return [_JacobienneT,]
422 else: return _JacobienneT
425 # Calcul de la valeur de l'adjoint en X appliqué à Y
426 # --------------------------------------------------
427 _Y = numpy.asmatrix(numpy.ravel( Y )).T
428 _HaY = numpy.dot(_JacobienneT, _Y)
429 if self.__mfEnabled: return [_HaY.A1,]
432 # ==============================================================================
444 Implémentation informatique de l'algorithme MMQR, basée sur la publication :
445 David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
446 Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
449 # Recuperation des donnees et informations initiales
450 # --------------------------------------------------
451 variables = numpy.ravel( x0 )
452 mesures = numpy.ravel( y )
453 increment = sys.float_info[0]
456 quantile = float(quantile)
458 # Calcul des parametres du MM
459 # ---------------------------
460 tn = float(toler) / n
461 e0 = -tn / math.log(tn)
462 epsilon = (e0-tn)/(1+math.log(e0))
464 # Calculs d'initialisation
465 # ------------------------
466 residus = mesures - numpy.ravel( func( variables ) )
467 poids = 1./(epsilon+numpy.abs(residus))
468 veps = 1. - 2. * quantile - residus * poids
469 lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus)
472 # Recherche iterative
473 # -------------------
474 while (increment > toler) and (iteration < maxfun) :
477 Derivees = numpy.array(fprime(variables))
478 Derivees = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
479 DeriveesT = Derivees.transpose()
480 M = numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) )
481 SM = numpy.transpose(numpy.dot( DeriveesT , veps ))
482 step = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0]
484 variables = variables + step
485 if bounds is not None:
486 while( (variables < numpy.ravel(numpy.asmatrix(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asmatrix(bounds)[:,1])).any() ):
488 variables = variables - step
489 residus = mesures - numpy.ravel( func(variables) )
490 surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
492 while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) :
494 variables = variables - step
495 residus = mesures - numpy.ravel( func(variables) )
496 surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
498 increment = lastsurrogate-surrogate
499 poids = 1./(epsilon+numpy.abs(residus))
500 veps = 1. - 2. * quantile - residus * poids
501 lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus)
505 Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] )
507 return variables, Ecart, [n,p,iteration,increment,0]
509 # ==============================================================================
511 def _BackgroundEnsembleGeneration( _bgcenter, _bgcovariance, _nbmembers, _withSVD = True):
512 "Génération d'un ensemble d'ébauche de taille _nbmembers-1"
513 # ~ numpy.random.seed(1234567)
515 raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
517 U, s, V = numpy.linalg.svd(_bgcovariance, full_matrices=False)
518 _nbctl = len(_bgcenter)
519 if _nbmembers > _nbctl:
520 _Z = numpy.concatenate((numpy.dot(
521 numpy.diag(numpy.sqrt(s[:_nbctl])), V[:_nbctl]),
522 numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1-_nbctl)), axis = 0)
524 _Z = numpy.dot(numpy.diag(numpy.sqrt(s[:_nbmembers-1])), V[:_nbmembers-1])
525 _Zca = _CenteredAnomalies(_Z, _nbmembers)
526 BackgroundEnsemble = (_bgcenter + _Zca.T).T
528 if max(abs(_bgcovariance.flatten())) > 0:
529 _nbctl = len(_bgcenter)
530 _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1)
531 _Zca = _CenteredAnomalies(_Z, _nbmembers)
532 BackgroundEnsemble = (_bgcenter + _Zca.T).T
534 BackgroundEnsemble = numpy.tile([_bgcenter],(_nbmembers,1)).T
535 return BackgroundEnsemble
537 def _CenteredAnomalies(Zr, N):
539 Génère une matrice d'anomalies centrées selon les notes manuscrites de MB
540 et conforme au code de PS avec eps = -1
543 Q = numpy.eye(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
544 Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
545 R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
550 def _IEnKF_cycle_Lag_1_SDA_GN(
556 variant = "IEnKF", # IEnKF or IEKF
561 nbPS = 0, # nbPreviousSteps
564 if logging.getLogger().level < logging.WARNING:
565 assert len(E0.shape) == 2, "Ensemble E0 is not well formed: not of shape 2!"
566 assert len(RIdemi.shape) == 2, "R^{-1/2} is not well formed: not of shape 2!"
567 assert variant in ("IEnKF", "IEKF"), "Variant has to be IEnKF or IEKF"
569 nbCtl, nbMbr = E0.shape
572 if logging.getLogger().level < logging.WARNING:
573 assert RIdemi.shape[0] == RIdemi.shape[1] == nbObs, "R^{-1} not of good size: not of size nbObs!"
575 yo = yObs.reshape((nbObs,1))
576 IN = numpy.identity(nbMbr)
577 if variant == "IEnKF":
578 T = numpy.identity(nbMbr)
579 Tinv = numpy.identity(nbMbr)
580 x00 = numpy.mean(E0, axis = 1)
582 Ap0 = numpy.linalg.pinv( Ah0.T.dot(Ah0) )
583 if logging.getLogger().level < logging.WARNING:
584 assert len(Ah0.shape) == 2, "Ensemble A0 is not well formed, of shape 2!"
585 assert Ah0.shape[0] == nbCtl and Ah0.shape[1] == nbMbr, "Ensemble A0 is not well shaped!"
586 assert abs(max(numpy.mean(Ah0, axis = 1))) < nbMbr*mpr, "Ensemble A0 seems not to be centered!"
588 def _convergence_condition(j, dx, JCurr, JPrev):
590 logging.debug("Convergence on maximum number of iterations per cycle, that reach the limit of %i."%iMaximum)
596 _deltaOnJ = abs(JCurr - JPrev) / JPrev
597 if _deltaOnJ <= jTolerance:
598 logging.debug("Convergence on cost decrement tolerance, that is below the threshold of %.1e."%jTolerance)
601 _deltaOnX = numpy.linalg.norm(dx)
602 if _deltaOnX <= sTolerance:
603 logging.debug("Convergence on norm of state correction, that is below the threshold of %.1e."%sTolerance)
604 return True # En correction de l'état
608 St = dict([(k,[]) for k in [
609 "CurrentState", "CurrentEnsemble",
610 "CostFunctionJb", "CostFunctionJo", "CostFunctionJ",
613 j, convergence, JPrev = 1, False, numpy.nan
615 while not convergence:
616 logging.debug("Internal IEnKS step number %i"%j)
617 St["CurrentState"].append( x1.squeeze() )
618 if variant == "IEnKF": # Transform
621 E1 = x1 + epsilonE * Ah0
622 St["CurrentEnsemble"].append( E1 )
623 E2 = numpy.array([Mnnpu(_x) for _x in E1.T]).reshape((nbCtl, nbMbr)) # Evolution 1->2
624 HEL = numpy.array([Hn(_x) for _x in E2.T]).T # Observation à 2
625 yLm = numpy.mean( HEL, axis = 1).reshape((nbObs,1))
627 if variant == "IEnKF":
631 RIdemidy = RIdemi.dot(yo - yLm)
632 xs = RIdemidy / math.sqrt(nbMbr-1)
633 ES = RIdemi.dot(HA2) / math.sqrt(nbMbr-1)
634 G = numpy.linalg.inv(IN + ES.T.dot(ES))
635 xb = G.dot(ES.T.dot(xs))
636 dx = Ah0.dot(xb) + Ah0.dot(G.dot(Ap0.dot(Ah0.T.dot(x00 - x1))))
638 Jb = float(dx.T.dot(dx))
639 Jo = float(RIdemidy.T.dot(RIdemidy))
641 logging.debug("Values for cost functions are: J = %.5e Jo = %.5e Jb = %.5e"%(J,Jo,Jb))
642 St["CostFunctionJb"].append( Jb )
643 St["CostFunctionJo"].append( Jo )
644 St["CostFunctionJ"].append( J )
648 convergence = _convergence_condition(j, dx, J, JPrev)
651 if variant == "IEnKF":
652 T = numpy.real_if_close(scipy.linalg.sqrtm(G))
653 Tinv = numpy.linalg.inv(T)
655 # Stocke le dernier pas
656 x2 = numpy.mean( E2, axis = 1)
657 if variant == "IEKF":
659 A2 = A2.dot(numpy.linalg.cholesky(G)) / epsilonE
661 St["CurrentState"].append( x2.squeeze() )
662 St["CurrentEnsemble"].append( E2 )
664 IndexMin = numpy.argmin( St["CostFunctionJ"][nbPS:] ) + nbPS
665 xa = St["CurrentState"][IndexMin]
666 Ea = St["CurrentEnsemble"][IndexMin]
671 xb = None, # Background (None si E0)
672 E0 = None, # Background ensemble (None si xb)
673 yObs = None, # Observation (série)
675 RIdemi = None, # R^(-1/2)
676 Mnnpu = None, # Evolution operator
677 Hn = None, # Observation operator
678 variant = "IEnKF", # IEnKF or IEKF
679 nMembers = 5, # Number of members
680 sMaximum = 0, # Number of spinup steps
681 cMaximum = 15000, # Number of steps or cycles
682 iMaximum = 15000, # Number of iterations per cycle
683 sTolerance = mfp, # State correction tolerance
684 jTolerance = mfp, # Cost decrement tolerance
687 nbPS = 0, # Number of previous steps
692 if setSeed is not None: numpy.random.seed(setSeed)
693 if E0 is None: E0 = _BackgroundEnsembleGeneration( xb, B, nMembers)
700 xa, Ea, Sa = [xb,], [E0,], [{}]
701 for step in range(cMaximum):
702 if hasattr(yObs,"store"): Ynpu = numpy.ravel( yObs[step+1] )
703 elif type(yObs) in [list, tuple]: Ynpu = numpy.ravel( yObs[step+1] )
704 else: Ynpu = numpy.ravel( yObs )
706 (xa_c, Ea_c, Sa_c) = _IEnKF_cycle_Lag_1_SDA_GN(
723 # Inflation for next cycle
724 E0 = xa_c + inflation * (Ea_c - xa_c)
728 def _IEnKS_cycle_Lag_L_SDA_GN(
734 method = "Transform",
743 if logging.getLogger().level < logging.WARNING:
744 assert len(E0.shape) == 2, "Ensemble E0 is not well formed: not of shape 2!"
745 assert len(RIdemi.shape) == 2, "R^{-1/2} is not well formed: not of shape 2!"
746 assert method in ("Transform", "Bundle"), "Method has to be Transform or Bundle"
748 nbCtl, nbMbr = E0.shape
751 if logging.getLogger().level < logging.WARNING:
752 assert RIdemi.shape[0] == RIdemi.shape[1] == nbObs, "R^{-1} not of good size: not of size nbObs!"
754 yo = yObs.reshape((nbObs,1))
755 IN = numpy.identity(nbMbr)
756 if method == "Transform":
757 T = numpy.identity(nbMbr)
758 Tinv = numpy.identity(nbMbr)
759 x00 = numpy.mean(E0, axis = 1)
761 Am0 = (1/math.sqrt(nbMbr - 1)) * Ah0
762 w = numpy.zeros((nbMbr,1))
763 if logging.getLogger().level < logging.WARNING:
764 assert len(Ah0.shape) == 2, "Ensemble A0 is not well formed, of shape 2!"
765 assert Ah0.shape[0] == nbCtl and Ah0.shape[1] == nbMbr, "Ensemble A0 is not well shaped!"
766 assert abs(max(numpy.mean(Ah0, axis = 1))) < nbMbr*mpr, "Ensemble A0 seems not to be centered!"
768 def _convergence_condition(j, dw, JCurr, JPrev):
770 logging.debug("Convergence on maximum number of iterations per cycle, that reach the limit of %i."%iMaximum)
776 _deltaOnJ = abs(JCurr - JPrev) / JPrev
777 if _deltaOnJ <= jTolerance:
778 logging.debug("Convergence on cost decrement tolerance, that is below the threshold of %.1e."%jTolerance)
781 _deltaOnW = numpy.sqrt(numpy.mean(dw.squeeze()**2))
782 if _deltaOnW <= sTolerance:
783 logging.debug("Convergence on norm of weights correction, that is below the threshold of %.1e."%sTolerance)
784 return True # En correction des poids
788 St = dict([(k,[]) for k in [
789 "CurrentState", "CurrentEnsemble", "CurrentWeights",
790 "CostFunctionJb", "CostFunctionJo", "CostFunctionJ",
793 j, convergence, JPrev = 1, False, numpy.nan
794 while not convergence:
795 logging.debug("Internal IEnKS step number %i"%j)
796 x0 = x00 + Am0.dot( w )
797 St["CurrentState"].append( x0.squeeze() )
798 if method == "Transform":
801 E0 = x0 + epsilon * Am0
802 St["CurrentEnsemble"].append( E0 )
804 yHmean = numpy.mean(E0, axis = 1)
805 for k in range(1, Lag+1):
806 Ek = numpy.array([Mnnpu(_x) for _x in Ek.T]).reshape((nbCtl, nbMbr)) # Evolution 0->L
807 if method == "Transform":
808 yHmean = Mnnpu(yHmean)
809 HEL = numpy.array([Hn(_x) for _x in Ek.T]).T # Observation à L
811 if method == "Transform":
812 yLm = Hn( yHmean ).reshape((nbObs,1))
813 YL = RIdemi.dot( (HEL - numpy.mean( HEL, axis = 1).reshape((nbObs,1))).dot(Tinv) ) / math.sqrt(nbMbr-1)
815 yLm = numpy.mean( HEL, axis = 1).reshape((nbObs,1))
816 YL = RIdemi.dot(HEL - yLm) / epsilon
817 dy = RIdemi.dot(yo - yLm)
819 Jb = float(w.T.dot(w))
820 Jo = float(dy.T.dot(dy))
822 logging.debug("Values for cost functions are: J = %.5e Jo = %.5e Jb = %.5e"%(J,Jo,Jb))
823 St["CurrentWeights"].append( w.squeeze() )
824 St["CostFunctionJb"].append( Jb )
825 St["CostFunctionJo"].append( Jo )
826 St["CostFunctionJ"].append( J )
827 if method == "Transform":
828 GradJ = w - YL.T.dot(dy)
829 HTild = IN + YL.T.dot(YL)
831 GradJ = (nbMbr - 1)*w - YL.T.dot(RIdemi.dot(dy))
832 HTild = (nbMbr - 1)*IN + YL.T.dot(RIdemi.dot(YL))
833 HTild = numpy.array(HTild, dtype=float)
834 dw = numpy.linalg.solve( HTild, numpy.array(GradJ, dtype=float) )
837 convergence = _convergence_condition(j, dw, J, JPrev)
840 if method == "Transform":
841 (U, s, _) = numpy.linalg.svd(HTild, full_matrices=False) # Hess = U s V
842 T = U.dot(numpy.diag(numpy.sqrt(1./s)).dot(U.T)) # T = Hess^(-1/2)
843 Tinv = U.dot(numpy.diag(numpy.sqrt(s)).dot(U.T)) # Tinv = T^(-1)
845 # Stocke le dernier pas
846 St["CurrentState"].append( numpy.mean( Ek, axis = 1).squeeze() )
847 St["CurrentEnsemble"].append( Ek )
849 IndexMin = numpy.argmin( St["CostFunctionJ"][nbPS:] ) + nbPS
850 xa = St["CurrentState"][IndexMin]
851 Ea = St["CurrentEnsemble"][IndexMin]
856 xb = None, # Background
857 yObs = None, # Observation (série)
858 E0 = None, # Background ensemble
860 RIdemi = None, # R^(-1/2)
861 Mnnpu = None, # Evolution operator
862 Hn = None, # Observation operator
863 method = "Transform", # Bundle ou Transform
864 nMembers = 5, # Number of members
865 cMaximum = 15000, # Number of steps or cycles
866 iMaximum = 15000, # Number of iterations per cycle
867 sTolerance = mfp, # Weights correction tolerance
868 jTolerance = mfp, # Cost decrement tolerance
869 Lag = 1, # Lenght of smoothing window
872 nbPS = 0, # Number of previous steps
877 if setSeed is not None: numpy.random.seed(setSeed)
878 if E0 is None: E0 = _BackgroundEnsembleGeneration( xb, B, nMembers)
885 xa, Ea, Sa = [], [], []
886 for i in range(Lag): # Lag void results
890 for i in range(Lag,cMaximum):
891 (xa_c, Ea_c, Sa_c) = _IEnKS_cycle_Lag_L_SDA_GN(
909 # Inflation for next cycle
910 E0 = xa_c + inflation * (Ea_c - xa_c)
914 # ==============================================================================
915 if __name__ == "__main__":
916 print('\n AUTODIAGNOSTIC\n')