#
return variables, Ecart, [n,p,iteration,increment,0]
-# ==============================================================================
-
-def _BackgroundEnsembleGeneration( _bgcenter, _bgcovariance, _nbmembers, _withSVD = True):
- "Génération d'un ensemble d'ébauche de taille _nbmembers-1"
- # ~ numpy.random.seed(1234567)
- if _nbmembers < 1:
- raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
- if _withSVD:
- U, s, V = numpy.linalg.svd(_bgcovariance, full_matrices=False)
- _nbctl = len(_bgcenter)
- if _nbmembers > _nbctl:
- _Z = numpy.concatenate((numpy.dot(
- numpy.diag(numpy.sqrt(s[:_nbctl])), V[:_nbctl]),
- numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1-_nbctl)), axis = 0)
- else:
- _Z = numpy.dot(numpy.diag(numpy.sqrt(s[:_nbmembers-1])), V[:_nbmembers-1])
- _Zca = _CenteredAnomalies(_Z, _nbmembers)
- BackgroundEnsemble = (_bgcenter + _Zca.T).T
- else:
- if max(abs(_bgcovariance.flatten())) > 0:
- _nbctl = len(_bgcenter)
- _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1)
- _Zca = _CenteredAnomalies(_Z, _nbmembers)
- BackgroundEnsemble = (_bgcenter + _Zca.T).T
- else:
- BackgroundEnsemble = numpy.tile([_bgcenter],(_nbmembers,1)).T
- return BackgroundEnsemble
-
-def _CenteredAnomalies(Zr, N):
- """
- Génère une matrice d'anomalies centrées selon les notes manuscrites de MB
- et conforme au code de PS avec eps = -1
- """
- eps = -1
- Q = numpy.eye(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
- Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
- R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
- Q = numpy.dot(Q,R)
- Zr = numpy.dot(Q,Zr)
- return Zr.T
-
-def _IEnKF_cycle_Lag_1_SDA_GN(
- E0 = None,
- yObs = None,
- RIdemi = None,
- Mnnpu = None,
- Hn = None,
- variant = "IEnKF", # IEnKF or IEKF
- iMaximum = 15000,
- sTolerance = mfp,
- jTolerance = mfp,
- epsilonE = 1e-5,
- nbPS = 0, # nbPreviousSteps
- ):
- # 201206
- if logging.getLogger().level < logging.WARNING:
- assert len(E0.shape) == 2, "Ensemble E0 is not well formed: not of shape 2!"
- assert len(RIdemi.shape) == 2, "R^{-1/2} is not well formed: not of shape 2!"
- assert variant in ("IEnKF", "IEKF"), "Variant has to be IEnKF or IEKF"
- #
- nbCtl, nbMbr = E0.shape
- nbObs = yObs.size
- #
- if logging.getLogger().level < logging.WARNING:
- assert RIdemi.shape[0] == RIdemi.shape[1] == nbObs, "R^{-1} not of good size: not of size nbObs!"
- #
- yo = yObs.reshape((nbObs,1))
- IN = numpy.identity(nbMbr)
- if variant == "IEnKF":
- T = numpy.identity(nbMbr)
- Tinv = numpy.identity(nbMbr)
- x00 = numpy.mean(E0, axis = 1)
- Ah0 = E0 - x00
- Ap0 = numpy.linalg.pinv( Ah0.T.dot(Ah0) )
- if logging.getLogger().level < logging.WARNING:
- assert len(Ah0.shape) == 2, "Ensemble A0 is not well formed, of shape 2!"
- assert Ah0.shape[0] == nbCtl and Ah0.shape[1] == nbMbr, "Ensemble A0 is not well shaped!"
- assert abs(max(numpy.mean(Ah0, axis = 1))) < nbMbr*mpr, "Ensemble A0 seems not to be centered!"
- #
- def _convergence_condition(j, dx, JCurr, JPrev):
- if j > iMaximum:
- logging.debug("Convergence on maximum number of iterations per cycle, that reach the limit of %i."%iMaximum)
- return True
- #---------
- if j == 1:
- _deltaOnJ = 1.
- else:
- _deltaOnJ = abs(JCurr - JPrev) / JPrev
- if _deltaOnJ <= jTolerance:
- logging.debug("Convergence on cost decrement tolerance, that is below the threshold of %.1e."%jTolerance)
- return True
- #---------
- _deltaOnX = numpy.linalg.norm(dx)
- if _deltaOnX <= sTolerance:
- logging.debug("Convergence on norm of state correction, that is below the threshold of %.1e."%sTolerance)
- return True # En correction de l'état
- #---------
- return False
- #
- St = dict([(k,[]) for k in [
- "CurrentState", "CurrentEnsemble",
- "CostFunctionJb", "CostFunctionJo", "CostFunctionJ",
- ]])
- #
- j, convergence, JPrev = 1, False, numpy.nan
- x1 = x00
- while not convergence:
- logging.debug("Internal IEnKS step number %i"%j)
- St["CurrentState"].append( x1.squeeze() )
- if variant == "IEnKF": # Transform
- E1 = x1 + Ah0.dot(T)
- else: # IEKF
- E1 = x1 + epsilonE * Ah0
- St["CurrentEnsemble"].append( E1 )
- E2 = numpy.array([Mnnpu(_x) for _x in E1.T]).reshape((nbCtl, nbMbr)) # Evolution 1->2
- HEL = numpy.array([Hn(_x) for _x in E2.T]).T # Observation à 2
- yLm = numpy.mean( HEL, axis = 1).reshape((nbObs,1))
- HA2 = HEL - yLm
- if variant == "IEnKF":
- HA2 = HA2.dot(Tinv)
- else:
- HA2 = HA2 / epsilonE
- RIdemidy = RIdemi.dot(yo - yLm)
- xs = RIdemidy / math.sqrt(nbMbr-1)
- ES = RIdemi.dot(HA2) / math.sqrt(nbMbr-1)
- G = numpy.linalg.inv(IN + ES.T.dot(ES))
- xb = G.dot(ES.T.dot(xs))
- dx = Ah0.dot(xb) + Ah0.dot(G.dot(Ap0.dot(Ah0.T.dot(x00 - x1))))
- #
- Jb = float(dx.T.dot(dx))
- Jo = float(RIdemidy.T.dot(RIdemidy))
- J = Jo + Jb
- logging.debug("Values for cost functions are: J = %.5e Jo = %.5e Jb = %.5e"%(J,Jo,Jb))
- St["CostFunctionJb"].append( Jb )
- St["CostFunctionJo"].append( Jo )
- St["CostFunctionJ"].append( J )
- #
- x1 = x1 + dx
- j = j + 1
- convergence = _convergence_condition(j, dx, J, JPrev)
- JPrev = J
- #
- if variant == "IEnKF":
- T = numpy.real_if_close(scipy.linalg.sqrtm(G))
- Tinv = numpy.linalg.inv(T)
- #
- # Stocke le dernier pas
- x2 = numpy.mean( E2, axis = 1)
- if variant == "IEKF":
- A2 = E2 - x2
- A2 = A2.dot(numpy.linalg.cholesky(G)) / epsilonE
- E2 = x2 + A2
- St["CurrentState"].append( x2.squeeze() )
- St["CurrentEnsemble"].append( E2 )
- #
- IndexMin = numpy.argmin( St["CostFunctionJ"][nbPS:] ) + nbPS
- xa = St["CurrentState"][IndexMin]
- Ea = St["CurrentEnsemble"][IndexMin]
- #
- return (xa, Ea, St)
-
-def ienkf(
- xb = None, # Background (None si E0)
- E0 = None, # Background ensemble (None si xb)
- yObs = None, # Observation (série)
- B = None, # B
- RIdemi = None, # R^(-1/2)
- Mnnpu = None, # Evolution operator
- Hn = None, # Observation operator
- variant = "IEnKF", # IEnKF or IEKF
- nMembers = 5, # Number of members
- sMaximum = 0, # Number of spinup steps
- cMaximum = 15000, # Number of steps or cycles
- iMaximum = 15000, # Number of iterations per cycle
- sTolerance = mfp, # State correction tolerance
- jTolerance = mfp, # Cost decrement tolerance
- epsilon = 1e-5,
- inflation = 1.,
- nbPS = 0, # Number of previous steps
- setSeed = None,
- ):
- #
- # Initial
- if setSeed is not None: numpy.random.seed(setSeed)
- if E0 is None: E0 = _BackgroundEnsembleGeneration( xb, B, nMembers)
- #
- # Spinup
- # ------
- #
- # Cycles
- # ------
- xa, Ea, Sa = [xb,], [E0,], [{}]
- for step in range(cMaximum):
- if hasattr(yObs,"store"): Ynpu = numpy.ravel( yObs[step+1] )
- elif type(yObs) in [list, tuple]: Ynpu = numpy.ravel( yObs[step+1] )
- else: Ynpu = numpy.ravel( yObs )
- #
- (xa_c, Ea_c, Sa_c) = _IEnKF_cycle_Lag_1_SDA_GN(
- E0,
- Ynpu,
- RIdemi,
- Mnnpu,
- Hn,
- variant,
- iMaximum,
- sTolerance,
- jTolerance,
- epsilon,
- nbPS,
- )
- xa.append( xa_c )
- Ea.append( Ea_c )
- Sa.append( Sa_c )
- #
- # Inflation for next cycle
- E0 = xa_c + inflation * (Ea_c - xa_c)
- #
- return (xa, Ea, Sa)
-
-def _IEnKS_cycle_Lag_L_SDA_GN(
- E0 = None,
- yObs = None,
- RIdemi = None,
- Mnnpu = None,
- Hn = None,
- method = "Transform",
- iMaximum = 15000,
- sTolerance = mfp,
- jTolerance = mfp,
- Lag = 1,
- epsilon = -1.,
- nbPS = 0,
- ):
- # 201407 & 201905
- if logging.getLogger().level < logging.WARNING:
- assert len(E0.shape) == 2, "Ensemble E0 is not well formed: not of shape 2!"
- assert len(RIdemi.shape) == 2, "R^{-1/2} is not well formed: not of shape 2!"
- assert method in ("Transform", "Bundle"), "Method has to be Transform or Bundle"
- #
- nbCtl, nbMbr = E0.shape
- nbObs = yObs.size
- #
- if logging.getLogger().level < logging.WARNING:
- assert RIdemi.shape[0] == RIdemi.shape[1] == nbObs, "R^{-1} not of good size: not of size nbObs!"
- #
- yo = yObs.reshape((nbObs,1))
- IN = numpy.identity(nbMbr)
- if method == "Transform":
- T = numpy.identity(nbMbr)
- Tinv = numpy.identity(nbMbr)
- x00 = numpy.mean(E0, axis = 1)
- Ah0 = E0 - x00
- Am0 = (1/math.sqrt(nbMbr - 1)) * Ah0
- w = numpy.zeros((nbMbr,1))
- if logging.getLogger().level < logging.WARNING:
- assert len(Ah0.shape) == 2, "Ensemble A0 is not well formed, of shape 2!"
- assert Ah0.shape[0] == nbCtl and Ah0.shape[1] == nbMbr, "Ensemble A0 is not well shaped!"
- assert abs(max(numpy.mean(Ah0, axis = 1))) < nbMbr*mpr, "Ensemble A0 seems not to be centered!"
- #
- def _convergence_condition(j, dw, JCurr, JPrev):
- if j > iMaximum:
- logging.debug("Convergence on maximum number of iterations per cycle, that reach the limit of %i."%iMaximum)
- return True
- #---------
- if j == 1:
- _deltaOnJ = 1.
- else:
- _deltaOnJ = abs(JCurr - JPrev) / JPrev
- if _deltaOnJ <= jTolerance:
- logging.debug("Convergence on cost decrement tolerance, that is below the threshold of %.1e."%jTolerance)
- return True
- #---------
- _deltaOnW = numpy.sqrt(numpy.mean(dw.squeeze()**2))
- if _deltaOnW <= sTolerance:
- logging.debug("Convergence on norm of weights correction, that is below the threshold of %.1e."%sTolerance)
- return True # En correction des poids
- #---------
- return False
- #
- St = dict([(k,[]) for k in [
- "CurrentState", "CurrentEnsemble", "CurrentWeights",
- "CostFunctionJb", "CostFunctionJo", "CostFunctionJ",
- ]])
- #
- j, convergence, JPrev = 1, False, numpy.nan
- while not convergence:
- logging.debug("Internal IEnKS step number %i"%j)
- x0 = x00 + Am0.dot( w )
- St["CurrentState"].append( x0.squeeze() )
- if method == "Transform":
- E0 = x0 + Ah0.dot(T)
- else:
- E0 = x0 + epsilon * Am0
- St["CurrentEnsemble"].append( E0 )
- Ek = E0
- yHmean = numpy.mean(E0, axis = 1)
- for k in range(1, Lag+1):
- Ek = numpy.array([Mnnpu(_x) for _x in Ek.T]).reshape((nbCtl, nbMbr)) # Evolution 0->L
- if method == "Transform":
- yHmean = Mnnpu(yHmean)
- HEL = numpy.array([Hn(_x) for _x in Ek.T]).T # Observation à L
- #
- if method == "Transform":
- yLm = Hn( yHmean ).reshape((nbObs,1))
- YL = RIdemi.dot( (HEL - numpy.mean( HEL, axis = 1).reshape((nbObs,1))).dot(Tinv) ) / math.sqrt(nbMbr-1)
- else:
- yLm = numpy.mean( HEL, axis = 1).reshape((nbObs,1))
- YL = RIdemi.dot(HEL - yLm) / epsilon
- dy = RIdemi.dot(yo - yLm)
- #
- Jb = float(w.T.dot(w))
- Jo = float(dy.T.dot(dy))
- J = Jo + Jb
- logging.debug("Values for cost functions are: J = %.5e Jo = %.5e Jb = %.5e"%(J,Jo,Jb))
- St["CurrentWeights"].append( w.squeeze() )
- St["CostFunctionJb"].append( Jb )
- St["CostFunctionJo"].append( Jo )
- St["CostFunctionJ"].append( J )
- if method == "Transform":
- GradJ = w - YL.T.dot(dy)
- HTild = IN + YL.T.dot(YL)
- else:
- GradJ = (nbMbr - 1)*w - YL.T.dot(RIdemi.dot(dy))
- HTild = (nbMbr - 1)*IN + YL.T.dot(RIdemi.dot(YL))
- HTild = numpy.array(HTild, dtype=float)
- dw = numpy.linalg.solve( HTild, numpy.array(GradJ, dtype=float) )
- w = w - dw
- j = j + 1
- convergence = _convergence_condition(j, dw, J, JPrev)
- JPrev = J
- #
- if method == "Transform":
- (U, s, _) = numpy.linalg.svd(HTild, full_matrices=False) # Hess = U s V
- T = U.dot(numpy.diag(numpy.sqrt(1./s)).dot(U.T)) # T = Hess^(-1/2)
- Tinv = U.dot(numpy.diag(numpy.sqrt(s)).dot(U.T)) # Tinv = T^(-1)
- #
- # Stocke le dernier pas
- St["CurrentState"].append( numpy.mean( Ek, axis = 1).squeeze() )
- St["CurrentEnsemble"].append( Ek )
- #
- IndexMin = numpy.argmin( St["CostFunctionJ"][nbPS:] ) + nbPS
- xa = St["CurrentState"][IndexMin]
- Ea = St["CurrentEnsemble"][IndexMin]
- #
- return (xa, Ea, St)
-
-def ienks(
- xb = None, # Background
- yObs = None, # Observation (série)
- E0 = None, # Background ensemble
- B = None, # B
- RIdemi = None, # R^(-1/2)
- Mnnpu = None, # Evolution operator
- Hn = None, # Observation operator
- method = "Transform", # Bundle ou Transform
- nMembers = 5, # Number of members
- cMaximum = 15000, # Number of steps or cycles
- iMaximum = 15000, # Number of iterations per cycle
- sTolerance = mfp, # Weights correction tolerance
- jTolerance = mfp, # Cost decrement tolerance
- Lag = 1, # Lenght of smoothing window
- epsilon = -1.,
- inflation = 1.,
- nbPS = 0, # Number of previous steps
- setSeed = None,
- ):
- #
- # Initial
- if setSeed is not None: numpy.random.seed(setSeed)
- if E0 is None: E0 = _BackgroundEnsembleGeneration( xb, B, nMembers)
- #
- # Spinup
- # ------
- #
- # Cycles
- # ------
- xa, Ea, Sa = [], [], []
- for i in range(Lag): # Lag void results
- xa.append([])
- Ea.append([])
- Sa.append([])
- for i in range(Lag,cMaximum):
- (xa_c, Ea_c, Sa_c) = _IEnKS_cycle_Lag_L_SDA_GN(
- E0,
- yObs[i-Lag:i],
- RIdemi,
- Mnnpu,
- Hn,
- method,
- iMaximum,
- sTolerance,
- jTolerance,
- Lag,
- epsilon,
- nbPS,
- )
- xa.append( xa_c )
- Ea.append( Ea_c )
- Sa.append( Sa_c )
- #
- # Inflation for next cycle
- E0 = xa_c + inflation * (Ea_c - xa_c)
- #
- return (xa, Ea, Sa)
-
# ==============================================================================
if __name__ == "__main__":
print('\n AUTODIAGNOSTIC\n')
