From eafbed91b9225a7576ec1d5ac22aa9b4e37facba Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Sat, 20 Jan 2018 18:46:32 +0100 Subject: [PATCH] Minor corrections and output variables checking --- .../daAlgorithms/EnsembleKalmanFilter.py | 85 +++++++++++++------ src/daComposant/daAlgorithms/KalmanFilter.py | 8 +- src/daComposant/daCore/Aidsm.py | 4 +- 3 files changed, 66 insertions(+), 31 deletions(-) diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py index b1605bb..7d74ed6 100644 --- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py +++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py @@ -25,8 +25,6 @@ from daCore import BasicObjects, PlatformInfo import numpy, math mfp = PlatformInfo.PlatformInfo().MaximumPrecision() -# Using "Ensemble Kalman filtering", L. HOUTEKAMER and HERSCHEL L. MITCHELL, QJRMS (2005), 131, pp. 3269–3289 - # ============================================================================== class ElementaryAlgorithm(BasicObjects.Algorithm): def __init__(self): @@ -38,20 +36,6 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): message = "Nombre de membres dans l'ensemble", minval = -1, ) - self.defineRequiredParameter( - name = "Minimizer", - default = "EnKF", - typecast = str, - message = "Schéma de mise a jour des informations d'ensemble", - listval = ["EnKF", "ETKF", "DEnKF"], - ) - self.defineRequiredParameter( - name = "ConstrainedBy", - default = "EstimateProjection", - typecast = str, - message = "Prise en compte des contraintes", - listval = ["EstimateProjection"], - ) self.defineRequiredParameter( name = "EstimationOf", default = "State", @@ -59,6 +43,11 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): message = "Estimation d'etat ou de parametres", listval = ["State", "Parameters"], ) + self.defineRequiredParameter( + name = "SetSeed", + typecast = numpy.random.seed, + message = "Graine fixée pour le générateur aléatoire", + ) self.defineRequiredParameter( name = "StoreInternalVariables", default = False, @@ -70,16 +59,18 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): default = [], typecast = tuple, message = "Liste de calculs supplémentaires à stocker et/ou effectuer", - listval = ["APosterioriCorrelations", "APosterioriCovariance", "APosterioriStandardDeviations", "APosterioriVariances", "BMA", "CurrentState", "CostFunctionJ", "CostFunctionJb", "CostFunctionJo", "Innovation"] - ) - self.defineRequiredParameter( # Pas de type - name = "Bounds", - message = "Liste des valeurs de bornes", - ) - self.defineRequiredParameter( - name = "SetSeed", - typecast = numpy.random.seed, - message = "Graine fixée pour le générateur aléatoire", + listval = [ + "APosterioriCorrelations", + "APosterioriCovariance", + "APosterioriStandardDeviations", + "APosterioriVariances", + "BMA", + "CostFunctionJ", + "CostFunctionJb", + "CostFunctionJo", + "CurrentState", + "Innovation", + ] ) self.requireInputArguments( mandatory= ("Xb", "Y", "HO", "R", "B"), @@ -115,7 +106,11 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): # # Précalcul des inversions de B et R # ---------------------------------- - if self._parameters["StoreInternalVariables"]: + if self._parameters["StoreInternalVariables"] \ + or "CostFunctionJ" in self._parameters["StoreSupplementaryCalculations"] \ + or "CostFunctionJb" in self._parameters["StoreSupplementaryCalculations"] \ + or "CostFunctionJo" in self._parameters["StoreSupplementaryCalculations"] \ + or "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]: BI = B.getI() RI = R.getI() BIdemi = B.choleskyI() @@ -198,6 +193,28 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): self.StoredVariables["Analysis"].store( Xa ) # del Yo, PfHT, HPfHT + if self._parameters["StoreInternalVariables"] \ + or "CostFunctionJ" in self._parameters["StoreSupplementaryCalculations"] \ + or "CostFunctionJb" in self._parameters["StoreSupplementaryCalculations"] \ + or "CostFunctionJo" in self._parameters["StoreSupplementaryCalculations"] \ + or "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"] \ + or "Innovation" in self._parameters["StoreSupplementaryCalculations"]: + d = Ynpu - numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T + self.StoredVariables["Innovation"].store( d ) + if self._parameters["StoreInternalVariables"] \ + or "CurrentState" in self._parameters["StoreSupplementaryCalculations"]: + self.StoredVariables["CurrentState"].store( Xn ) + if self._parameters["StoreInternalVariables"] \ + or "CostFunctionJ" in self._parameters["StoreSupplementaryCalculations"] \ + or "CostFunctionJb" in self._parameters["StoreSupplementaryCalculations"] \ + or "CostFunctionJo" in self._parameters["StoreSupplementaryCalculations"] \ + or "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]: + Jb = 0.5 * (Xa - Xb).T * BI * (Xa - Xb) + Jo = 0.5 * d.T * RI * d + J = float( Jb ) + float( Jo ) + self.StoredVariables["CostFunctionJb"].store( Jb ) + self.StoredVariables["CostFunctionJo"].store( Jo ) + self.StoredVariables["CostFunctionJ" ].store( J ) if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]: Ht = HO["Tangent"].asMatrix(ValueForMethodForm = Xa) Ht = Ht.reshape(__p,__n) # ADAO & check shape @@ -207,6 +224,20 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): Pf = (1./(__m-1)) * Pf Pn = (1. - K * Ht) * Pf self.StoredVariables["APosterioriCovariance"].store( Pn ) + if J < previousJMinimum: + previousJMinimum = J + Xa = Xn + covarianceXa = Pn + # + # Stockage supplementaire de l'optimum en estimation de parametres + # ---------------------------------------------------------------- + if self._parameters["EstimationOf"] == "Parameters": + self.StoredVariables["Analysis"].store( Xa.A1 ) + if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]: + self.StoredVariables["APosterioriCovariance"].store( covarianceXa ) + # + if "BMA" in self._parameters["StoreSupplementaryCalculations"]: + self.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) ) # self._post_run(HO) return 0 diff --git a/src/daComposant/daAlgorithms/KalmanFilter.py b/src/daComposant/daAlgorithms/KalmanFilter.py index fa413cd..b6bed33 100644 --- a/src/daComposant/daAlgorithms/KalmanFilter.py +++ b/src/daComposant/daAlgorithms/KalmanFilter.py @@ -82,7 +82,10 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): # # Précalcul des inversions de B et R # ---------------------------------- - if self._parameters["StoreInternalVariables"]: + if self._parameters["StoreInternalVariables"] \ + or "CostFunctionJ" in self._parameters["StoreSupplementaryCalculations"] \ + or "CostFunctionJb" in self._parameters["StoreSupplementaryCalculations"] \ + or "CostFunctionJo" in self._parameters["StoreSupplementaryCalculations"]: BI = B.getI() RI = R.getI() # @@ -143,7 +146,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): self.StoredVariables["APosterioriCovariance"].store( Pn ) if "Innovation" in self._parameters["StoreSupplementaryCalculations"]: self.StoredVariables["Innovation"].store( numpy.ravel( d.A1 ) ) - if self._parameters["StoreInternalVariables"] or "CurrentState" in self._parameters["StoreSupplementaryCalculations"]: + if self._parameters["StoreInternalVariables"] \ + or "CurrentState" in self._parameters["StoreSupplementaryCalculations"]: self.StoredVariables["CurrentState"].store( Xn ) if self._parameters["StoreInternalVariables"] \ or "CostFunctionJ" in self._parameters["StoreSupplementaryCalculations"] \ diff --git a/src/daComposant/daCore/Aidsm.py b/src/daComposant/daCore/Aidsm.py index cc7736f..54e9ca3 100644 --- a/src/daComposant/daCore/Aidsm.py +++ b/src/daComposant/daCore/Aidsm.py @@ -689,10 +689,10 @@ class Aidsm(object): except Exception as e: if isinstance(e, SyntaxError): msg = "at %s: %s"%(e.offset, e.text) else: msg = "" - raise ValueError("during execution, the following error occurs:\n"+\ + raise ValueError(("during execution, the following error occurs:\n"+\ "\n%s %s\n\nSee also the potential messages, "+\ "which can show the origin of the above error, "+\ - "in the launching terminal.\n"%(str(e),msg)) + "in the launching terminal.\n")%(str(e),msg)) return 0 def __executePythonScheme(self, FileName=None): -- 2.39.2