Setting a posteriori covariance calculations optionnal

diff --git a/src/daComposant/daAlgorithms/3DVAR.py b/src/daComposant/daAlgorithms/3DVAR.py
index cee5c22b553ee25c9aec84af73b8993b8e873fe5..6210227822e4cbe018a6386b3d29f0ea533284cb 100644 (file)
--- a/src/daComposant/daAlgorithms/3DVAR.py
+++ b/src/daComposant/daAlgorithms/3DVAR.py
@@ -1,6 +1,6 @@
 #-*-coding:iso-8859-1-*-
 #
-#  Copyright (C) 2008-2011  EDF R&D
+#  Copyright (C) 2008-2012 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
@@ -127,7 +127,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Paramètres de pilotage
         # ----------------------
-        # Potentiels : "Bounds", "Minimizer", "MaximumNumberOfSteps", "ProjectedGradientTolerance", "GradientNormTolerance", "InnerMinimizer"
+        # Potentiels : "Bounds", "Minimizer", "MaximumNumberOfSteps", "ProjectedGradientTolerance", "GradientNormTolerance", "InnerMinimizer", "CalculateAPosterioriCovariance"
         if Parameters.has_key("Bounds") and (type(Parameters["Bounds"]) is type([]) or type(Parameters["Bounds"]) is type(())) and (len(Parameters["Bounds"]) > 0):
             Bounds = Parameters["Bounds"]
         else:
@@ -167,6 +167,11 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         else:
             InnerMinimizer = "BFGS"
         logging.debug("%s Minimiseur interne utilisé = %s"%(self._name, InnerMinimizer))
+        if Parameters.has_key("CalculateAPosterioriCovariance"):
+            CalculateAPosterioriCovariance = bool(Parameters["CalculateAPosterioriCovariance"])
+        else:
+            CalculateAPosterioriCovariance = False
+        logging.debug("%s Calcul de la covariance a posteriori = %s"%(self._name, CalculateAPosterioriCovariance))
         #
         # Minimisation de la fonctionnelle
         # --------------------------------
@@ -252,6 +257,19 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         self.StoredVariables["Analysis"].store( Xa.A1 )
         self.StoredVariables["Innovation"].store( d.A1 )
         #
+        # Calcul de la covariance d'analyse
+        # ---------------------------------
+        if CalculateAPosterioriCovariance:
+            Hessienne = []
+            nb = len(Xini)
+            for i in range(nb):
+                ee = numpy.matrix(numpy.zeros(nb)).T
+                ee[i] = 1.
+                Hessienne.append( ( BI*ee + Ht((Xa,RI*Hm(ee))) ).A1 )
+            Hessienne = numpy.matrix( Hessienne )
+            A = Hessienne.I
+            self.StoredVariables["APosterioriCovariance"].store( A )
+        #
         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("MB")))
         logging.debug("%s Terminé"%self._name)
         #

diff --git a/src/daComposant/daAlgorithms/Blue.py b/src/daComposant/daAlgorithms/Blue.py
index 8bc884257787ae9dd191387a88e694fdb310609a..4ab6a00c4e8f0f5204097ff056653ef0e2ad5916 100644 (file)
--- a/src/daComposant/daAlgorithms/Blue.py
+++ b/src/daComposant/daAlgorithms/Blue.py
@@ -1,6 +1,6 @@
 #-*-coding:iso-8859-1-*-
 #
-#  Copyright (C) 2008-2011  EDF R&D
+#  Copyright (C) 2008-2012 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
@@ -74,6 +74,15 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         d  = Y - HXb
         logging.debug("%s Innovation d = %s"%(self._name, d))
         #
+        # Paramètres de pilotage
+        # ----------------------
+        # Potentiels : "CalculateAPosterioriCovariance"
+        if Parameters.has_key("CalculateAPosterioriCovariance"):
+            CalculateAPosterioriCovariance = bool(Parameters["CalculateAPosterioriCovariance"])
+        else:
+            CalculateAPosterioriCovariance = False
+        logging.debug("%s Calcul de la covariance a posteriori = %s"%(self._name, CalculateAPosterioriCovariance))
+        #
         # Calcul de la matrice de gain dans l'espace le plus petit et de l'analyse
         # ------------------------------------------------------------------------
         if Y.size <= Xb.size:
@@ -100,6 +109,12 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         self.StoredVariables["CostFunctionJo"].store( Jo )
         self.StoredVariables["CostFunctionJ" ].store( J )
         #
+        # Calcul de la covariance d'analyse
+        # ---------------------------------
+        if CalculateAPosterioriCovariance:
+            A = ( 1.0 -  K * Hm ) * B
+            self.StoredVariables["APosterioriCovariance"].store( A )
+        #
         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("MB")))
         logging.debug("%s Terminé"%self._name)
         #

