From d4f0a8f335b824294ea92f5064c43bf9b7a010dc Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Wed, 30 Dec 2020 18:28:03 +0100 Subject: [PATCH] Improvement and extension of EnKF algorithm (KFF) --- .../daAlgorithms/EnsembleKalmanFilter.py | 10 ++++-- src/daComposant/daCore/BasicObjects.py | 24 +++++++++++++ src/daComposant/daCore/NumericObjects.py | 35 +++++++++++-------- 3 files changed, 52 insertions(+), 17 deletions(-) diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py index 4a3e6ad..d242b6b 100644 --- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py +++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py @@ -33,7 +33,11 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): default = "StochasticEnKF", typecast = str, message = "Minimiseur utilisé", - listval = ["StochasticEnKF", "DeterministicEnKF", "ETKF"], + listval = [ + "StochasticEnKF", + "ETKF", + "ETKF-KFF", + ], ) self.defineRequiredParameter( name = "NumberOfMembers", @@ -152,8 +156,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): # if self._parameters["Minimizer"] == "StochasticEnKF": NumericObjects.senkf(self, Xb, Y, U, HO, EM, CM, R, B, Q) - elif self._parameters["Minimizer"] in ["DeterministicEnKF", "ETKF"]: - NumericObjects.etkf(self, Xb, Y, U, HO, EM, CM, R, B, Q) + elif self._parameters["Minimizer"] in ["ETKF", "ETKF-KFF"]: + NumericObjects.etkf(self, Xb, Y, U, HO, EM, CM, R, B, Q, KorV="KalmanFilterFormula") else: raise ValueError("Error in Minimizer name: %s"%self._parameters["Minimizer"]) # diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index d3b060f..429c2bc 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -1826,6 +1826,30 @@ class Covariance(object): elif self.isobject() and hasattr(self.__C,"choleskyI"): return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() ) + def sqrtm(self): + "Racine carrée matricielle" + if self.ismatrix(): + import scipy + return Covariance(self.__name+"C", asCovariance = 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() ) + + 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 ) + 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() ) + def diag(self, msize=None): "Diagonale de la matrice" if self.ismatrix(): diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index 4ff8cb5..15300a9 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -791,7 +791,7 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): return 0 # ============================================================================== -def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): +def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, KorV="KalmanFilterFormula"): """ Ensemble-Transform EnKF (ETKF or Deterministic EnKF: Bishop 2001, Hunt 2007) @@ -832,7 +832,6 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): or selfA._toStore("APosterioriCovariance"): BI = B.getI() RI = R.getI() - RIdemi = R.choleskyI() # # Initialisation # -------------- @@ -898,18 +897,25 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float') Hfm = HX_predicted.mean(axis=1, dtype=mfp).astype('float') # - EaX = (Xn_predicted - Xfm.reshape((__n,-1))) / numpy.sqrt(__m-1) - EaHX = (HX_predicted - Hfm.reshape((__p,-1))) / numpy.sqrt(__m-1) - # - mS = RIdemi * EaHX - delta = RIdemi * ( Ynpu.reshape((__p,-1)) - Hfm.reshape((__p,-1)) ) - mT = numpy.linalg.inv( numpy.eye(__m) + mS.T @ mS ) - vw = mT @ mS.transpose() @ delta - # - Tdemi = numpy.linalg.cholesky(mT) - mU = numpy.eye(__m) - # - Xn = Xfm.reshape((__n,-1)) + EaX @ ( vw.reshape((__m,-1)) + numpy.sqrt(__m-1) * Tdemi @ mU ) + EaX = (Xn_predicted - Xfm.reshape((__n,-1))) + EaHX = (HX_predicted - Hfm.reshape((__p,-1))) + # + #-------------------------- + if KorV == "KalmanFilterFormula": + EaX = EaX / numpy.sqrt(__m-1) + RIdemi = R.choleskyI() + mS = RIdemi * EaHX / numpy.sqrt(__m-1) + delta = RIdemi * ( Ynpu.reshape((__p,-1)) - Hfm.reshape((__p,-1)) ) + mT = numpy.linalg.inv( numpy.eye(__m) + mS.T @ mS ) + vw = mT @ mS.transpose() @ delta + # + Tdemi = numpy.real(scipy.linalg.sqrtm(mT)) + mU = numpy.eye(__m) + # + Xn = Xfm.reshape((__n,-1)) + EaX @ ( vw.reshape((__m,-1)) + numpy.sqrt(__m-1) * Tdemi @ mU ) + #-------------------------- + else: + raise ValueError("KorV has to be chosen as either \"KalmanFilterFormula\" or \"Variational\".") # if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies": Xn = CovarianceInflation( Xn, @@ -918,6 +924,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): ) # Xa = Xn.mean(axis=1, dtype=mfp).astype('float') + #-------------------------- # if selfA._parameters["StoreInternalVariables"] \ or selfA._toStore("CostFunctionJ") \ -- 2.39.2