# -*- coding: utf-8 -*-
#
-# Copyright (C) 2008-2021 EDF R&D
+# Copyright (C) 2008-2022 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
import os, time, copy, types, sys, logging
import math, numpy, scipy, scipy.optimize, scipy.version
-from daCore.BasicObjects import Operator
+from daCore.BasicObjects import Operator, Covariance, PartialAlgorithm
from daCore.PlatformInfo import PlatformInfo
mpr = PlatformInfo().MachinePrecision()
mfp = PlatformInfo().MaximumPrecision()
raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
#
if _bgcovariance is None:
- BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers)
+ _Perturbations = numpy.tile( _bgcenter, _nbmembers)
else:
_Z = numpy.random.multivariate_normal(numpy.zeros(_bgcenter.size), _bgcovariance, size=_nbmembers).T
- BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers) + _Z
+ _Perturbations = numpy.tile( _bgcenter, _nbmembers) + _Z
#
- return BackgroundEnsemble
+ return _Perturbations
# ==============================================================================
def EnsembleOfBackgroundPerturbations( _bgcenter, _bgcovariance, _nbmembers, _withSVD = True):
if _nbmembers < 1:
raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
if _bgcovariance is None:
- BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers)
+ _Perturbations = numpy.tile( _bgcenter, _nbmembers)
else:
if _withSVD:
- U, s, V = numpy.linalg.svd(_bgcovariance, full_matrices=False)
+ _U, _s, _V = numpy.linalg.svd(_bgcovariance, full_matrices=False)
_nbctl = _bgcenter.size
if _nbmembers > _nbctl:
_Z = numpy.concatenate((numpy.dot(
- numpy.diag(numpy.sqrt(s[:_nbctl])), V[:_nbctl]),
+ 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])
+ _Z = numpy.dot(numpy.diag(numpy.sqrt(_s[:_nbmembers-1])), _V[:_nbmembers-1])
_Zca = __CenteredRandomAnomalies(_Z, _nbmembers)
- BackgroundEnsemble = _bgcenter + _Zca
+ _Perturbations = _bgcenter + _Zca
else:
if max(abs(_bgcovariance.flatten())) > 0:
_nbctl = _bgcenter.size
_Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1)
_Zca = __CenteredRandomAnomalies(_Z, _nbmembers)
- BackgroundEnsemble = _bgcenter + _Zca
+ _Perturbations = _bgcenter + _Zca
else:
- BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers)
+ _Perturbations = numpy.tile( _bgcenter, _nbmembers)
#
- return BackgroundEnsemble
+ return _Perturbations
# ==============================================================================
def EnsembleMean( __Ensemble ):
if LBounds is not None: # "EstimateProjection" par défaut
Xr = numpy.max(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,0].reshape((-1,1)))),axis=1)
Xr = numpy.min(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,1].reshape((-1,1)))),axis=1)
- Yr = Hm( Xr )
+ Yr = numpy.asarray(Hm( Xr ))
else:
raise ValueError("Quantile simulations has only to be Linear or NonLinear.")
#
# Conserve une valeur par défaut à None s'il n'y a pas de bornes
if __Bounds is None: return None
# Recentre les valeurs numériques de bornes
- return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).transpose((-1,1))
+ return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1,1))
# ==============================================================================
def ApplyBounds( __Vector, __Bounds, __newClip = True):
if not isinstance(__Bounds, numpy.ndarray): # Is an array
raise ValueError("Incorrect array definition of bounds data")
if 2*__Vector.size != __Bounds.size: # Is a 2 column array of vector lenght
- raise ValueError("Incorrect bounds number to be applied for this vector")
+ raise ValueError("Incorrect bounds number (%i) to be applied for this vector (of size %i)"%(__Bounds.size,__Vector.size))
if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
raise ValueError("Incorrectly shaped bounds data")
#
#
return __Vector
+# ==============================================================================
+def Apply3DVarRecentringOnEnsemble(__EnXn, __EnXf, __Ynpu, __HO, __R, __B, __Betaf):
+ "Recentre l'ensemble Xn autour de l'analyse 3DVAR"
+ #
+ Xf = EnsembleMean( __EnXf )
+ Pf = Covariance( asCovariance=EnsembleErrorCovariance(__EnXf) )
+ Pf = (1 - __Betaf) * __B + __Betaf * Pf
+ #
+ selfB = PartialAlgorithm("3DVAR")
+ selfB._parameters["Minimizer"] = "LBFGSB"
+ selfB._parameters["MaximumNumberOfSteps"] = 15000
+ selfB._parameters["CostDecrementTolerance"] = 1.e-7
+ selfB._parameters["ProjectedGradientTolerance"] = -1
+ selfB._parameters["GradientNormTolerance"] = 1.e-05
+ selfB._parameters["StoreInternalVariables"] = False
+ selfB._parameters["optiprint"] = -1
+ selfB._parameters["optdisp"] = 0
+ selfB._parameters["Bounds"] = None
+ selfB._parameters["InitializationPoint"] = Xf
+ std3dvar(selfB, Xf, __Ynpu, None, __HO, None, None, __R, Pf, None)
+ Xa = selfB.get("Analysis")[-1].reshape((-1,1))
+ del selfB
+ #
+ return Xa + EnsembleOfAnomalies( __EnXn )
+
# ==============================================================================
def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
"""
RI = R.getI()
#
__n = Xb.size
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
Xn = Xb
XEtnnpi = numpy.asarray( Mm( (Xnp[:,point], Un) ) ).reshape((-1,1))
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape
- XEtnnpi = XEtnnpi + Cm * Un
+ XEtnnpi = XEtnnpi + Cm @ Un
if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
XEtnnpi = ApplyBounds( XEtnnpi, selfA._parameters["Bounds"] )
elif selfA._parameters["EstimationOf"] == "Parameters":
_Innovation = Ynpu - Yncm.reshape((-1,1))
if selfA._parameters["EstimationOf"] == "Parameters":
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
- _Innovation = _Innovation - Cm * Un
+ _Innovation = _Innovation - Cm @ Un
#
Kn = Pxyn * Pyyn.I
Xn = Xncm.reshape((-1,1)) + Kn * _Innovation
or selfA._toStore("CostFunctionJo") \
or selfA._toStore("CurrentOptimum") \
or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
+ Jb = float( 0.5 * (Xa - Xb).T * (BI * (Xa - Xb)) )
Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
J = Jb + Jo
selfA.StoredVariables["CostFunctionJb"].store( Jb )
RI = R.getI()
#
__n = Xb.size
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
Xn = Xb
Xn_predicted = numpy.ravel( M( (Xn, Un) ) ).reshape((__n,1))
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
+ Xn_predicted = Xn_predicted + Cm @ Un
Pn_predicted = Q + Mt * (Pn * Ma)
elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
# --- > Par principe, M = Id, Q = 0
HX_predicted = numpy.ravel( H( (Xn_predicted, Un) ) ).reshape((__p,1))
_Innovation = Ynpu - HX_predicted
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
- _Innovation = _Innovation - Cm * Un
+ _Innovation = _Innovation - Cm @ Un
#
Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha))
Xn = Xn_predicted + Kn * _Innovation
returnSerieAsArrayMatrix = True )
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- EZ = EZ + Cm * Un
+ EZ = EZ + Cm @ Un
elif selfA._parameters["EstimationOf"] == "Parameters":
# --- > Par principe, M = Id, Q = 0
EZ = H( [(EL[:,i], Un) for i in range(__m)],
return 0
# ==============================================================================
-def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"):
+def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q,
+ VariantM="KalmanFilterFormula",
+ Hybrid=None,
+ ):
"""
Ensemble-Transform EnKF
"""
#
__n = Xb.size
__m = selfA._parameters["NumberOfMembers"]
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
+ previousJMinimum = numpy.finfo(float).max
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
elif selfA._parameters["nextStep"]:
Xn = selfA._getInternalState("Xn")
#
- previousJMinimum = numpy.finfo(float).max
- #
for step in range(duration-1):
numpy.random.set_state(selfA._getInternalState("seed"))
if hasattr(Y,"store"):
returnSerieAsArrayMatrix = True )
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
+ Xn_predicted = Xn_predicted + Cm @ Un
elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
# --- > Par principe, M = Id, Q = 0
Xn_predicted = EMX = Xn
returnSerieAsArrayMatrix = True )
#
# Mean of forecast and observation of forecast
- Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
- Hfm = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
+ Xfm = EnsembleMean( Xn_predicted )
+ Hfm = EnsembleMean( HX_predicted )
#
# Anomalies
EaX = EnsembleOfAnomalies( Xn_predicted, Xfm )
HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
def CostFunction(w):
_A = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
- _Jo = 0.5 * _A.T * RI * _A
+ _Jo = 0.5 * _A.T * (RI * _A)
_Jb = 0.5 * (__m+1) * math.log(1 + 1/__m + w.T @ w)
_J = _Jo + _Jb
return float(_J)
selfA._parameters["InflationFactor"],
)
#
- Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
+ if Hybrid == "E3DVAR":
+ betaf = selfA._parameters["HybridCovarianceEquilibrium"]
+ Xn = Apply3DVarRecentringOnEnsemble(Xn, EMX, Ynpu, HO, R, B, betaf)
+ #
+ Xa = EnsembleMean( Xn )
#--------------------------
selfA._setInternalState("Xn", Xn)
selfA._setInternalState("seed", numpy.random.get_state())
or selfA._toStore("InnovationAtCurrentAnalysis") \
or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
+ _HXa = numpy.ravel( H((Xa, Un)) ).reshape((-1,1))
_Innovation = Ynpu - _HXa
#
selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
or selfA._toStore("CostFunctionJo") \
or selfA._toStore("CurrentOptimum") \
or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
+ Jb = float( 0.5 * (Xa - Xb).T * (BI * (Xa - Xb)) )
+ Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
J = Jb + Jo
selfA.StoredVariables["CostFunctionJb"].store( Jb )
selfA.StoredVariables["CostFunctionJo"].store( Jo )
RI = R.getI()
#
__n = Xb.size
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
Xn = Xb
Xn_predicted = numpy.ravel( M( (Xn, Un) ) ).reshape((__n,1))
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
+ Xn_predicted = Xn_predicted + Cm @ Un
Pn_predicted = Q + Mt * (Pn * Ma)
elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
# --- > Par principe, M = Id, Q = 0
HX_predicted = numpy.ravel( H( (Xn_predicted, Un) ) ).reshape((__p,1))
_Innovation = Ynpu - HX_predicted
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
- _Innovation = _Innovation - Cm * Un
+ _Innovation = _Innovation - Cm @ Un
#
Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha))
Xn = Xn_predicted + Kn * _Innovation
#
__n = Xb.size
__m = selfA._parameters["NumberOfMembers"]
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
+ previousJMinimum = numpy.finfo(float).max
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
elif selfA._parameters["nextStep"]:
Xn = selfA._getInternalState("Xn")
#
- previousJMinimum = numpy.finfo(float).max
- #
for step in range(duration-1):
numpy.random.set_state(selfA._getInternalState("seed"))
if hasattr(Y,"store"):
#
if U is not None:
if hasattr(U,"store") and len(U)>1:
- Un = numpy.asmatrix(numpy.ravel( U[step] )).T
+ Un = numpy.ravel( U[step] ).reshape((-1,1))
elif hasattr(U,"store") and len(U)==1:
- Un = numpy.asmatrix(numpy.ravel( U[0] )).T
+ Un = numpy.ravel( U[0] ).reshape((-1,1))
else:
- Un = numpy.asmatrix(numpy.ravel( U )).T
+ Un = numpy.ravel( U ).reshape((-1,1))
else:
Un = None
#
selfA._parameters["InflationFactor"],
)
#
- Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
+ Xa = EnsembleMean( Xn )
#--------------------------
selfA._setInternalState("Xn", Xn)
selfA._setInternalState("seed", numpy.random.get_state())
or selfA._toStore("InnovationAtCurrentAnalysis") \
or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
+ _HXa = numpy.ravel( H((Xa, Un)) ).reshape((-1,1))
_Innovation = Ynpu - _HXa
#
selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
or selfA._toStore("CostFunctionJo") \
or selfA._toStore("CurrentOptimum") \
or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
+ Jb = float( 0.5 * (Xa - Xb).T * (BI * (Xa - Xb)) )
+ Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
J = Jb + Jo
selfA.StoredVariables["CostFunctionJb"].store( Jb )
selfA.StoredVariables["CostFunctionJo"].store( Jo )
BI = B.getI()
RI = R.getI()
#
- Xini = selfA._parameters["InitializationPoint"]
- #
- HXb = numpy.asmatrix(numpy.ravel( Hm( Xb ) )).T
+ HXb = numpy.asarray(Hm( Xb )).reshape((-1,1))
Innovation = Y - HXb
#
# Outer Loop
iOuter = 0
J = 1./mpr
DeltaJ = 1./mpr
- Xr = Xini.reshape((-1,1))
+ Xr = numpy.asarray(selfA._parameters["InitializationPoint"]).reshape((-1,1))
while abs(DeltaJ) >= selfA._parameters["CostDecrementTolerance"] and iOuter <= selfA._parameters["MaximumNumberOfSteps"]:
#
# Inner Loop
# Définition de la fonction-coût
# ------------------------------
def CostFunction(dx):
- _dX = numpy.asmatrix(numpy.ravel( dx )).T
+ _dX = numpy.asarray(dx).reshape((-1,1))
if selfA._parameters["StoreInternalVariables"] or \
selfA._toStore("CurrentState") or \
selfA._toStore("CurrentOptimum"):
selfA.StoredVariables["CurrentState"].store( Xb + _dX )
- _HdX = Ht * _dX
- _HdX = numpy.asmatrix(numpy.ravel( _HdX )).T
+ _HdX = (Ht @ _dX).reshape((-1,1))
_dInnovation = Innovation - _HdX
if selfA._toStore("SimulatedObservationAtCurrentState") or \
selfA._toStore("SimulatedObservationAtCurrentOptimum"):
if selfA._toStore("InnovationAtCurrentState"):
selfA.StoredVariables["InnovationAtCurrentState"].store( _dInnovation )
#
- Jb = float( 0.5 * _dX.T * BI * _dX )
- Jo = float( 0.5 * _dInnovation.T * RI * _dInnovation )
+ Jb = float( 0.5 * _dX.T * (BI * _dX) )
+ Jo = float( 0.5 * _dInnovation.T * (RI * _dInnovation) )
J = Jb + Jo
#
selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
return J
#
def GradientOfCostFunction(dx):
- _dX = numpy.asmatrix(numpy.ravel( dx )).T
- _HdX = Ht * _dX
- _HdX = numpy.asmatrix(numpy.ravel( _HdX )).T
+ _dX = numpy.ravel( dx )
+ _HdX = (Ht @ _dX).reshape((-1,1))
_dInnovation = Innovation - _HdX
- GradJb = BI * _dX
+ GradJb = BI @ _dX
GradJo = - Ht.T @ (RI * _dInnovation)
GradJ = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
return GradJ
import scipy.optimize as optimiseur
Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
func = CostFunction,
- x0 = numpy.zeros(Xini.size),
+ x0 = numpy.zeros(Xb.size),
fprime = GradientOfCostFunction,
args = (),
bounds = RecentredBounds(selfA._parameters["Bounds"], Xb),
elif selfA._parameters["Minimizer"] == "TNC":
Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
func = CostFunction,
- x0 = numpy.zeros(Xini.size),
+ x0 = numpy.zeros(Xb.size),
fprime = GradientOfCostFunction,
args = (),
bounds = RecentredBounds(selfA._parameters["Bounds"], Xb),
elif selfA._parameters["Minimizer"] == "CG":
Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
f = CostFunction,
- x0 = numpy.zeros(Xini.size),
+ x0 = numpy.zeros(Xb.size),
fprime = GradientOfCostFunction,
args = (),
maxiter = selfA._parameters["MaximumNumberOfSteps"],
elif selfA._parameters["Minimizer"] == "NCG":
Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
f = CostFunction,
- x0 = numpy.zeros(Xini.size),
+ x0 = numpy.zeros(Xb.size),
fprime = GradientOfCostFunction,
args = (),
maxiter = selfA._parameters["MaximumNumberOfSteps"],
elif selfA._parameters["Minimizer"] == "BFGS":
Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
f = CostFunction,
- x0 = numpy.zeros(Xini.size),
+ x0 = numpy.zeros(Xb.size),
fprime = GradientOfCostFunction,
args = (),
maxiter = selfA._parameters["MaximumNumberOfSteps"],
return 0
# ==============================================================================
-def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13",
- BnotT=False, _epsilon=1.e-3, _e=1.e-7, _jmax=15000):
+def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q,
+ VariantM="MLEF13", BnotT=False, _epsilon=1.e-3, _e=1.e-7, _jmax=15000,
+ Hybrid=None,
+ ):
"""
Maximum Likelihood Ensemble Filter
"""
#
__n = Xb.size
__m = selfA._parameters["NumberOfMembers"]
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
+ previousJMinimum = numpy.finfo(float).max
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
elif selfA._parameters["nextStep"]:
Xn = selfA._getInternalState("Xn")
#
- previousJMinimum = numpy.finfo(float).max
- #
for step in range(duration-1):
numpy.random.set_state(selfA._getInternalState("seed"))
if hasattr(Y,"store"):
#
if U is not None:
if hasattr(U,"store") and len(U)>1:
- Un = numpy.asmatrix(numpy.ravel( U[step] )).T
+ Un = numpy.ravel( U[step] ).reshape((-1,1))
elif hasattr(U,"store") and len(U)==1:
- Un = numpy.asmatrix(numpy.ravel( U[0] )).T
+ Un = numpy.ravel( U[0] ).reshape((-1,1))
else:
- Un = numpy.asmatrix(numpy.ravel( U )).T
+ Un = numpy.ravel( U ).reshape((-1,1))
else:
Un = None
#
Xn_predicted = EnsemblePerturbationWithGivenCovariance( EMX, Q )
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
+ Xn_predicted = Xn_predicted + Cm @ Un
elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
# --- > Par principe, M = Id, Q = 0
Xn_predicted = EMX = Xn
selfA._parameters["InflationFactor"],
)
#
- Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
+ if Hybrid == "E3DVAR":
+ betaf = selfA._parameters["HybridCovarianceEquilibrium"]
+ Xn = Apply3DVarRecentringOnEnsemble(Xn, EMX, Ynpu, HO, R, B, betaf)
+ #
+ Xa = EnsembleMean( Xn )
#--------------------------
selfA._setInternalState("Xn", Xn)
selfA._setInternalState("seed", numpy.random.get_state())
or selfA._toStore("InnovationAtCurrentAnalysis") \
or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
+ _HXa = numpy.ravel( H((Xa, Un)) ).reshape((-1,1))
_Innovation = Ynpu - _HXa
#
selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
or selfA._toStore("CostFunctionJo") \
or selfA._toStore("CurrentOptimum") \
or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
+ Jb = float( 0.5 * (Xa - Xb).T * (BI * (Xa - Xb)) )
+ Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
J = Jb + Jo
selfA.StoredVariables["CostFunctionJb"].store( Jb )
selfA.StoredVariables["CostFunctionJo"].store( Jo )
iteration += 1
#
Derivees = numpy.array(fprime(variables))
- Derivees = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
+ Derivees = Derivees.reshape(n,p) # ADAO & check shape
DeriveesT = Derivees.transpose()
M = numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) )
SM = numpy.transpose(numpy.dot( DeriveesT , veps ))
# Initialisation
# --------------
if selfA._parameters["EstimationOf"] == "State":
- M = EM["Direct"].appliedTo
+ M = EM["Direct"].appliedControledFormTo
+ if CM is not None and "Tangent" in CM and U is not None:
+ Cm = CM["Tangent"].asMatrix(Xb)
+ else:
+ Cm = None
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
Xn = numpy.ravel(Xb).reshape((-1,1))
else:
Ynpu = numpy.ravel( Y ).reshape((-1,1))
#
+ if U is not None:
+ if hasattr(U,"store") and len(U)>1:
+ Un = numpy.ravel( U[step] ).reshape((-1,1))
+ elif hasattr(U,"store") and len(U)==1:
+ Un = numpy.ravel( U[0] ).reshape((-1,1))
+ else:
+ Un = numpy.ravel( U ).reshape((-1,1))
+ else:
+ Un = None
+ #
if selfA._parameters["EstimationOf"] == "State": # Forecast
- Xn_predicted = M( Xn )
+ Xn_predicted = M( (Xn, Un) )
if selfA._toStore("ForecastState"):
selfA.StoredVariables["ForecastState"].store( Xn_predicted )
+ if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
+ Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
+ Xn_predicted = Xn_predicted + Cm @ Un
elif selfA._parameters["EstimationOf"] == "Parameters": # No forecast
# --- > Par principe, M = Id, Q = 0
Xn_predicted = Xn
Xn_predicted = numpy.ravel(Xn_predicted).reshape((-1,1))
#
- oneCycle(selfA, Xn_predicted, Ynpu, U, HO, None, None, R, B, None)
+ oneCycle(selfA, Xn_predicted, Ynpu, None, HO, None, None, R, B, None)
#
Xn = selfA.StoredVariables["Analysis"][-1]
#--------------------------
Hm = HO["Direct"].appliedTo
#
if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]:
- HXb = Hm( Xb, HO["AppliedInX"]["HXb"] )
+ HXb = numpy.asarray(Hm( Xb, HO["AppliedInX"]["HXb"] ))
else:
- HXb = Hm( Xb )
- HXb = numpy.asmatrix(numpy.ravel( HXb )).T
+ HXb = numpy.asarray(Hm( Xb ))
+ HXb = numpy.ravel( HXb ).reshape((-1,1))
if Y.size != HXb.size:
raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size))
if max(Y.shape) != max(HXb.shape):
HBHTpR = R + Ht * BHT
Innovation = Y - HXb
#
- Xini = numpy.zeros(Xb.shape)
+ Xini = numpy.zeros(Y.size)
#
# Définition de la fonction-coût
# ------------------------------
def CostFunction(w):
- _W = w.reshape((-1,1))
+ _W = numpy.asarray(w).reshape((-1,1))
if selfA._parameters["StoreInternalVariables"] or \
selfA._toStore("CurrentState") or \
selfA._toStore("CurrentOptimum"):
return J
#
def GradientOfCostFunction(w):
- _W = w.reshape((-1,1))
+ _W = numpy.asarray(w).reshape((-1,1))
GradJb = HBHTpR @ _W
GradJo = - Innovation
GradJ = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
x0 = Xini,
fprime = GradientOfCostFunction,
args = (),
- bounds = RecentredBounds(selfA._parameters["Bounds"], Xb),
maxfun = selfA._parameters["MaximumNumberOfSteps"]-1,
factr = selfA._parameters["CostDecrementTolerance"]*1.e14,
pgtol = selfA._parameters["ProjectedGradientTolerance"],
x0 = Xini,
fprime = GradientOfCostFunction,
args = (),
- bounds = RecentredBounds(selfA._parameters["Bounds"], Xb),
maxfun = selfA._parameters["MaximumNumberOfSteps"],
pgtol = selfA._parameters["ProjectedGradientTolerance"],
ftol = selfA._parameters["CostDecrementTolerance"],
return 0
# ==============================================================================
-def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16"):
+def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q,
+ VariantM="KalmanFilterFormula16",
+ Hybrid=None,
+ ):
"""
Stochastic EnKF
"""
#
__n = Xb.size
__m = selfA._parameters["NumberOfMembers"]
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
+ previousJMinimum = numpy.finfo(float).max
#
if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
else: Rn = R
elif selfA._parameters["nextStep"]:
Xn = selfA._getInternalState("Xn")
#
- previousJMinimum = numpy.finfo(float).max
- #
for step in range(duration-1):
numpy.random.set_state(selfA._getInternalState("seed"))
if hasattr(Y,"store"):
#
if U is not None:
if hasattr(U,"store") and len(U)>1:
- Un = numpy.asmatrix(numpy.ravel( U[step] )).T
+ Un = numpy.ravel( U[step] ).reshape((-1,1))
elif hasattr(U,"store") and len(U)==1:
- Un = numpy.asmatrix(numpy.ravel( U[0] )).T
+ Un = numpy.ravel( U[0] ).reshape((-1,1))
else:
- Un = numpy.asmatrix(numpy.ravel( U )).T
+ Un = numpy.ravel( U ).reshape((-1,1))
else:
Un = None
#
returnSerieAsArrayMatrix = True )
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
+ Xn_predicted = Xn_predicted + Cm @ Un
elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
# --- > Par principe, M = Id, Q = 0
Xn_predicted = EMX = Xn
returnSerieAsArrayMatrix = True )
#
# Mean of forecast and observation of forecast
- Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
- Hfm = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
+ Xfm = EnsembleMean( Xn_predicted )
+ Hfm = EnsembleMean( HX_predicted )
#
#--------------------------
if VariantM == "KalmanFilterFormula05":
selfA._parameters["InflationFactor"],
)
#
- Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
+ if Hybrid == "E3DVAR":
+ betaf = selfA._parameters["HybridCovarianceEquilibrium"]
+ Xn = Apply3DVarRecentringOnEnsemble(Xn, EMX, Ynpu, HO, R, B, betaf)
+ #
+ Xa = EnsembleMean( Xn )
#--------------------------
selfA._setInternalState("Xn", Xn)
selfA._setInternalState("seed", numpy.random.get_state())
or selfA._toStore("InnovationAtCurrentAnalysis") \
or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
+ _HXa = numpy.ravel( H((Xa, Un)) ).reshape((-1,1))
_Innovation = Ynpu - _HXa
#
selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
or selfA._toStore("CostFunctionJo") \
or selfA._toStore("CurrentOptimum") \
or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
+ Jb = float( 0.5 * (Xa - Xb).T * (BI * (Xa - Xb)) )
+ Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
J = Jb + Jo
selfA.StoredVariables["CostFunctionJb"].store( Jb )
selfA.StoredVariables["CostFunctionJo"].store( Jo )
Ha = HO["Adjoint"].appliedInXTo
#
if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]:
- HXb = Hm( Xb, HO["AppliedInX"]["HXb"] )
+ HXb = numpy.asarray(Hm( Xb, HO["AppliedInX"]["HXb"] ))
else:
- HXb = Hm( Xb )
+ HXb = numpy.asarray(Hm( Xb ))
HXb = HXb.reshape((-1,1))
if Y.size != HXb.size:
raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size))
# Définition de la fonction-coût
# ------------------------------
def CostFunction(x):
- _X = numpy.ravel( x ).reshape((-1,1))
+ _X = numpy.asarray(x).reshape((-1,1))
if selfA._parameters["StoreInternalVariables"] or \
selfA._toStore("CurrentState") or \
selfA._toStore("CurrentOptimum"):
selfA.StoredVariables["CurrentState"].store( _X )
- _HX = Hm( _X ).reshape((-1,1))
+ _HX = numpy.asarray(Hm( _X )).reshape((-1,1))
_Innovation = Y - _HX
if selfA._toStore("SimulatedObservationAtCurrentState") or \
selfA._toStore("SimulatedObservationAtCurrentOptimum"):
return J
#
def GradientOfCostFunction(x):
- _X = x.reshape((-1,1))
- _HX = Hm( _X ).reshape((-1,1))
+ _X = numpy.asarray(x).reshape((-1,1))
+ _HX = numpy.asarray(Hm( _X )).reshape((-1,1))
GradJb = BI * (_X - Xb)
GradJo = - Ha( (_X, RI * (Y - _HX)) )
GradJ = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
def Un(_step):
if U is not None:
if hasattr(U,"store") and 1<=_step<len(U) :
- _Un = numpy.asmatrix(numpy.ravel( U[_step] )).T
+ _Un = numpy.ravel( U[_step] ).reshape((-1,1))
elif hasattr(U,"store") and len(U)==1:
- _Un = numpy.asmatrix(numpy.ravel( U[0] )).T
+ _Un = numpy.ravel( U[0] ).reshape((-1,1))
else:
- _Un = numpy.asmatrix(numpy.ravel( U )).T
+ _Un = numpy.ravel( U ).reshape((-1,1))
else:
_Un = None
return _Un
def CmUn(_xn,_un):
if Cm is not None and _un is not None: # Attention : si Cm est aussi dans M, doublon !
_Cm = Cm.reshape(_xn.size,_un.size) # ADAO & check shape
- _CmUn = _Cm * _un
+ _CmUn = (_Cm @ _un).reshape((-1,1))
else:
_CmUn = 0.
return _CmUn
selfA.DirectCalculation = [None,] # Le pas 0 n'est pas observé
selfA.DirectInnovation = [None,] # Le pas 0 n'est pas observé
def CostFunction(x):
- _X = numpy.asmatrix(numpy.ravel( x )).T
+ _X = numpy.asarray(x).reshape((-1,1))
if selfA._parameters["StoreInternalVariables"] or \
selfA._toStore("CurrentState") or \
selfA._toStore("CurrentOptimum"):
selfA.StoredVariables["CurrentState"].store( _X )
- Jb = float( 0.5 * (_X - Xb).T * BI * (_X - Xb) )
+ Jb = float( 0.5 * (_X - Xb).T * (BI * (_X - Xb)) )
selfA.DirectCalculation = [None,]
selfA.DirectInnovation = [None,]
Jo = 0.
_Xn = _X
for step in range(0,duration-1):
if hasattr(Y,"store"):
- _Ynpu = numpy.asmatrix(numpy.ravel( Y[step+1] )).T
+ _Ynpu = numpy.ravel( Y[step+1] ).reshape((-1,1))
else:
- _Ynpu = numpy.asmatrix(numpy.ravel( Y )).T
+ _Ynpu = numpy.ravel( Y ).reshape((-1,1))
_Un = Un(step)
#
# Etape d'évolution
if selfA._parameters["EstimationOf"] == "State":
- _Xn = Mm( (_Xn, _Un) ) + CmUn(_Xn, _Un)
+ _Xn = Mm( (_Xn, _Un) ).reshape((-1,1)) + CmUn(_Xn, _Un)
elif selfA._parameters["EstimationOf"] == "Parameters":
pass
#
#
# Etape de différence aux observations
if selfA._parameters["EstimationOf"] == "State":
- _YmHMX = _Ynpu - numpy.asmatrix(numpy.ravel( Hm( (_Xn, None) ) )).T
+ _YmHMX = _Ynpu - numpy.ravel( Hm( (_Xn, None) ) ).reshape((-1,1))
elif selfA._parameters["EstimationOf"] == "Parameters":
- _YmHMX = _Ynpu - numpy.asmatrix(numpy.ravel( Hm( (_Xn, _Un) ) )).T - CmUn(_Xn, _Un)
+ _YmHMX = _Ynpu - numpy.ravel( Hm( (_Xn, _Un) ) ).reshape((-1,1)) - CmUn(_Xn, _Un)
#
# Stockage de l'état
selfA.DirectCalculation.append( _Xn )
selfA.DirectInnovation.append( _YmHMX )
#
# Ajout dans la fonctionnelle d'observation
- Jo = Jo + 0.5 * float( _YmHMX.T * RI * _YmHMX )
+ Jo = Jo + 0.5 * float( _YmHMX.T * (RI * _YmHMX) )
J = Jb + Jo
#
selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
return J
#
def GradientOfCostFunction(x):
- _X = numpy.asmatrix(numpy.ravel( x )).T
+ _X = numpy.asarray(x).reshape((-1,1))
GradJb = BI * (_X - Xb)
GradJo = 0.
for step in range(duration-1,0,-1):
RI = R.getI()
#
__n = Xb.size
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
Xn = Xb
#
if U is not None:
if hasattr(U,"store") and len(U)>1:
- Un = numpy.asmatrix(numpy.ravel( U[step] )).T
+ Un = numpy.ravel( U[step] ).reshape((-1,1))
elif hasattr(U,"store") and len(U)==1:
- Un = numpy.asmatrix(numpy.ravel( U[0] )).T
+ Un = numpy.ravel( U[0] ).reshape((-1,1))
else:
- Un = numpy.asmatrix(numpy.ravel( U )).T
+ Un = numpy.ravel( U ).reshape((-1,1))
else:
Un = None
#
if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
- Xn_predicted = Mt * Xn
+ Xn_predicted = Mt @ Xn
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
+ Xn_predicted = Xn_predicted + Cm @ Un
Pn_predicted = Q + Mt * (Pn * Ma)
elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
# --- > Par principe, M = Id, Q = 0
Pn_predicted = Pn
#
if selfA._parameters["EstimationOf"] == "State":
- HX_predicted = Ht * Xn_predicted
+ HX_predicted = Ht @ Xn_predicted
_Innovation = Ynpu - HX_predicted
elif selfA._parameters["EstimationOf"] == "Parameters":
- HX_predicted = Ht * Xn_predicted
+ HX_predicted = Ht @ Xn_predicted
_Innovation = Ynpu - HX_predicted
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
- _Innovation = _Innovation - Cm * Un
+ _Innovation = _Innovation - Cm @ Un
#
Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha))
Xn = Xn_predicted + Kn * _Innovation
or selfA._toStore("CostFunctionJo") \
or selfA._toStore("CurrentOptimum") \
or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
+ Jb = float( 0.5 * (Xa - Xb).T * (BI * (Xa - Xb)) )
+ Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
J = Jb + Jo
selfA.StoredVariables["CostFunctionJb"].store( Jb )
selfA.StoredVariables["CostFunctionJo"].store( Jo )
RI = R.getI()
#
__n = Xb.size
+ nbPreviousSteps = len(selfA.StoredVariables["Analysis"])
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
Xn = Xb
XEtnnpi = numpy.asarray( Mm( (Xnp[:,point], Un) ) ).reshape((-1,1))
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape
- XEtnnpi = XEtnnpi + Cm * Un
+ XEtnnpi = XEtnnpi + Cm @ Un
elif selfA._parameters["EstimationOf"] == "Parameters":
# --- > Par principe, M = Id, Q = 0
XEtnnpi = Xnp[:,point]
_Innovation = Ynpu - Yncm.reshape((-1,1))
if selfA._parameters["EstimationOf"] == "Parameters":
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
- _Innovation = _Innovation - Cm * Un
+ _Innovation = _Innovation - Cm @ Un
#
Kn = Pxyn * Pyyn.I
Xn = Xncm.reshape((-1,1)) + Kn * _Innovation
or selfA._toStore("CostFunctionJo") \
or selfA._toStore("CurrentOptimum") \
or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
+ Jb = float( 0.5 * (Xa - Xb).T * (BI * (Xa - Xb)) )
Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
J = Jb + Jo
selfA.StoredVariables["CostFunctionJb"].store( Jb )
BT = B.getT()
RI = R.getI()
#
- Xini = numpy.zeros(Xb.shape)
+ Xini = numpy.zeros(Xb.size)
#
# Définition de la fonction-coût
# ------------------------------
def CostFunction(v):
- _V = numpy.asmatrix(numpy.ravel( v )).T
- _X = Xb + B * _V
+ _V = numpy.asarray(v).reshape((-1,1))
+ _X = Xb + (B @ _V).reshape((-1,1))
if selfA._parameters["StoreInternalVariables"] or \
selfA._toStore("CurrentState") or \
selfA._toStore("CurrentOptimum"):
selfA.StoredVariables["CurrentState"].store( _X )
- _HX = Hm( _X )
- _HX = numpy.asmatrix(numpy.ravel( _HX )).T
+ _HX = numpy.asarray(Hm( _X )).reshape((-1,1))
_Innovation = Y - _HX
if selfA._toStore("SimulatedObservationAtCurrentState") or \
selfA._toStore("SimulatedObservationAtCurrentOptimum"):
if selfA._toStore("InnovationAtCurrentState"):
selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
#
- Jb = float( 0.5 * _V.T * BT * _V )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
+ Jb = float( 0.5 * _V.T * (BT * _V) )
+ Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
J = Jb + Jo
#
selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
return J
#
def GradientOfCostFunction(v):
- _V = v.reshape((-1,1))
+ _V = numpy.asarray(v).reshape((-1,1))
_X = Xb + (B @ _V).reshape((-1,1))
- _HX = Hm( _X ).reshape((-1,1))
+ _HX = numpy.asarray(Hm( _X )).reshape((-1,1))
GradJb = BT * _V
GradJo = - Ha( (_X, RI * (Y - _HX)) )
GradJ = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
