1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2013 EDF R&D
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.
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.
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
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 Définit les outils généraux élémentaires.
26 Ce module est destiné à etre appelée par AssimilationStudy pour constituer
27 les objets élémentaires de l'algorithme.
29 __author__ = "Jean-Philippe ARGAUD"
35 # ==============================================================================
38 Classe générale d'interface de type opérateur
40 def __init__(self, fromMethod=None, fromMatrix=None):
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.
45 - fromMethod : argument de type fonction Python
46 - fromMatrix : argument adapté au constructeur numpy.matrix
48 if fromMethod is not None:
49 self.__Method = fromMethod
51 self.__Type = "Method"
52 elif fromMatrix is not None:
54 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
55 self.__Type = "Matrix"
64 def appliedTo(self, xValue):
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.
70 - xValue : argument adapté pour appliquer l'opérateur
72 if self.__Matrix is not None:
73 return self.__Matrix * xValue
75 return self.__Method( xValue )
77 def appliedControledFormTo(self, (xValue, uValue) ):
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.
84 - xValue : argument X adapté pour appliquer l'opérateur
85 - uValue : argument U adapté pour appliquer l'opérateur
87 if self.__Matrix is not None:
88 return self.__Matrix * xValue
89 elif uValue is not None:
90 return self.__Method( (xValue, uValue) )
92 return self.__Method( xValue )
94 def appliedInXTo(self, (xNominal, xValue) ):
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
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
107 if self.__Matrix is not None:
108 return self.__Matrix * xValue
110 return self.__Method( (xNominal, xValue) )
112 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
114 Permet de renvoyer l'opérateur sous la forme d'une matrice
116 if self.__Matrix is not None:
118 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
119 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
121 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
125 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
126 la forme d'une matrice
128 if self.__Matrix is not None:
129 return self.__Matrix.shape
131 raise ValueError("Matrix form of the operator is not available, nor the shape")
133 # ==============================================================================
136 Classe générale d'interface de type algorithme
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.
142 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
144 def __init__(self, name):
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".
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
171 logging.debug("%s Initialisation"%str(name))
172 self._name = str( name )
173 self._parameters = {}
174 self.__required_parameters = {}
175 self.StoredVariables = {}
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 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
195 def get(self, key=None):
197 Renvoie l'une des variables stockées identifiée par la clé, ou le
198 dictionnaire de l'ensemble des variables disponibles en l'absence de
199 clé. Ce sont directement les variables sous forme objet qui sont
200 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
201 des classes de persistance.
204 return self.StoredVariables[key]
206 return self.StoredVariables
208 def has_key(self, key=None):
210 Vérifie si l'une des variables stockées est identifiée par la clé.
212 return self.StoredVariables.has_key(key)
216 Renvoie la liste des clés de variables stockées.
218 return self.StoredVariables.keys()
220 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
222 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
223 sa forme mathématique la plus naturelle possible.
225 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
227 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
229 Permet de définir dans l'algorithme des paramètres requis et leurs
230 caractéristiques par défaut.
233 raise ValueError("A name is mandatory to define a required parameter.")
235 self.__required_parameters[name] = {
237 "typecast" : typecast,
243 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
245 def getRequiredParameters(self, noDetails=True):
247 Renvoie la liste des noms de paramètres requis ou directement le
248 dictionnaire des paramètres requis.
251 ks = self.__required_parameters.keys()
255 return self.__required_parameters
257 def setParameterValue(self, name=None, value=None):
259 Renvoie la valeur d'un paramètre requis de manière contrôlée
261 default = self.__required_parameters[name]["default"]
262 typecast = self.__required_parameters[name]["typecast"]
263 minval = self.__required_parameters[name]["minval"]
264 maxval = self.__required_parameters[name]["maxval"]
265 listval = self.__required_parameters[name]["listval"]
267 if value is None and default is None:
269 elif value is None and default is not None:
270 if typecast is None: __val = default
271 else: __val = typecast( default )
273 if typecast is None: __val = value
274 else: __val = typecast( value )
276 if minval is not None and __val < minval:
277 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
278 if maxval is not None and __val > maxval:
279 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
280 if listval is not None:
281 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
284 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
285 elif __val not in listval:
286 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
289 def setParameters(self, fromDico={}):
291 Permet de stocker les paramètres reçus dans le dictionnaire interne.
293 self._parameters.update( fromDico )
294 for k in self.__required_parameters.keys():
295 if k in fromDico.keys():
296 self._parameters[k] = self.setParameterValue(k,fromDico[k])
298 self._parameters[k] = self.setParameterValue(k)
299 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
301 # ==============================================================================
304 Classe générale d'interface de type diagnostic
306 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
307 même temps que l'une des classes de persistance, comme "OneScalar" par
310 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
311 méthode "_formula" pour écrire explicitement et proprement la formule pour
312 l'écriture mathématique du calcul du diagnostic (méthode interne non
313 publique), et "calculate" pour activer la précédente tout en ayant vérifié
314 et préparé les données, et pour stocker les résultats à chaque pas (méthode
315 externe d'activation).
317 def __init__(self, name = "", parameters = {}):
318 self.name = str(name)
319 self.parameters = dict( parameters )
321 def _formula(self, *args):
323 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
324 mathématique la plus naturelle possible.
326 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
328 def calculate(self, *args):
330 Active la formule de calcul avec les arguments correctement rangés
332 raise NotImplementedError("Diagnostic activation method has not been implemented!")
334 # ==============================================================================
337 Classe générale d'interface de type covariance
340 name = "GenericCovariance",
342 asEyeByScalar = None,
343 asEyeByVector = None,
346 Permet de définir une covariance :
347 - asCovariance : entrée des données, comme une matrice compatible avec
348 le constructeur de numpy.matrix
349 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
350 multiplicatif d'une matrice de corrélation identité, aucune matrice
351 n'étant donc explicitement à donner
352 - asEyeByVector : entrée des données comme un seul vecteur de variance,
353 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
354 n'étant donc explicitement à donner
356 self.__name = str(name)
359 self.__is_scalar = False
360 self.__is_vector = False
361 self.__is_matrix = False
362 if asEyeByScalar is not None:
363 self.__is_scalar = True
364 self.__B = numpy.abs( float(asEyeByScalar) )
367 elif asEyeByVector is not None:
368 self.__is_vector = True
369 self.__B = numpy.abs( numpy.array( numpy.ravel( asEyeByVector ), float ) )
370 self.shape = (self.__B.size,self.__B.size)
371 self.size = self.__B.size**2
372 elif asCovariance is not None:
373 self.__is_matrix = True
374 self.__B = numpy.matrix( asCovariance, float )
375 self.shape = self.__B.shape
376 self.size = self.__B.size
379 # raise ValueError("The %s covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix."%self.__name)
383 def __validate(self):
384 if self.ismatrix() and min(self.shape) != max(self.shape):
385 raise ValueError("The given matrix for %s is not a square one, its shape is %s. Please check your matrix input."%(self.__name,self.shape))
386 if self.isscalar() and self.__B <= 0:
387 raise ValueError("The %s covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__B_scalar))
388 if self.isvector() and (self.__B <= 0).any():
389 raise ValueError("The %s covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
390 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
392 L = numpy.linalg.cholesky( self.__B )
394 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
397 return self.__is_scalar
400 return self.__is_vector
403 return self.__is_matrix
407 return Covariance(self.__name+"I", asCovariance = self.__B.I )
408 elif self.isvector():
409 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__B )
410 elif self.isscalar():
411 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__B )
417 return Covariance(self.__name+"T", asCovariance = self.__B.T )
418 elif self.isvector():
419 return Covariance(self.__name+"T", asEyeByVector = self.__B )
420 elif self.isscalar():
421 return Covariance(self.__name+"T", asEyeByScalar = self.__B )
425 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__B) )
426 elif self.isvector():
427 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__B ) )
428 elif self.isscalar():
429 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__B ) )
433 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__B).I )
434 elif self.isvector():
435 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__B ) )
436 elif self.isscalar():
437 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__B ) )
439 def diag(self, msize=None):
441 return numpy.diag(self.__B)
442 elif self.isvector():
444 elif self.isscalar():
446 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
448 return self.__B * numpy.ones(int(msize))
450 def asfullmatrix(self, msize=None):
453 elif self.isvector():
454 return numpy.matrix( numpy.diag(self.__B), float )
455 elif self.isscalar():
457 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
459 return numpy.matrix( self.__B * numpy.eye(int(msize)), float )
461 def trace(self, msize=None):
463 return numpy.trace(self.__B)
464 elif self.isvector():
465 return float(numpy.sum(self.__B))
466 elif self.isscalar():
468 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
470 return self.__B * int(msize)
473 return repr(self.__B)
478 def __add__(self, other):
480 return self.__B + numpy.asmatrix(other)
481 elif self.isvector() or self.isscalar():
482 _A = numpy.asarray(other)
483 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__B
484 return numpy.asmatrix(_A)
486 def __radd__(self, other):
487 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
489 def __sub__(self, other):
491 return self.__B - numpy.asmatrix(other)
492 elif self.isvector() or self.isscalar():
493 _A = numpy.asarray(other)
494 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__B - _A.reshape(_A.size)[::_A.shape[1]+1]
495 return numpy.asmatrix(_A)
497 def __rsub__(self, other):
498 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
503 def __mul__(self, other):
504 if self.ismatrix() and isinstance(other,numpy.matrix):
505 return self.__B * other
506 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
507 or isinstance(other,list) \
508 or isinstance(other,tuple)):
509 if numpy.ravel(other).size == self.shape[1]: # Vecteur
510 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
511 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
512 return self.__B * numpy.asmatrix(other)
514 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
515 elif self.isvector() and (isinstance(other,numpy.matrix) \
516 or isinstance(other,numpy.ndarray) \
517 or isinstance(other,list) \
518 or isinstance(other,tuple)):
519 if numpy.ravel(other).size == self.shape[1]: # Vecteur
520 return numpy.asmatrix(self.__B * numpy.ravel(other)).T
521 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
522 return numpy.asmatrix((self.__B * (numpy.asarray(other).transpose())).transpose())
524 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
525 elif self.isscalar() and isinstance(other,numpy.matrix):
526 return self.__B * other
527 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
528 or isinstance(other,list) \
529 or isinstance(other,tuple)):
530 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
531 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
533 return self.__B * numpy.asmatrix(other)
535 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
537 def __rmul__(self, other):
538 if self.ismatrix() and isinstance(other,numpy.matrix):
539 return other * self.__B
540 elif self.isvector() and isinstance(other,numpy.matrix):
541 if numpy.ravel(other).size == self.shape[0]: # Vecteur
542 return numpy.asmatrix(numpy.ravel(other) * self.__B)
543 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
544 return numpy.asmatrix(numpy.array(other) * self.__B)
546 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
547 elif self.isscalar() and isinstance(other,numpy.matrix):
548 return other * self.__B
550 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
555 # ==============================================================================
556 if __name__ == "__main__":
557 print '\n AUTODIAGNOSTIC \n'