Improvement and extension of EnKF algorithm (MLEF-T)

-------------

diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
index 313ea75ecddcbe9f21e3be7f03c4a8101ada2221..8a56f1fe0e33dce15cae54bc6966ad04cb558c5d 100644 (file)
--- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
+++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
@@ -44,6 +44,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
                 "ETKF-N-16",
                 "MLEF",
                 "MLEF-B",
+                "MLEF-T",
                 ],
             )
         self.defineRequiredParameter(
@@ -189,6 +190,9 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         elif self._parameters["Minimizer"] in ["MLEF-B", "MLEF"]:
             NumericObjects.mlef(self, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=False)
         #
+        elif self._parameters["Minimizer"] == "MLEF-T":
+            NumericObjects.mlef(self, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=True)
+        #
         #--------------------------
         else:
             raise ValueError("Error in Minimizer name: %s"%self._parameters["Minimizer"])

diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py
index 9d6ed56c9a9e6e7b9f6dcb95725a9415970e1dec..7103817264ffa61f5f097f3d65475750a7f9368e 100644 (file)
--- a/src/daComposant/daCore/NumericObjects.py
+++ b/src/daComposant/daCore/NumericObjects.py
@@ -1262,6 +1262,8 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=False, _epsilon=1.e-1, _e=1
             #
             if BnotT:
                 EE = vx.reshape((__n,-1)) + _epsilon * EaX # 7:
+            else:
+                EE = vx.reshape((__n,-1)) + numpy.sqrt(__m-1) * EaX @ Ta # 8:
             #
             EZ = H( [(EE[:,i], Un) for i in range(__m)],
                 argsAsSerie = True,
@@ -1271,6 +1273,8 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=False, _epsilon=1.e-1, _e=1
             #
             if BnotT:
                 EY = (EZ - ybar) / _epsilon # 11:
+            else:
+                EY = ( (EZ - ybar) @ numpy.linalg.inv(Ta) ) / numpy.sqrt(__m-1) # 12:
             #
             GradJ = numpy.ravel(vw.reshape((__m,1)) - EY.transpose() @ (RI * (Ynpu - ybar))) # 13:
             mH = numpy.eye(__m) + EY.transpose() @ (RI * EY) # 14: