From: Jean-Philippe ARGAUD Date: Fri, 26 Mar 2021 21:05:08 +0000 (+0100) Subject: Minor documentation, improvements and fixes for internal variables X-Git-Tag: V9_7_0b1~11 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=fb340ea27ba3b6cbc6e581eda73ea9c460244a4f;p=modules%2Fadao.git Minor documentation, improvements and fixes for internal variables --- diff --git a/doc/en/ref_algorithm_EnsembleKalmanFilter.rst b/doc/en/ref_algorithm_EnsembleKalmanFilter.rst index 5f7efef..78d3b74 100644 --- a/doc/en/ref_algorithm_EnsembleKalmanFilter.rst +++ b/doc/en/ref_algorithm_EnsembleKalmanFilter.rst @@ -91,6 +91,10 @@ Without being a universal recommandation, one recommend to use EnKF as a referen .. include:: snippets/EstimationOf.rst +.. include:: snippets/InflationFactor.rst + +.. include:: snippets/InflationType.rst + .. include:: snippets/NumberOfMembers.rst .. include:: snippets/SetSeed.rst diff --git a/doc/en/snippets/InflationFactor.rst b/doc/en/snippets/InflationFactor.rst new file mode 100644 index 0000000..e6009f9 --- /dev/null +++ b/doc/en/snippets/InflationFactor.rst @@ -0,0 +1,12 @@ +.. 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.}`` diff --git a/doc/en/snippets/InflationType.rst b/doc/en/snippets/InflationType.rst new file mode 100644 index 0000000..4cbd8be --- /dev/null +++ b/doc/en/snippets/InflationType.rst @@ -0,0 +1,17 @@ +.. 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"}`` diff --git a/doc/fr/ref_algorithm_EnsembleKalmanFilter.rst b/doc/fr/ref_algorithm_EnsembleKalmanFilter.rst index f981840..467f589 100644 --- a/doc/fr/ref_algorithm_EnsembleKalmanFilter.rst +++ b/doc/fr/ref_algorithm_EnsembleKalmanFilter.rst @@ -92,6 +92,10 @@ Sans pouvoir prétendre à l'universalité, on recommande d'utiliser l'EnKF comm .. include:: snippets/EstimationOf.rst +.. include:: snippets/InflationFactor.rst + +.. include:: snippets/InflationType.rst + .. include:: snippets/NumberOfMembers.rst .. include:: snippets/SetSeed.rst diff --git a/doc/fr/snippets/ConstrainedBy.rst b/doc/fr/snippets/ConstrainedBy.rst index f13a6c4..5df59e0 100644 --- a/doc/fr/snippets/ConstrainedBy.rst +++ b/doc/fr/snippets/ConstrainedBy.rst @@ -3,7 +3,7 @@ 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"}`` diff --git a/doc/fr/snippets/InflationFactor.rst b/doc/fr/snippets/InflationFactor.rst new file mode 100644 index 0000000..2253319 --- /dev/null +++ b/doc/fr/snippets/InflationFactor.rst @@ -0,0 +1,12 @@ +.. 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.}`` diff --git a/doc/fr/snippets/InflationType.rst b/doc/fr/snippets/InflationType.rst new file mode 100644 index 0000000..a3923f1 --- /dev/null +++ b/doc/fr/snippets/InflationType.rst @@ -0,0 +1,17 @@ +.. 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"}`` diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py index 5eb7d5f..cd43a13 100644 --- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py +++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py @@ -72,17 +72,20 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): ) 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( @@ -93,29 +96,18 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): 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", @@ -156,7 +148,11 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): "SimulatedObservationAtCurrentAnalysis", "SimulatedObservationAtCurrentOptimum", "SimulatedObservationAtCurrentState", - ] + ], + listadv = [ + "CurrentEnsembleState", + "LastEnsembleForecastState", + ], ) self.requireInputArguments( mandatory= ("Xb", "Y", "HO", "R", "B"), diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index 27d7695..c5f2cb0 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -661,7 +661,8 @@ class Algorithm(object): 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") @@ -676,6 +677,7 @@ class Algorithm(object): 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") @@ -1835,12 +1837,12 @@ class Covariance(object): 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 @@ -1853,8 +1855,10 @@ class Covariance(object): 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" @@ -1866,41 +1870,49 @@ class Covariance(object): 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" @@ -1913,8 +1925,10 @@ class Covariance(object): 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" @@ -1929,6 +1943,8 @@ class Covariance(object): 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" @@ -1943,6 +1959,8 @@ class Covariance(object): 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 diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index 30c05e5..2468049 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -546,14 +546,14 @@ def EnsembleOfBackgroundPerturbations( _bgcenter, _bgcovariance, _nbmembers, _wi 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: @@ -582,14 +582,28 @@ def EnsembleOfBackgroundPerturbations( _bgcenter, _bgcovariance, _nbmembers, _wi 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( @@ -632,7 +646,7 @@ 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.: @@ -2181,8 +2195,8 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): 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) # @@ -2266,10 +2280,7 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): 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 @@ -2331,7 +2342,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): elif VariantM != "KalmanFilterFormula": RI = R.getI() if VariantM == "KalmanFilterFormula": - RIdemi = R.choleskyI() + RIdemi = R.sqrtmI() # # Initialisation # -------------- @@ -2344,6 +2355,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): 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 ) @@ -2404,16 +2416,16 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): # #-------------------------- 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 @@ -2438,7 +2450,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): ) # 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 ) @@ -2470,7 +2482,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): # 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 # @@ -2503,7 +2515,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): # 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 # @@ -2536,7 +2548,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): # 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 # @@ -2622,16 +2634,18 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): 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 # ---------------------------------------------------------------------- @@ -2744,12 +2758,12 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", #-------------------------- 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)) @@ -2757,7 +2771,7 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", 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, @@ -2767,10 +2781,10 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", 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 @@ -2783,7 +2797,7 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", 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.") @@ -2862,10 +2876,7 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", 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 @@ -2971,11 +2982,11 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", #-------------------------- 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)) @@ -2983,7 +2994,7 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", 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)], @@ -3003,10 +3014,10 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", 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 @@ -3020,7 +3031,7 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", # 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 #-------------------------- @@ -3101,10 +3112,7 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", 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