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")
194 def get(self, key=None):
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.
203 return self.StoredVariables[key]
205 return self.StoredVariables
207 def has_key(self, key=None):
209 Vérifie si l'une des variables stockées est identifiée par la clé.
211 return self.StoredVariables.has_key(key)
215 Renvoie la liste des clés de variables stockées.
217 return self.StoredVariables.keys()
219 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
221 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
222 sa forme mathématique la plus naturelle possible.
224 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
226 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
228 Permet de définir dans l'algorithme des paramètres requis et leurs
229 caractéristiques par défaut.
232 raise ValueError("A name is mandatory to define a required parameter.")
234 self.__required_parameters[name] = {
236 "typecast" : typecast,
242 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
244 def getRequiredParameters(self, noDetails=True):
246 Renvoie la liste des noms de paramètres requis ou directement le
247 dictionnaire des paramètres requis.
250 ks = self.__required_parameters.keys()
254 return self.__required_parameters
256 def setParameterValue(self, name=None, value=None):
258 Renvoie la valeur d'un paramètre requis de manière contrôlée
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"]
266 if value is None and default is None:
268 elif value is None and default is not None:
269 if typecast is None: __val = default
270 else: __val = typecast( default )
272 if typecast is None: __val = value
273 else: __val = typecast( value )
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:
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))
288 def setParameters(self, fromDico={}):
290 Permet de stocker les paramètres reçus dans le dictionnaire interne.
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])
297 self._parameters[k] = self.setParameterValue(k)
298 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
300 # ==============================================================================
303 Classe générale d'interface de type diagnostic
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
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).
316 def __init__(self, name = "", parameters = {}):
317 self.name = str(name)
318 self.parameters = dict( parameters )
320 def _formula(self, *args):
322 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
323 mathématique la plus naturelle possible.
325 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
327 def calculate(self, *args):
329 Active la formule de calcul avec les arguments correctement rangés
331 raise NotImplementedError("Diagnostic activation method has not been implemented!")
333 # ==============================================================================
336 Classe générale d'interface de type covariance
339 name = "GenericCovariance",
341 asEyeByScalar = None,
342 asEyeByVector = None,
345 Permet de définir une covariance :
346 - asCovariance : entrée des données, comme une matrice compatible avec
347 le constructeur de numpy.matrix
348 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
349 multiplicatif d'une matrice de corrélation identité, aucune matrice
350 n'étant donc explicitement à donner
351 - asEyeByVector : entrée des données comme un seul vecteur de variance,
352 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
353 n'étant donc explicitement à donner
355 self.__name = str(name)
358 self.__is_scalar = False
359 self.__is_vector = False
360 self.__is_matrix = False
361 if asEyeByScalar is not None:
362 self.__is_scalar = True
363 self.__B = float(asEyeByScalar)
366 elif asEyeByVector is not None:
367 self.__is_vector = True
368 self.__B = numpy.array( numpy.ravel( asEyeByVector ), float )
369 self.shape = (self.__B.size,self.__B.size)
370 self.size = self.__B.size**2
371 elif asCovariance is not None:
372 self.__is_matrix = True
373 self.__B = numpy.matrix( asCovariance, float )
374 self.shape = self.__B.shape
375 self.size = self.__B.size
378 # 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)
382 def __validate(self):
383 if self.ismatrix() and min(self.shape) != max(self.shape):
384 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))
385 if self.isscalar() and self.__B <= 0:
386 raise ValueError("The %s covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__B_scalar))
387 if self.isvector() and (self.__B <= 0).any():
388 raise ValueError("The %s covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
389 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
391 L = numpy.linalg.cholesky( self.__B )
393 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
396 return self.__is_scalar
399 return self.__is_vector
402 return self.__is_matrix
406 return Covariance(self.__name+"I", asCovariance = self.__B.I )
407 elif self.isvector():
408 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__B )
409 elif self.isscalar():
410 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__B )
416 return Covariance(self.__name+"T", asCovariance = self.__B.T )
417 elif self.isvector():
418 return Covariance(self.__name+"T", asEyeByVector = self.__B )
419 elif self.isscalar():
420 return Covariance(self.__name+"T", asEyeByScalar = self.__B )
424 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__B) )
425 elif self.isvector():
426 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__B ) )
427 elif self.isscalar():
428 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__B ) )
432 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__B).I )
433 elif self.isvector():
434 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__B ) )
435 elif self.isscalar():
436 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__B ) )
438 def diag(self, msize=None):
440 return numpy.diag(self.__B)
441 elif self.isvector():
443 elif self.isscalar():
445 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,))
447 return self.__B * numpy.ones(int(msize))
449 def trace(self, msize=None):
451 return numpy.trace(self.__B)
452 elif self.isvector():
453 return float(numpy.sum(self.__B))
454 elif self.isscalar():
456 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,))
458 return self.__B * int(msize)
461 return repr(self.__B)
466 def __add__(self, other):
468 return self.__B + numpy.asmatrix(other)
469 elif self.isvector() or self.isscalar():
470 _A = numpy.asarray(other)
471 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__B
472 return numpy.asmatrix(_A)
474 def __radd__(self, other):
475 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
477 def __sub__(self, other):
479 return self.__B - numpy.asmatrix(other)
480 elif self.isvector() or self.isscalar():
481 _A = numpy.asarray(other)
482 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__B - _A.reshape(_A.size)[::_A.shape[1]+1]
483 return numpy.asmatrix(_A)
485 def __rsub__(self, other):
486 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
491 def __mul__(self, other):
492 if self.ismatrix() and isinstance(other,numpy.matrix):
493 return self.__B * other
494 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
495 or isinstance(other,list) \
496 or isinstance(other,tuple)):
497 if numpy.ravel(other).size == self.shape[1]: # Vecteur
498 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
499 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
500 return self.__B * numpy.asmatrix(other)
502 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
503 elif self.isvector() and (isinstance(other,numpy.matrix) \
504 or isinstance(other,numpy.ndarray) \
505 or isinstance(other,list) \
506 or isinstance(other,tuple)):
507 if numpy.ravel(other).size == self.shape[1]: # Vecteur
508 return numpy.asmatrix(self.__B * numpy.ravel(other)).T
509 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
510 return numpy.asmatrix((self.__B * (numpy.asarray(other).transpose())).transpose())
512 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
513 elif self.isscalar() and isinstance(other,numpy.matrix):
514 return self.__B * other
515 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
516 or isinstance(other,list) \
517 or isinstance(other,tuple)):
518 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
519 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
521 return self.__B * numpy.asmatrix(other)
523 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
525 def __rmul__(self, other):
526 if self.ismatrix() and isinstance(other,numpy.matrix):
527 return other * self.__B
528 elif self.isvector() and isinstance(other,numpy.matrix):
529 if numpy.ravel(other).size == self.shape[0]: # Vecteur
530 return numpy.asmatrix(numpy.ravel(other) * self.__B)
531 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
532 return numpy.asmatrix(numpy.array(other) * self.__B)
534 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
535 elif self.isscalar() and isinstance(other,numpy.matrix):
536 return other * self.__B
538 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
543 # ==============================================================================
544 if __name__ == "__main__":
545 print '\n AUTODIAGNOSTIC \n'