diff --git a/src/daComposant/daAlgorithms/EnsembleBlue.py b/src/daComposant/daAlgorithms/EnsembleBlue.py
index d4c13fba2f5c4abdd2e729137c555ec098f89bc0..678609afa7f78107b7e1fc91bbe98e09d3b55be2 100644 (file)
--- a/src/daComposant/daAlgorithms/EnsembleBlue.py
+++ b/src/daComposant/daAlgorithms/EnsembleBlue.py
@@ -1,6 +1,6 @@
 #-*-coding:iso-8859-1-*-
 #
-#  Copyright (C) 2008-2011  EDF R&D
+#  Copyright (C) 2008-2012 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
@@ -60,9 +60,9 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         # -----------------------------------------------------------------
         Hm = H["Direct"].asMatrix()
         Ht = H["Adjoint"].asMatrix()
-        
+        #
         K  = B * Ht * (Hm * B * Ht + R).I
-        
+        #
         # Calcul du BLUE pour chaque membre de l'ensemble
         # -----------------------------------------------
         for iens in range(nb_ens):

diff --git a/src/daComposant/daAlgorithms/KalmanFilter.py b/src/daComposant/daAlgorithms/KalmanFilter.py
index e7b92840a41a34834a55369c974b4eae696a8dc0..683ac6677ef475967a942c84bd5292a8d3883962 100644 (file)
--- a/src/daComposant/daAlgorithms/KalmanFilter.py
+++ b/src/daComposant/daAlgorithms/KalmanFilter.py
@@ -1,6 +1,6 @@
 #-*-coding:iso-8859-1-*-
 #
-#  Copyright (C) 2008-2011  EDF R&D
+#  Copyright (C) 2008-2012 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
@@ -52,6 +52,16 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         Mm = M["Direct"].asMatrix()
         Mt = M["Adjoint"].asMatrix()
         #
+        # Paramètres de pilotage
+        # ----------------------
+        if Parameters.has_key("CalculateAPosterioriCovariance"):
+            CalculateAPosterioriCovariance = bool(Parameters["CalculateAPosterioriCovariance"])
+        else:
+            CalculateAPosterioriCovariance = False
+        logging.debug("%s Calcul de la covariance a posteriori = %s"%(self._name, CalculateAPosterioriCovariance))
+        #
+        # Nombre de pas du Kalman identique au nombre de pas d'observations
+        # -----------------------------------------------------------------
         duration = Y.stepnumber()
         #
         # Initialisation
@@ -59,7 +69,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         Xn = Xb
         Pn = B
         self.StoredVariables["Analysis"].store( Xn.A1 )
-        self.StoredVariables["CovarianceAPosteriori"].store( Pn )
+        if CalculateAPosterioriCovariance:
+            self.StoredVariables["APosterioriCovariance"].store( Pn )
         #
         for step in range(duration-1):
             logging.debug("%s Etape de Kalman %i (i.e. %i->%i) sur un total de %i"%(self._name, step+1, step,step+1, duration-1))
@@ -77,8 +88,9 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             Pn = Pn_predicted - K * Hm * Pn_predicted
             #
             self.StoredVariables["Analysis"].store( Xn.A1 )
-            self.StoredVariables["CovarianceAPosteriori"].store( Pn )
             self.StoredVariables["Innovation"].store( d.A1 )
+            if CalculateAPosterioriCovariance:
+                self.StoredVariables["APosterioriCovariance"].store( Pn )
         #
         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
         logging.debug("%s Terminé"%self._name)

diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py
index b73660146d6ab287fab8460b0890c4cc90fc216c..4f5a216d547ace5e38eda16bd861cf3ba90174e3 100644 (file)
--- a/src/daComposant/daCore/BasicObjects.py
+++ b/src/daComposant/daCore/BasicObjects.py
@@ -1,6 +1,6 @@
 #-*-coding:iso-8859-1-*-
 #
-#  Copyright (C) 2008-2011  EDF R&D
+#  Copyright (C) 2008-2012 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
@@ -137,7 +137,7 @@ class Algorithm:
             - OMA : Observation moins Analysis : Y - Xa
             - OMB : Observation moins Background : Y - Xb
             - AMB : Analysis moins Background : Xa - Xb
-            - CovarianceAPosteriori : matrice A
+            - APosterioriCovariance : matrice A
         On peut rajouter des variables à stocker dans l'initialisation de
         l'algorithme élémentaire qui va hériter de cette classe
         """
@@ -158,7 +158,7 @@ class Algorithm:
         self.StoredVariables["OMA"]                      = Persistence.OneVector(name = "OMA")
         self.StoredVariables["OMB"]                      = Persistence.OneVector(name = "OMB")
         self.StoredVariables["BMA"]                      = Persistence.OneVector(name = "BMA")
-        self.StoredVariables["CovarianceAPosteriori"]    = Persistence.OneMatrix(name = "CovarianceAPosteriori")
+        self.StoredVariables["APosterioriCovariance"]    = Persistence.OneMatrix(name = "APosterioriCovariance")
 
     def get(self, key=None):
         """