Salome HOME
Improving KF and EKF for parameter estimation
[modules/adao.git] / src / daComposant / daCore / BasicObjects.py
1 #-*-coding:iso-8859-1-*-
2 #
3 # Copyright (C) 2008-2013 EDF R&D
4 #
5 # This library is free software; you can redistribute it and/or
6 # modify it under the terms of the GNU Lesser General Public
7 # License as published by the Free Software Foundation; either
8 # version 2.1 of the License.
9 #
10 # This library is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13 # Lesser General Public License for more details.
14 #
15 # You should have received a copy of the GNU Lesser General Public
16 # License along with this library; if not, write to the Free Software
17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
18 #
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
22
23 __doc__ = """
24     Définit les outils généraux élémentaires.
25     
26     Ce module est destiné à etre appelée par AssimilationStudy pour constituer
27     les objets élémentaires de l'algorithme.
28 """
29 __author__ = "Jean-Philippe ARGAUD"
30
31 import logging
32 import numpy
33 import Persistence
34
35 # ==============================================================================
36 class Operator:
37     """
38     Classe générale d'interface de type opérateur
39     """
40     def __init__(self, fromMethod=None, fromMatrix=None):
41         """
42         On construit un objet de ce type en fournissant à l'aide de l'un des
43         deux mots-clé, soit une fonction python, soit matrice.
44         Arguments :
45         - fromMethod : argument de type fonction Python
46         - fromMatrix : argument adapté au constructeur numpy.matrix
47         """
48         if   fromMethod is not None:
49             self.__Method = fromMethod
50             self.__Matrix = None
51             self.__Type   = "Method"
52         elif fromMatrix is not None:
53             self.__Method = None
54             self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
55             self.__Type   = "Matrix"
56         else:
57             self.__Method = None
58             self.__Matrix = None
59             self.__Type   = None
60
61     def isType(self):
62         return self.__Type
63
64     def appliedTo(self, xValue):
65         """
66         Permet de restituer le résultat de l'application de l'opérateur à un
67         argument xValue. Cette méthode se contente d'appliquer, son argument
68         devant a priori être du bon type.
69         Arguments :
70         - xValue : argument adapté pour appliquer l'opérateur
71         """
72         if self.__Matrix is not None:
73             return self.__Matrix * xValue
74         else:
75             return self.__Method( xValue )
76
77     def appliedControledFormTo(self, (xValue, uValue) ):
78         """
79         Permet de restituer le résultat de l'application de l'opérateur à une
80         paire (xValue, uValue). Cette méthode se contente d'appliquer, son
81         argument devant a priori être du bon type. Si la uValue est None,
82         on suppose que l'opérateur ne s'applique qu'à xValue.
83         Arguments :
84         - xValue : argument X adapté pour appliquer l'opérateur
85         - uValue : argument U adapté pour appliquer l'opérateur
86         """
87         if self.__Matrix is not None:
88             return self.__Matrix * xValue
89         elif uValue is not None:
90             return self.__Method( (xValue, uValue) )
91         else:
92             return self.__Method( xValue )
93
94     def appliedInXTo(self, (xNominal, xValue) ):
95         """
96         Permet de restituer le résultat de l'application de l'opérateur à un
97         argument xValue, sachant que l'opérateur est valable en xNominal.
98         Cette méthode se contente d'appliquer, son argument devant a priori
99         être du bon type. Si l'opérateur est linéaire car c'est une matrice,
100         alors il est valable en tout point nominal et il n'est pas nécessaire
101         d'utiliser xNominal.
102         Arguments : une liste contenant
103         - xNominal : argument permettant de donner le point où l'opérateur
104           est construit pour etre ensuite appliqué
105         - xValue : argument adapté pour appliquer l'opérateur
106         """
107         if self.__Matrix is not None:
108             return self.__Matrix * xValue
109         else:
110             return self.__Method( (xNominal, xValue) )
111
112     def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
113         """
114         Permet de renvoyer l'opérateur sous la forme d'une matrice
115         """
116         if self.__Matrix is not None:
117             return self.__Matrix
118         elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
119             return self.__Method( (ValueForMethodForm, None) )
120         else:
121             raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
122
123     def shape(self):
124         """
125         Renvoie la taille sous forme numpy si l'opérateur est disponible sous
126         la forme d'une matrice
127         """
128         if self.__Matrix is not None:
129             return self.__Matrix.shape
130         else:
131             raise ValueError("Matrix form of the operator is not available, nor the shape")
132
133 # ==============================================================================
134 class Algorithm:
135     """
136     Classe générale d'interface de type algorithme
137     
138     Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
139     d'assimilation, en fournissant un container (dictionnaire) de variables
140     persistantes initialisées, et des méthodes d'accès à ces variables stockées.
141     
142     Une classe élémentaire d'algorithme doit implémenter la méthode "run".
143     """
144     def __init__(self, name):
145         """
146         L'initialisation présente permet de fabriquer des variables de stockage
147         disponibles de manière générique dans les algorithmes élémentaires. Ces
148         variables de stockage sont ensuite conservées dans un dictionnaire
149         interne à l'objet, mais auquel on accède par la méthode "get".
150         
151         Les variables prévues sont :
152             - CostFunctionJ  : fonction-cout globale, somme des deux parties suivantes
153             - CostFunctionJb : partie ébauche ou background de la fonction-cout
154             - CostFunctionJo : partie observations de la fonction-cout
155             - GradientOfCostFunctionJ  : gradient de la fonction-cout globale
156             - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
157             - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
158             - CurrentState : état courant lors d'itérations
159             - Analysis : l'analyse
160             - Innovation : l'innovation : d = Y - H Xb
161             - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
162             - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
163             - MahalanobisConsistency : indicateur de consistance des covariances
164             - OMA : Observation moins Analysis : Y - Xa
165             - OMB : Observation moins Background : Y - Xb
166             - AMB : Analysis moins Background : Xa - Xb
167             - APosterioriCovariance : matrice A
168         On peut rajouter des variables à stocker dans l'initialisation de
169         l'algorithme élémentaire qui va hériter de cette classe
170         """
171         logging.debug("%s Initialisation"%str(name))
172         self._name = str( name )
173         self._parameters = {}
174         self.__required_parameters = {}
175         self.StoredVariables = {}
176         #
177         self.StoredVariables["CostFunctionJ"]            = Persistence.OneScalar(name = "CostFunctionJ")
178         self.StoredVariables["CostFunctionJb"]           = Persistence.OneScalar(name = "CostFunctionJb")
179         self.StoredVariables["CostFunctionJo"]           = Persistence.OneScalar(name = "CostFunctionJo")
180         self.StoredVariables["GradientOfCostFunctionJ"]  = Persistence.OneVector(name = "GradientOfCostFunctionJ")
181         self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
182         self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
183         self.StoredVariables["CurrentState"]             = Persistence.OneVector(name = "CurrentState")
184         self.StoredVariables["Analysis"]                 = Persistence.OneVector(name = "Analysis")
185         self.StoredVariables["Innovation"]               = Persistence.OneVector(name = "Innovation")
186         self.StoredVariables["SigmaObs2"]                = Persistence.OneScalar(name = "SigmaObs2")
187         self.StoredVariables["SigmaBck2"]                = Persistence.OneScalar(name = "SigmaBck2")
188         self.StoredVariables["MahalanobisConsistency"]   = Persistence.OneScalar(name = "MahalanobisConsistency")
189         self.StoredVariables["OMA"]                      = Persistence.OneVector(name = "OMA")
190         self.StoredVariables["OMB"]                      = Persistence.OneVector(name = "OMB")
191         self.StoredVariables["BMA"]                      = Persistence.OneVector(name = "BMA")
192         self.StoredVariables["APosterioriCovariance"]    = Persistence.OneMatrix(name = "APosterioriCovariance")
193
194     def get(self, key=None):
195         """
196         Renvoie l'une des variables stockées identifiée par la clé, ou le
197         dictionnaire de l'ensemble des variables disponibles en l'absence de
198         clé. Ce sont directement les variables sous forme objet qui sont
199         renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
200         des classes de persistance.
201         """
202         if key is not None:
203             return self.StoredVariables[key]
204         else:
205             return self.StoredVariables
206
207     def has_key(self, key=None):
208         """
209         Vérifie si l'une des variables stockées est identifiée par la clé.
210         """
211         return self.StoredVariables.has_key(key)
212
213     def keys(self):
214         """
215         Renvoie la liste des clés de variables stockées.
216         """
217         return self.StoredVariables.keys()
218
219     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
220         """
221         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
222         sa forme mathématique la plus naturelle possible.
223         """
224         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
225
226     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
227         """
228         Permet de définir dans l'algorithme des paramètres requis et leurs
229         caractéristiques par défaut.
230         """
231         if name is None:
232             raise ValueError("A name is mandatory to define a required parameter.")
233         #
234         self.__required_parameters[name] = {
235             "default"  : default,
236             "typecast" : typecast,
237             "minval"   : minval,
238             "maxval"   : maxval,
239             "listval"  : listval,
240             "message"  : message,
241             }
242         logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
243
244     def getRequiredParameters(self, noDetails=True):
245         """
246         Renvoie la liste des noms de paramètres requis ou directement le
247         dictionnaire des paramètres requis.
248         """
249         if noDetails:
250             ks = self.__required_parameters.keys()
251             ks.sort()
252             return ks
253         else:
254             return self.__required_parameters
255
256     def setParameterValue(self, name=None, value=None):
257         """
258         Renvoie la valeur d'un paramètre requis de manière contrôlée
259         """
260         default  = self.__required_parameters[name]["default"]
261         typecast = self.__required_parameters[name]["typecast"]
262         minval   = self.__required_parameters[name]["minval"]
263         maxval   = self.__required_parameters[name]["maxval"]
264         listval  = self.__required_parameters[name]["listval"]
265         #
266         if value is None and default is None:
267             __val = None
268         elif value is None and default is not None:
269             if typecast is None: __val = default
270             else:                __val = typecast( default )
271         else:
272             if typecast is None: __val = value
273             else:                __val = typecast( value )
274         #
275         if minval is not None and __val < minval:
276             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
277         if maxval is not None and __val > maxval:
278             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
279         if listval is not None:
280             if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
281                 for v in __val:
282                     if v not in listval:
283                         raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
284             elif __val not in listval:
285                 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
286         return __val
287
288     def setParameters(self, fromDico={}):
289         """
290         Permet de stocker les paramètres reçus dans le dictionnaire interne.
291         """
292         self._parameters.update( fromDico )
293         for k in self.__required_parameters.keys():
294             if k in fromDico.keys():
295                 self._parameters[k] = self.setParameterValue(k,fromDico[k])
296             else:
297                 self._parameters[k] = self.setParameterValue(k)
298             logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
299
300 # ==============================================================================
301 class Diagnostic:
302     """
303     Classe générale d'interface de type diagnostic
304         
305     Ce template s'utilise de la manière suivante : il sert de classe "patron" en
306     même temps que l'une des classes de persistance, comme "OneScalar" par
307     exemple.
308     
309     Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
310     méthode "_formula" pour écrire explicitement et proprement la formule pour
311     l'écriture mathématique du calcul du diagnostic (méthode interne non
312     publique), et "calculate" pour activer la précédente tout en ayant vérifié
313     et préparé les données, et pour stocker les résultats à chaque pas (méthode
314     externe d'activation).
315     """
316     def __init__(self, name = "", parameters = {}):
317         self.name       = str(name)
318         self.parameters = dict( parameters )
319
320     def _formula(self, *args):
321         """
322         Doit implémenter l'opération élémentaire de diagnostic sous sa forme
323         mathématique la plus naturelle possible.
324         """
325         raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
326
327     def calculate(self, *args):
328         """
329         Active la formule de calcul avec les arguments correctement rangés
330         """
331         raise NotImplementedError("Diagnostic activation method has not been implemented!")
332
333 # ==============================================================================
334 if __name__ == "__main__":
335     print '\n AUTODIAGNOSTIC \n'