Improvement and extension of EnKF algorithm (loc)

diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
index c24bc9400226e917e17cf1c8c00764e60c236460..4a3e6ad658f79419793a5e796f14163b45c01df9 100644 (file)
--- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
+++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py
@@ -49,6 +49,51 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             message  = "Estimation d'etat ou de parametres",
             listval  = ["State", "Parameters"],
             )
+        self.defineRequiredParameter(
+            name     = "InflationType",
+            default  = "MultiplicativeOnAnalysisCovariance",
+            typecast = str,
+            message  = "Méthode d'inflation d'ensemble",
+            listval  = [
+                "MultiplicativeOnAnalysisCovariance",
+                "MultiplicativeOnBackgroundCovariance",
+                "MultiplicativeOnAnalysisAnomalies",
+                "MultiplicativeOnBackgroundAnomalies",
+                "AdditiveOnBackgroundCovariance",
+                "AdditiveOnAnalysisCovariance",
+                "HybridOnBackgroundCovariance",
+                ],
+            )
+        self.defineRequiredParameter(
+            name     = "InflationFactor",
+            default  = 1.,
+            typecast = float,
+            message  = "Facteur d'inflation",
+            minval   = 0.,
+            )
+        self.defineRequiredParameter(
+            name     = "LocalizationType",
+            default  = "SchurLocalization",
+            typecast = str,
+            message  = "Méthode d'inflation d'ensemble",
+            listval  = [
+                "CovarianceLocalization",
+                "DomainLocalization",
+                "SchurLocalization",
+                "GaspariCohnLocalization",
+                ],
+            )
+        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",
+            )
         self.defineRequiredParameter(
             name     = "SetSeed",
             typecast = numpy.random.seed,

diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py
index 8acdd77e8332d00420f7e64016d26140427b2e21..4ff8cb553b43917fe936bda35f465e30975b1441 100644 (file)
--- a/src/daComposant/daCore/NumericObjects.py
+++ b/src/daComposant/daCore/NumericObjects.py
@@ -507,6 +507,71 @@ def mmqr(
     #
     return variables, Ecart, [n,p,iteration,increment,0]
 
+# ==============================================================================
+def CovarianceInflation(
+        InputCovOrEns,
+        InflationType   = None,
+        InflationFactor = None,
+        BackgroundCov   = None,
+        ):
+    """
+    Inflation applicable soit sur Pb ou Pa, soit sur les ensembles EXb ou EXa
+
+    Synthèse : Hunt 2007, section 2.3.5
+    """
+    if InflationFactor is None:
+        return InputCovOrEns
+    else:
+        InflationFactor = float(InflationFactor)
+    #
+    if InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
+        if InflationFactor < 1.:
+            raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
+        if InflationFactor < 1.+mpr:
+            return InputCovOrEns
+        OutputCovOrEns = InflationFactor**2 * InputCovOrEns
+    #
+    elif InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
+        if InflationFactor < 1.:
+            raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
+        if InflationFactor < 1.+mpr:
+            return InputCovOrEns
+        InputCovOrEnsMean = InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
+        OutputCovOrEns = InputCovOrEnsMean[:,numpy.newaxis] \
+            + InflationFactor * (InputCovOrEns - InputCovOrEnsMean[:,numpy.newaxis])
+    #
+    elif InflationType in ["AdditiveOnBackgroundCovariance", "AdditiveOnAnalysisCovariance"]:
+        if InflationFactor < 0.:
+            raise ValueError("Inflation factor for additive inflation has to be greater or equal than 0.")
+        if InflationFactor < mpr:
+            return InputCovOrEns
+        __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)
+    #
+    elif InflationType == "HybridOnBackgroundCovariance":
+        if InflationFactor < 0.:
+            raise ValueError("Inflation factor for hybrid inflation has to be greater or equal than 0.")
+        if InflationFactor < mpr:
+            return InputCovOrEns
+        __n, __m = numpy.asarray(InputCovOrEns).shape
+        if __n != __m:
+            raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
+        if BackgroundCov is None:
+            raise ValueError("Background covariance matrix B has to be given for hybrid inflation.")
+        if InputCovOrEns.shape != BackgroundCov.shape:
+            raise ValueError("Ensemble covariance matrix has to be of same size than background covariance matrix B.")
+        OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * BackgroundCov
+    #
+    elif InflationType == "Relaxation":
+        raise NotImplementedError("InflationType Relaxation")
+    #
+    else:
+        raise ValueError("Error in inflation type, '%s' is not a valid keyword."%InflationType)
+    #
+    return OutputCovOrEns
+
 # ==============================================================================
 def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
     """
@@ -590,6 +655,12 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
         else:
             Un = None
         #
+        if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
+            Xn = CovarianceInflation( Xn,
+                selfA._parameters["InflationType"],
+                selfA._parameters["InflationFactor"],
+                )
+        #
         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
             for i in range(__m):
                 qi = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__n), Qn, (1,1,1))).T
@@ -623,6 +694,12 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
             ri = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__p), Rn, (1,1,1))).T
             Xn[:,i] = Xn_predicted[:,i] + K * (Ynpu + ri - HX_predicted[:,i])
         #
+        if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
+            Xn = CovarianceInflation( Xn,
+                selfA._parameters["InflationType"],
+                selfA._parameters["InflationFactor"],
+                )
+        #
         Xa = Xn.mean(axis=1, dtype=mfp).astype('float')
         #
         if selfA._parameters["StoreInternalVariables"] \
@@ -797,6 +874,12 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
         else:
             Un = None
         #
+        if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
+            Xn = CovarianceInflation( Xn,
+                selfA._parameters["InflationType"],
+                selfA._parameters["InflationFactor"],
+                )
+        #
         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
             for i in range(__m):
                 qi = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__n), Qn, (1,1,1))).T
@@ -828,6 +911,12 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
         #
         Xn = Xfm.reshape((__n,-1)) + EaX @ ( vw.reshape((__m,-1)) + numpy.sqrt(__m-1) * Tdemi @ mU )
         #
+        if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
+            Xn = CovarianceInflation( Xn,
+                selfA._parameters["InflationType"],
+                selfA._parameters["InflationFactor"],
+                )
+        #
         Xa = Xn.mean(axis=1, dtype=mfp).astype('float')
         #
         if selfA._parameters["StoreInternalVariables"] \