From 57fa783892b3ab06fb2a82a49ee6e9937da03009 Mon Sep 17 00:00:00 2001
From: Jean-Philippe ARGAUD <jean-philippe.argaud@edf.fr>
Date: Wed, 6 Jan 2021 17:52:31 +0100
Subject: [PATCH] Improvement of EnKF algorithm and multifonction

---
 .../daAlgorithms/EnsembleKalmanFilter.py      |  8 ++++---
 src/daComposant/daCore/BasicObjects.py        | 22 +++++++++----------
 src/daComposant/daCore/NumericObjects.py      | 10 ++++++---
 3 files changed, 23 insertions(+), 17 deletions(-)

diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
index 9171340..0da6ba8 100644
--- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
+++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
@@ -35,7 +35,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             message  = "Minimiseur utilisé",
             listval  = [
                 "StochasticEnKF",
-                "ETKF",
+                "ETKF",     # Default ETKF
                 "ETKF-KFF",
                 "ETKF-VAR",
                 ],
@@ -45,7 +45,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             default  = 100,
             typecast = int,
             message  = "Nombre de membres dans l'ensemble",
-            minval   = -1,
+            minval   = 2,
             )
         self.defineRequiredParameter(
             name     = "EstimationOf",
@@ -157,8 +157,10 @@ 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 ["ETKF", "ETKF-KFF"]:
+        #
+        elif self._parameters["Minimizer"] in ["ETKF-KFF", "ETKF"]:
             NumericObjects.etkf(self, Xb, Y, U, HO, EM, CM, R, B, Q, KorV="KalmanFilterFormula")
+        #
         elif self._parameters["Minimizer"] == "ETKF-VAR":
             NumericObjects.etkf(self, Xb, Y, U, HO, EM, CM, R, B, Q, KorV="Variational")
         else:
diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py
index 7a97326..74f9bfb 100644
--- a/src/daComposant/daCore/BasicObjects.py
+++ b/src/daComposant/daCore/BasicObjects.py
@@ -266,28 +266,28 @@ class Operator(object):
         PlatformInfo.isIterable( _xuValue, True, " in Operator.appliedControledFormTo" )
         #
         if self.__Matrix is not None:
-            HxValue = []
+            _HxValue = []
             for paire in _xuValue:
                 _xValue, _uValue = paire
                 self.__addOneMatrixCall()
-                HxValue.append( self.__Matrix * _xValue )
+                _HxValue.append( self.__Matrix * _xValue )
         else:
-            HxValue = []
+            _HxValue = []
+            _xuArgs = []
             for paire in _xuValue:
-                _xuValue = []
                 _xValue, _uValue = paire
                 if _uValue is not None:
-                    _xuValue.append( paire )
+                    _xuArgs.append( paire )
                 else:
-                    _xuValue.append( _xValue )
-            self.__addOneMethodCall( len(_xuValue) )
+                    _xuArgs.append( _xValue )
+            self.__addOneMethodCall( len(_xuArgs) )
             if self.__extraArgs is None:
-                HxValue = self.__Method( _xuValue ) # Calcul MF
+                _HxValue = self.__Method( _xuArgs ) # Calcul MF
             else:
-                HxValue = self.__Method( _xuValue, self.__extraArgs ) # Calcul MF
+                _HxValue = self.__Method( _xuArgs, self.__extraArgs ) # Calcul MF
         #
-        if argsAsSerie: return HxValue
-        else:           return HxValue[-1]
+        if argsAsSerie: return _HxValue
+        else:           return _HxValue[-1]
 
     def appliedInXTo(self, paires, argsAsSerie = False):
         """
diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py
index 208e736..5f59d11 100644
--- a/src/daComposant/daCore/NumericObjects.py
+++ b/src/daComposant/daCore/NumericObjects.py
@@ -880,18 +880,22 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, KorV="KalmanFilterFormula"):
                 )
         #
         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
+            EMX = M( [(Xn[:,i], Un) for i in range(__m)], argsAsSerie = True )
             for i in range(__m):
                 qi = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__n), Qn, (1,1,1))).T
-                Xn_predicted[:,i] = numpy.asmatrix(numpy.ravel( M((Xn[:,i], Un)) )).T + qi
-                HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( H((Xn_predicted[:,i], Un)) )).T
+                Xn_predicted[:,i] = numpy.asmatrix(numpy.ravel( EMX[i] )).T + qi
+            EHX = H( [(Xn_predicted[:,i], Un) for i in range(__m)], argsAsSerie = True )
+            for i in range(__m):
+                HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( EHX[i] )).T
             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
                 Xn_predicted = Xn_predicted + Cm * Un
         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
             # --- > Par principe, M = Id, Q = 0
             Xn_predicted = Xn
+            EHX = H( [(Xn_predicted[:,i], Un) for i in range(__m)], argsAsSerie = True )
             for i in range(__m):
-                HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( H((Xn_predicted[:,i], Un)) )).T
+                HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( EHX[i] )).T
         #
         # Mean of forecast and observation of forecast
         Xfm  = Xn_predicted.mean(axis=1, dtype=mfp).astype('float')
-- 
2.39.2