.. include:: snippets/EstimationOf.rst
+.. include:: snippets/InflationFactor.rst
+
+.. include:: snippets/InflationType.rst
+
.. include:: snippets/NumberOfMembers.rst
.. include:: snippets/SetSeed.rst
--- /dev/null
+.. index:: single: InflationFactor
+
+InflationFactor
+ *Real value*. This key specifies the inflation factor in the ensemble
+ methods, to be applied on the covariance or the anomalies depending on the
+ choice of the type of inflation. Its value must be positive if the inflation
+ is additive, or greater than 1 if the inflation is multiplicative. The
+ default value is 1, which leads to an absence of multiplicative inflation.
+ The absence of additive inflation is obtained by entering a value of 0.
+
+ Example :
+ ``{"InflationFactor":1.}``
--- /dev/null
+.. index:: single: InflationType
+
+InflationType
+ *Predefined name*. This key is used to set the inflation method in ensemble
+ methods, for those that require such a technique. Inflation can be applied in
+ various ways, according to the following options: multiplicative or additive
+ by the specified inflation factor, applied on the background or on the
+ analysis, applied on covariances or on anomalies. Only one type of inflation
+ is applied at the same time, and the default value is
+ "MultiplicativeOnAnalysisAnomalies". The possible names are in the following
+ list: [
+ "MultiplicativeOnAnalysisAnomalies",
+ "MultiplicativeOnBackgroundAnomalies",
+ ].
+
+ Example :
+ ``{"InflationType":"MultiplicativeOnAnalysisAnomalies"}``
.. include:: snippets/EstimationOf.rst
+.. include:: snippets/InflationFactor.rst
+
+.. include:: snippets/InflationType.rst
+
.. include:: snippets/NumberOfMembers.rst
.. include:: snippets/SetSeed.rst
ConstrainedBy
*Nom prédéfini*. Cette clé permet d'indiquer la méthode de prise en compte
des contraintes de bornes. La seule disponible est "EstimateProjection", qui
- projete l'estimation de l'état courant sur les contraintes de bornes.
+ projette l'estimation de l'état courant sur les contraintes de bornes.
Exemple :
``{"ConstrainedBy":"EstimateProjection"}``
--- /dev/null
+.. index:: single: InflationFactor
+
+InflationFactor
+ *Valeur réelle*. Cette clé indique le facteur d'inflation dans les méthodes
+ d'ensemble, à appliquer sur la covariance ou les anomalies selon le choix du
+ type d'inflation. Sa valeur doit être positive si l'inflation est additive,
+ ou supérieure à 1 si l'inflation est multiplicative. La valeur par défaut est
+ 1, qui conduit à une absence d'inflation multiplicative. L'absence
+ d'inflation additive est obtenue en indiquant une valeur de 0.
+
+ Exemple :
+ ``{"InflationFactor":1.}``
--- /dev/null
+.. index:: single: InflationType
+
+InflationType
+ *Nom prédéfini*. Cette clé permet d'indiquer la méthode d'inflation dans les
+ méthodes d'ensemble, pour celles qui nécessitent une telle technique.
+ L'inflation peut être appliquée de diverses manières, selon les options
+ suivantes : multiplicative ou additive du facteur d'inflation spécifié,
+ appliquée sur l'ébauche ou sur l'analyse, appliquée sur les covariances ou
+ sur les anomalies. Un seul type d'inflation est appliqué à la fois, et la
+ valeur par défaut est "MultiplicativeOnAnalysisAnomalies". Les noms possibles
+ sont dans la liste suivante : [
+ "MultiplicativeOnAnalysisAnomalies",
+ "MultiplicativeOnBackgroundAnomalies",
+ ].
+
+ Exemple :
+ ``{"InflationType":"MultiplicativeOnAnalysisAnomalies"}``
)
self.defineRequiredParameter(
name = "InflationType",
- default = "MultiplicativeOnAnalysisCovariance",
+ default = "MultiplicativeOnAnalysisAnomalies",
typecast = str,
message = "Méthode d'inflation d'ensemble",
listval = [
- "MultiplicativeOnAnalysisCovariance",
- "MultiplicativeOnBackgroundCovariance",
"MultiplicativeOnAnalysisAnomalies",
"MultiplicativeOnBackgroundAnomalies",
+ ],
+ listadv = [
+ "MultiplicativeOnAnalysisCovariance",
+ "MultiplicativeOnBackgroundCovariance",
"AdditiveOnAnalysisCovariance",
"AdditiveOnBackgroundCovariance",
"HybridOnBackgroundCovariance",
+ "Relaxation",
],
)
self.defineRequiredParameter(
minval = 0.,
)
self.defineRequiredParameter(
- name = "LocalizationType",
- default = "SchurLocalization",
- typecast = str,
- message = "Méthode d'inflation d'ensemble",
- listval = [
- "SchurLocalization",
- ],
- listadv = [
- "CovarianceLocalization",
- "DomainLocalization",
- "GaspariCohnLocalization",
- ],
+ name = "SmootherLagL",
+ default = 1,
+ typecast = int,
+ message = "Nombre d'intervalles de temps de lissage dans le passé",
+ minval = 1,
)
self.defineRequiredParameter(
- name = "LocalizationFactor",
- default = 1.,
- typecast = float,
- message = "Facteur de localisation",
- minval = 0.,
- )
- self.defineRequiredParameter( # Pas de type
- name = "LocalizationMatrix",
- message = "Matrice de localisation ou de distances",
+ name = "SmootherShiftS",
+ default = 1,
+ typecast = int,
+ message = "Nombre d'intervalles de temps de décalage de lissage",
+ minval = 1,
)
self.defineRequiredParameter(
name = "SetSeed",
"SimulatedObservationAtCurrentAnalysis",
"SimulatedObservationAtCurrentOptimum",
"SimulatedObservationAtCurrentState",
- ]
+ ],
+ listadv = [
+ "CurrentEnsembleState",
+ "LastEnsembleForecastState",
+ ],
)
self.requireInputArguments(
mandatory= ("Xb", "Y", "HO", "R", "B"),
self.StoredVariables["CostFunctionJbAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJbAtCurrentOptimum")
self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
self.StoredVariables["CostFunctionJoAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJoAtCurrentOptimum")
- self.StoredVariables["CurrentIterationNumber"] = Persistence.OneIndex(name = "CurrentIterationNumber")
+ self.StoredVariables["CurrentEnsembleState"] = Persistence.OneMatrix(name = "CurrentEnsembleState")
+ self.StoredVariables["CurrentIterationNumber"] = Persistence.OneIndex(name = "CurrentIterationNumber")
self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
self.StoredVariables["ForecastState"] = Persistence.OneVector(name = "ForecastState")
self.StoredVariables["JacobianMatrixAtCurrentState"] = Persistence.OneMatrix(name = "JacobianMatrixAtCurrentState")
self.StoredVariables["JacobianMatrixAtOptimum"] = Persistence.OneMatrix(name = "JacobianMatrixAtOptimum")
self.StoredVariables["KalmanGainAtOptimum"] = Persistence.OneMatrix(name = "KalmanGainAtOptimum")
+ self.StoredVariables["LastEnsembleForecastState"] = Persistence.OneMatrix(name = "LastEnsembleForecastState")
self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
def getI(self):
"Inversion"
if self.ismatrix():
- return Covariance(self.__name+"I", asCovariance = self.__C.I )
+ return Covariance(self.__name+"I", asCovariance = numpy.linalg.inv(self.__C) )
elif self.isvector():
return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
elif self.isscalar():
return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
- elif self.isobject():
+ elif self.isobject() and hasattr(self.__C,"getI"):
return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
else:
return None # Indispensable
return Covariance(self.__name+"T", asEyeByVector = self.__C )
elif self.isscalar():
return Covariance(self.__name+"T", asEyeByScalar = self.__C )
- elif self.isobject():
+ elif self.isobject() and hasattr(self.__C,"getT"):
return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
+ else:
+ raise AttributeError("the %s covariance matrix has no getT attribute."%(self.__name,))
def cholesky(self):
"Décomposition de Cholesky"
return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
elif self.isobject() and hasattr(self.__C,"cholesky"):
return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
+ else:
+ raise AttributeError("the %s covariance matrix has no cholesky attribute."%(self.__name,))
def choleskyI(self):
"Inversion de la décomposition de Cholesky"
if self.ismatrix():
- return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
+ return Covariance(self.__name+"H", asCovariance = numpy.linalg.inv(numpy.linalg.cholesky(self.__C)) )
elif self.isvector():
return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
elif self.isscalar():
return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
elif self.isobject() and hasattr(self.__C,"choleskyI"):
return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
+ else:
+ raise AttributeError("the %s covariance matrix has no choleskyI attribute."%(self.__name,))
def sqrtm(self):
"Racine carrée matricielle"
if self.ismatrix():
import scipy
- return Covariance(self.__name+"C", asCovariance = scipy.linalg.sqrtm(self.__C) )
+ return Covariance(self.__name+"C", asCovariance = numpy.real(scipy.linalg.sqrtm(self.__C)) )
elif self.isvector():
return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
elif self.isscalar():
return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
- elif self.isobject() and hasattr(self.__C,"sqrt"):
- return Covariance(self.__name+"C", asCovObject = self.__C.sqrt() )
+ elif self.isobject() and hasattr(self.__C,"sqrtm"):
+ return Covariance(self.__name+"C", asCovObject = self.__C.sqrtm() )
+ else:
+ raise AttributeError("the %s covariance matrix has no sqrtm attribute."%(self.__name,))
def sqrtmI(self):
"Inversion de la racine carrée matricielle"
if self.ismatrix():
import scipy
- return Covariance(self.__name+"H", asCovariance = scipy.linalg.sqrtm(self.__C).I )
+ return Covariance(self.__name+"H", asCovariance = numpy.linalg.inv(numpy.real(scipy.linalg.sqrtm(self.__C))) )
elif self.isvector():
return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
elif self.isscalar():
return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
- elif self.isobject() and hasattr(self.__C,"sqrtI"):
- return Covariance(self.__name+"H", asCovObject = self.__C.sqrtI() )
+ elif self.isobject() and hasattr(self.__C,"sqrtmI"):
+ return Covariance(self.__name+"H", asCovObject = self.__C.sqrtmI() )
+ else:
+ raise AttributeError("the %s covariance matrix has no sqrtmI attribute."%(self.__name,))
def diag(self, msize=None):
"Diagonale de la matrice"
raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
else:
return self.__C * numpy.ones(int(msize))
- elif self.isobject():
+ elif self.isobject() and hasattr(self.__C,"diag"):
return self.__C.diag()
+ else:
+ raise AttributeError("the %s covariance matrix has no diag attribute."%(self.__name,))
def asfullmatrix(self, msize=None):
"Matrice pleine"
return numpy.asarray( self.__C * numpy.eye(int(msize)), float )
elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
return self.__C.asfullmatrix()
+ else:
+ raise AttributeError("the %s covariance matrix has no asfullmatrix attribute."%(self.__name,))
def trace(self, msize=None):
"Trace de la matrice"
return self.__C * int(msize)
elif self.isobject():
return self.__C.trace()
+ else:
+ raise AttributeError("the %s covariance matrix has no trace attribute."%(self.__name,))
def getO(self):
return self
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.identity(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
#
- _bgcenter = numpy.ravel(_bgcenter)[:,None]
+ _bgcenter = numpy.ravel(_bgcenter).reshape((-1,1))
if _nbmembers < 1:
raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
if _bgcovariance is None:
return BackgroundEnsemble
# ==============================================================================
-def EnsembleOfAnomalies( _ensemble, _optmean = None):
+def EnsembleOfAnomalies( Ensemble, OptMean = None, Normalisation = 1.):
"Renvoie les anomalies centrées à partir d'un ensemble TailleEtat*NbMembres"
- if _optmean is None:
- Em = numpy.asarray(_ensemble).mean(axis=1, dtype=mfp).astype('float')[:,numpy.newaxis]
+ if OptMean is None:
+ __Em = numpy.asarray(Ensemble).mean(axis=1, dtype=mfp).astype('float').reshape((-1,1))
else:
- Em = numpy.ravel(_optmean)[:,numpy.newaxis]
+ __Em = numpy.ravel(OptMean).reshape((-1,1))
#
- return numpy.asarray(_ensemble) - Em
+ return Normalisation * (numpy.asarray(Ensemble) - __Em)
+
+# ==============================================================================
+def EnsembleErrorCovariance( Ensemble ):
+ "Renvoie la covariance d'ensemble"
+ __Anomalies = EnsembleOfAnomalies( Ensemble )
+ __n, __m = numpy.asarray(__Anomalies).shape
+ __Covariance = (__Anomalies @ __Anomalies.T) / (__m-1)
+ # Assure la symétrie
+ __Covariance = (__Covariance + __Covariance.T) * 0.5
+ # Assure la positivité
+ __epsilon = mpr*numpy.trace(__Covariance)
+ __Covariance = __Covariance + __epsilon * numpy.identity(__n)
+ #
+ return __Covariance
# ==============================================================================
def CovarianceInflation(
__n, __m = numpy.asarray(InputCovOrEns).shape
if __n != __m:
raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
- OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * numpy.eye(__n)
+ OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * numpy.identity(__n)
#
elif InflationType == "HybridOnBackgroundCovariance":
if InflationFactor < 0.:
EpY = EnsembleOfCenteredPerturbations(Ynpu, Rn, __m)
EpYm = EpY.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
#
- EaX = EnsembleOfAnomalies( Xn_predicted ) / numpy.sqrt(__m-1)
- EaY = (HX_predicted - Hfm - EpY + EpYm) / numpy.sqrt(__m-1)
+ EaX = EnsembleOfAnomalies( Xn_predicted ) / math.sqrt(__m-1)
+ EaY = (HX_predicted - Hfm - EpY + EpYm) / math.sqrt(__m-1)
#
Kn = EaX @ EaY.T @ numpy.linalg.inv( EaY @ EaY.T)
#
if selfA._toStore("CostFunctionJAtCurrentOptimum"):
selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
if selfA._toStore("APosterioriCovariance"):
- Eai = EnsembleOfAnomalies( Xn ) / numpy.sqrt(__m-1) # Anomalies
- Pn = Eai @ Eai.T
- Pn = 0.5 * (Pn + Pn.T)
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
+ selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
if selfA._parameters["EstimationOf"] == "Parameters" \
and J < previousJMinimum:
previousJMinimum = J
elif VariantM != "KalmanFilterFormula":
RI = R.getI()
if VariantM == "KalmanFilterFormula":
- RIdemi = R.choleskyI()
+ RIdemi = R.sqrtmI()
#
# Initialisation
# --------------
if hasattr(Q,"asfullmatrix"): Qn = Q.asfullmatrix(__n)
else: Qn = Q
Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
+ #~ Xn = EnsembleOfBackgroundPerturbations( Xb, Pn, __m )
#
if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
selfA.StoredVariables["Analysis"].store( Xb )
#
#--------------------------
if VariantM == "KalmanFilterFormula":
- mS = RIdemi * EaHX / numpy.sqrt(__m-1)
+ mS = RIdemi * EaHX / math.sqrt(__m-1)
delta = RIdemi * ( Ynpu - Hfm )
- mT = numpy.linalg.inv( numpy.eye(__m) + mS.T @ mS )
- vw = mT @ mS.transpose() @ delta
+ mT = numpy.linalg.inv( numpy.identity(__m) + mS.T @ mS )
+ vw = mT @ mS.T @ delta
#
Tdemi = numpy.real(scipy.linalg.sqrtm(mT))
- mU = numpy.eye(__m)
+ mU = numpy.identity(__m)
#
- EaX = EaX / numpy.sqrt(__m-1)
- Xn = Xfm + EaX @ ( vw.reshape((__m,1)) + numpy.sqrt(__m-1) * Tdemi @ mU )
+ EaX = EaX / math.sqrt(__m-1)
+ Xn = Xfm + EaX @ ( vw.reshape((__m,1)) + math.sqrt(__m-1) * Tdemi @ mU )
#--------------------------
elif VariantM == "Variational":
HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
)
#
Hto = EaHX.T @ (RI * EaHX)
- Htb = (__m-1) * numpy.eye(__m)
+ Htb = (__m-1) * numpy.identity(__m)
Hta = Hto + Htb
#
Pta = numpy.linalg.inv( Hta )
#
Hto = EaHX.T @ (RI * EaHX)
Htb = __m * \
- ( (1 + 1/__m + vw.T @ vw) * numpy.eye(__m) - 2 * vw @ vw.T ) \
+ ( (1 + 1/__m + vw.T @ vw) * numpy.identity(__m) - 2 * vw @ vw.T ) \
/ (1 + 1/__m + vw.T @ vw)**2
Hta = Hto + Htb
#
#
Hto = EaHX.T @ (RI * EaHX)
Htb = (__m+1) * \
- ( (1 + 1/__m + vw.T @ vw) * numpy.eye(__m) - 2 * vw @ vw.T ) \
+ ( (1 + 1/__m + vw.T @ vw) * numpy.identity(__m) - 2 * vw @ vw.T ) \
/ (1 + 1/__m + vw.T @ vw)**2
Hta = Hto + Htb
#
#
Hto = EaHX.T @ (RI * EaHX)
Htb = ((__m+1) / (__m-1)) * \
- ( (1 + 1/__m + vw.T @ vw / (__m-1)) * numpy.eye(__m) - 2 * vw @ vw.T / (__m-1) ) \
+ ( (1 + 1/__m + vw.T @ vw / (__m-1)) * numpy.identity(__m) - 2 * vw @ vw.T / (__m-1) ) \
/ (1 + 1/__m + vw.T @ vw / (__m-1))**2
Hta = Hto + Htb
#
if selfA._toStore("CostFunctionJAtCurrentOptimum"):
selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
if selfA._toStore("APosterioriCovariance"):
- Eai = EnsembleOfAnomalies( Xn ) / numpy.sqrt(__m-1) # Anomalies
- Pn = Eai @ Eai.T
- Pn = 0.5 * (Pn + Pn.T)
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
+ selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
if selfA._parameters["EstimationOf"] == "Parameters" \
and J < previousJMinimum:
previousJMinimum = J
XaMin = Xa
if selfA._toStore("APosterioriCovariance"):
covarianceXaMin = Pn
+ # ---> Pour les smoothers
+ if selfA._toStore("CurrentEnsembleState"):
+ selfA.StoredVariables["CurrentEnsembleState"].store( Xn )
+ if selfA._toStore("LastEnsembleForecastState"):
+ selfA.StoredVariables["LastEnsembleForecastState"].store( EMX )
#
# Stockage final supplémentaire de l'optimum en estimation de paramètres
# ----------------------------------------------------------------------
#--------------------------
if VariantM == "MLEF13":
Xfm = numpy.ravel(Xn_predicted.mean(axis=1, dtype=mfp).astype('float'))
- EaX = EnsembleOfAnomalies( Xn_predicted ) / numpy.sqrt(__m-1)
- Ua = numpy.eye(__m)
+ EaX = EnsembleOfAnomalies( Xn_predicted, Xfm, 1./math.sqrt(__m-1) )
+ Ua = numpy.identity(__m)
__j = 0
Deltaw = 1
if not BnotT:
- Ta = numpy.eye(__m)
+ Ta = numpy.identity(__m)
vw = numpy.zeros(__m)
while numpy.linalg.norm(Deltaw) >= _e and __j <= _jmax:
vx1 = (Xfm + EaX @ vw).reshape((__n,1))
if BnotT:
E1 = vx1 + _epsilon * EaX
else:
- E1 = vx1 + numpy.sqrt(__m-1) * EaX @ Ta
+ E1 = vx1 + math.sqrt(__m-1) * EaX @ Ta
#
HE2 = H( [(E1[:,i,numpy.newaxis], Un) for i in range(__m)],
argsAsSerie = True,
if BnotT:
EaY = (HE2 - vy2) / _epsilon
else:
- EaY = ( (HE2 - vy2) @ numpy.linalg.inv(Ta) ) / numpy.sqrt(__m-1)
+ EaY = ( (HE2 - vy2) @ numpy.linalg.inv(Ta) ) / math.sqrt(__m-1)
#
GradJ = numpy.ravel(vw[:,None] - EaY.transpose() @ (RI * ( Ynpu - vy2 )))
- mH = numpy.eye(__m) + EaY.transpose() @ (RI * EaY)
+ mH = numpy.identity(__m) + EaY.transpose() @ (RI * EaY)
Deltaw = - numpy.linalg.solve(mH,GradJ)
#
vw = vw + Deltaw
if BnotT:
Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
#
- Xn = vx1 + numpy.sqrt(__m-1) * EaX @ Ta @ Ua
+ Xn = vx1 + math.sqrt(__m-1) * EaX @ Ta @ Ua
#--------------------------
else:
raise ValueError("VariantM has to be chosen in the authorized methods list.")
if selfA._toStore("CostFunctionJAtCurrentOptimum"):
selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
if selfA._toStore("APosterioriCovariance"):
- Eai = EnsembleOfAnomalies( Xn ) / numpy.sqrt(__m-1) # Anomalies
- Pn = Eai @ Eai.T
- Pn = 0.5 * (Pn + Pn.T)
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
+ selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
if selfA._parameters["EstimationOf"] == "Parameters" \
and J < previousJMinimum:
previousJMinimum = J
#--------------------------
if VariantM == "IEnKF12":
Xfm = numpy.ravel(Xn.mean(axis=1, dtype=mfp).astype('float'))
- EaX = EnsembleOfAnomalies( Xn ) / numpy.sqrt(__m-1)
+ EaX = EnsembleOfAnomalies( Xn ) / math.sqrt(__m-1)
__j = 0
Deltaw = 1
if not BnotT:
- Ta = numpy.eye(__m)
+ Ta = numpy.identity(__m)
vw = numpy.zeros(__m)
while numpy.linalg.norm(Deltaw) >= _e and __j <= _jmax:
vx1 = (Xfm + EaX @ vw).reshape((__n,1))
if BnotT:
E1 = vx1 + _epsilon * EaX
else:
- E1 = vx1 + numpy.sqrt(__m-1) * EaX @ Ta
+ E1 = vx1 + math.sqrt(__m-1) * EaX @ Ta
#
if selfA._parameters["EstimationOf"] == "State": # Forecast + Q
E2 = M( [(E1[:,i,numpy.newaxis], Un) for i in range(__m)],
if BnotT:
EaY = (HE2 - vy2) / _epsilon
else:
- EaY = ( (HE2 - vy2) @ numpy.linalg.inv(Ta) ) / numpy.sqrt(__m-1)
+ EaY = ( (HE2 - vy2) @ numpy.linalg.inv(Ta) ) / math.sqrt(__m-1)
#
GradJ = numpy.ravel(vw[:,None] - EaY.transpose() @ (RI * ( Ynpu - vy1 )))
- mH = numpy.eye(__m) + EaY.transpose() @ (RI * EaY)
+ mH = numpy.identity(__m) + EaY.transpose() @ (RI * EaY)
Deltaw = - numpy.linalg.solve(mH,GradJ)
#
vw = vw + Deltaw
#
if BnotT:
Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
- A2 = numpy.sqrt(__m-1) * A2 @ Ta / _epsilon
+ A2 = math.sqrt(__m-1) * A2 @ Ta / _epsilon
#
Xn = vx2 + A2
#--------------------------
if selfA._toStore("CostFunctionJAtCurrentOptimum"):
selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
if selfA._toStore("APosterioriCovariance"):
- Eai = EnsembleOfAnomalies( Xn ) / numpy.sqrt(__m-1) # Anomalies
- Pn = Eai @ Eai.T
- Pn = 0.5 * (Pn + Pn.T)
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
+ selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
if selfA._parameters["EstimationOf"] == "Parameters" \
and J < previousJMinimum:
previousJMinimum = J
