Salome HOME
Documentation corrections for outputs
[modules/adao.git] / src / daComposant / daAlgorithms / KalmanFilter.py
index 9d6a0fdd2f3994d8ff7843b83afe74c8cb10d1f5..cc267e5a8965dbb8fa5ee2804dbb5e2957ae2d84 100644 (file)
@@ -1,6 +1,6 @@
 #-*-coding:iso-8859-1-*-
 #
-#  Copyright (C) 2008-2013 EDF R&D
+#  Copyright (C) 2008-2015 EDF R&D
 #
 #  This library is free software; you can redistribute it and/or
 #  modify it under the terms of the GNU Lesser General Public
@@ -21,8 +21,7 @@
 #  Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
 
 import logging
-from daCore import BasicObjects, PlatformInfo
-m = PlatformInfo.SystemUsage()
+from daCore import BasicObjects
 import numpy
 
 # ==============================================================================
@@ -30,14 +29,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
     def __init__(self):
         BasicObjects.Algorithm.__init__(self, "KALMANFILTER")
         self.defineRequiredParameter(
-            name     = "StoreSupplementaryCalculations",
-            default  = [],
-            typecast = tuple,
-            message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
-            listval  = ["APosterioriCovariance", "BMA", "Innovation"]
-            )
-        self.defineRequiredParameter(
-            name     = "EstimationType",
+            name     = "EstimationOf",
             default  = "State",
             typecast = str,
             message  = "Estimation d'etat ou de parametres",
@@ -49,16 +41,22 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             typecast = bool,
             message  = "Stockage des variables internes ou intermédiaires du calcul",
             )
+        self.defineRequiredParameter(
+            name     = "StoreSupplementaryCalculations",
+            default  = [],
+            typecast = tuple,
+            message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
+            listval  = ["APosterioriCovariance", "BMA", "CurrentState", "CostFunctionJ", "Innovation"]
+            )
 
     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
-        logging.debug("%s Lancement"%self._name)
-        logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
+        self._pre_run()
         #
         # Paramètres de pilotage
         # ----------------------
         self.setParameters(Parameters)
         #
-        if self._parameters["EstimationType"] == "Parameters":
+        if self._parameters["EstimationOf"] == "Parameters":
             self._parameters["StoreInternalVariables"] = True
         #
         # Opérateurs
@@ -71,7 +69,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         Ht = HO["Tangent"].asMatrix(Xb)
         Ha = HO["Adjoint"].asMatrix(Xb)
         #
-        if self._parameters["EstimationType"] == "State":
+        if self._parameters["EstimationOf"] == "State":
             Mt = EM["Tangent"].asMatrix(Xb)
             Ma = EM["Adjoint"].asMatrix(Xb)
         #
@@ -80,8 +78,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         else:
             Cm = None
         #
-        # Nombre de pas du Kalman identique au nombre de pas d'observations
-        # -----------------------------------------------------------------
+        # Nombre de pas identique au nombre de pas d'observations
+        # -------------------------------------------------------
         if hasattr(Y,"stepnumber"):
             duration = Y.stepnumber()
         else:
@@ -90,15 +88,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         # Précalcul des inversions de B et R
         # ----------------------------------
         if self._parameters["StoreInternalVariables"]:
-            if B is not None:
-                BI = B.I
-            elif self._parameters["B_scalar"] is not None:
-                BI = 1.0 / self._parameters["B_scalar"]
-            #
-            if R is not None:
-                RI = R.I
-            elif self._parameters["R_scalar"] is not None:
-                RI = 1.0 / self._parameters["R_scalar"]
+            BI = B.getI()
+            RI = R.getI()
         #
         # Initialisation
         # --------------
@@ -128,26 +119,27 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             else:
                 Un = None
             #
-            if self._parameters["EstimationType"] == "State":
+            if self._parameters["EstimationOf"] == "State":
                 Xn_predicted = Mt * Xn
                 if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
+                    Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape
                     Xn_predicted = Xn_predicted + Cm * Un
-                Pn_predicted = Mt * Pn * Ma + Q
-            elif self._parameters["EstimationType"] == "Parameters":
+                Pn_predicted = Q + Mt * Pn * Ma
+            elif self._parameters["EstimationOf"] == "Parameters":
                 # --- > Par principe, M = Id, Q = 0
                 Xn_predicted = Xn
                 Pn_predicted = Pn
             #
-            if self._parameters["EstimationType"] == "State":
+            if self._parameters["EstimationOf"] == "State":
                 d  = Ynpu - Ht * Xn_predicted
-            elif self._parameters["EstimationType"] == "Parameters":
+            elif self._parameters["EstimationOf"] == "Parameters":
                 d  = Ynpu - Ht * Xn_predicted
                 if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
                     d = d - Cm * Un
             #
-            K  = Pn_predicted * Ha * (Ht * Pn_predicted * Ha + R).I
-            Xn = Xn_predicted + K * d
-            Pn = Pn_predicted - K * Ht * Pn_predicted
+            Kn = Pn_predicted * Ha * (R + Ht * Pn_predicted * Ha).I
+            Xn = Xn_predicted + Kn * d
+            Pn = Pn_predicted - Kn * Ht * Pn_predicted
             #
             self.StoredVariables["Analysis"].store( Xn.A1 )
             if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
@@ -158,7 +150,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
                 Jb  = 0.5 * (Xn - Xb).T * BI * (Xn - Xb)
                 Jo  = 0.5 * d.T * RI * d
                 J   = float( Jb ) + float( Jo )
-                self.StoredVariables["CurrentState"].store( Xn.A1 )
+                if self._parameters["StoreInternalVariables"] or "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
+                    self.StoredVariables["CurrentState"].store( Xn )
                 self.StoredVariables["CostFunctionJb"].store( Jb )
                 self.StoredVariables["CostFunctionJo"].store( Jo )
                 self.StoredVariables["CostFunctionJ" ].store( J )
@@ -173,7 +166,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Stockage supplementaire de l'optimum en estimation de parametres
         # ----------------------------------------------------------------
-        if self._parameters["EstimationType"] == "Parameters":
+        if self._parameters["EstimationOf"] == "Parameters":
             self.StoredVariables["Analysis"].store( Xa.A1 )
             if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
                 self.StoredVariables["APosterioriCovariance"].store( covarianceXa )
@@ -181,9 +174,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         if "BMA" in self._parameters["StoreSupplementaryCalculations"]:
             self.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
         #
-        logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
-        logging.debug("%s Terminé"%self._name)
-        #
+        self._post_run(HO)
         return 0
 
 # ==============================================================================