Improvement and extension of EnKF algorithm (KFF)

diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
index 4a3e6ad658f79419793a5e796f14163b45c01df9..d242b6b6017a508110ddaba44e40fb74a68dd8c4 100644 (file)
--- 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 d3b060fa44357cdfca3d90726f38749ce4d0f630..429c2bc957fa3635cca588d06d476265bcb521bd 100644 (file)
--- 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 4ff8cb553b43917fe936bda35f465e30975b1441..15300a9704989247adc111518786507c4e0a1b19 100644 (file)
--- 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") \