From 02ca9780759176fdd699254889b1138012d026dc Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Tue, 3 Apr 2012 15:30:32 +0200 Subject: [PATCH] Setting a posteriori covariance calculations optionnal --- src/daComposant/daAlgorithms/3DVAR.py | 22 ++++++++++++++++++-- src/daComposant/daAlgorithms/Blue.py | 17 ++++++++++++++- src/daComposant/daAlgorithms/EnsembleBlue.py | 6 +++--- src/daComposant/daAlgorithms/KalmanFilter.py | 18 +++++++++++++--- src/daComposant/daCore/BasicObjects.py | 6 +++--- 5 files changed, 57 insertions(+), 12 deletions(-) diff --git a/src/daComposant/daAlgorithms/3DVAR.py b/src/daComposant/daAlgorithms/3DVAR.py index cee5c22..6210227 100644 --- 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 8bc8842..4ab6a00 100644 --- 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 d4c13fb..678609a 100644 --- 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 e7b9284..683ac66 100644 --- 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 b736601..4f5a216 100644 --- 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): """ -- 2.39.2