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 self.__NbCallsAsMatrix, self.__NbCallsAsMethod = 0, 0
49 if fromMethod is not None:
50 self.__Method = fromMethod
52 self.__Type = "Method"
53 elif fromMatrix is not None:
55 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
56 self.__Type = "Matrix"
65 def appliedTo(self, xValue):
67 Permet de restituer le résultat de l'application de l'opérateur à un
68 argument xValue. Cette méthode se contente d'appliquer, son argument
69 devant a priori être du bon type.
71 - xValue : argument adapté pour appliquer l'opérateur
73 if self.__Matrix is not None:
74 self.__NbCallsAsMatrix += 1
75 return self.__Matrix * xValue
77 self.__NbCallsAsMethod += 1
78 return self.__Method( xValue )
80 def appliedControledFormTo(self, (xValue, uValue) ):
82 Permet de restituer le résultat de l'application de l'opérateur à une
83 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
84 argument devant a priori être du bon type. Si la uValue est None,
85 on suppose que l'opérateur ne s'applique qu'à xValue.
87 - xValue : argument X adapté pour appliquer l'opérateur
88 - uValue : argument U adapté pour appliquer l'opérateur
90 if self.__Matrix is not None:
91 self.__NbCallsAsMatrix += 1
92 return self.__Matrix * xValue
93 elif uValue is not None:
94 self.__NbCallsAsMethod += 1
95 return self.__Method( (xValue, uValue) )
97 self.__NbCallsAsMethod += 1
98 return self.__Method( xValue )
100 def appliedInXTo(self, (xNominal, xValue) ):
102 Permet de restituer le résultat de l'application de l'opérateur à un
103 argument xValue, sachant que l'opérateur est valable en xNominal.
104 Cette méthode se contente d'appliquer, son argument devant a priori
105 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
106 alors il est valable en tout point nominal et il n'est pas nécessaire
108 Arguments : une liste contenant
109 - xNominal : argument permettant de donner le point où l'opérateur
110 est construit pour etre ensuite appliqué
111 - xValue : argument adapté pour appliquer l'opérateur
113 if self.__Matrix is not None:
114 self.__NbCallsAsMatrix += 1
115 return self.__Matrix * xValue
117 self.__NbCallsAsMethod += 1
118 return self.__Method( (xNominal, xValue) )
120 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
122 Permet de renvoyer l'opérateur sous la forme d'une matrice
124 if self.__Matrix is not None:
125 self.__NbCallsAsMatrix += 1
127 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
128 self.__NbCallsAsMethod += 1
129 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
131 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
135 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
136 la forme d'une matrice
138 if self.__Matrix is not None:
139 return self.__Matrix.shape
141 raise ValueError("Matrix form of the operator is not available, nor the shape")
145 Renvoie le nombre d'évaluations de l'opérateurs (total, matrice, méthode)
147 return (self.__NbCallsAsMatrix+self.__NbCallsAsMethod,self.__NbCallsAsMatrix,self.__NbCallsAsMethod)
149 # ==============================================================================
152 Classe générale d'interface de type algorithme
154 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
155 d'assimilation, en fournissant un container (dictionnaire) de variables
156 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
158 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
160 def __init__(self, name):
162 L'initialisation présente permet de fabriquer des variables de stockage
163 disponibles de manière générique dans les algorithmes élémentaires. Ces
164 variables de stockage sont ensuite conservées dans un dictionnaire
165 interne à l'objet, mais auquel on accède par la méthode "get".
167 Les variables prévues sont :
168 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
169 - CostFunctionJb : partie ébauche ou background de la fonction-cout
170 - CostFunctionJo : partie observations de la fonction-cout
171 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
172 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
173 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
174 - CurrentState : état courant lors d'itérations
175 - Analysis : l'analyse
176 - Innovation : l'innovation : d = Y - H Xb
177 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
178 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
179 - MahalanobisConsistency : indicateur de consistance des covariances
180 - OMA : Observation moins Analysis : Y - Xa
181 - OMB : Observation moins Background : Y - Xb
182 - AMB : Analysis moins Background : Xa - Xb
183 - APosterioriCovariance : matrice A
184 On peut rajouter des variables à stocker dans l'initialisation de
185 l'algorithme élémentaire qui va hériter de cette classe
187 logging.debug("%s Initialisation"%str(name))
188 self._name = str( name )
189 self._parameters = {}
190 self.__required_parameters = {}
191 self.StoredVariables = {}
193 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
194 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
195 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
196 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
197 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
198 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
199 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
200 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
201 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
202 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
203 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
204 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
205 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
206 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
207 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
208 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
209 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
211 def get(self, key=None):
213 Renvoie l'une des variables stockées identifiée par la clé, ou le
214 dictionnaire de l'ensemble des variables disponibles en l'absence de
215 clé. Ce sont directement les variables sous forme objet qui sont
216 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
217 des classes de persistance.
220 return self.StoredVariables[key]
222 return self.StoredVariables
224 def has_key(self, key=None):
226 Vérifie si l'une des variables stockées est identifiée par la clé.
228 return self.StoredVariables.has_key(key)
232 Renvoie la liste des clés de variables stockées.
234 return self.StoredVariables.keys()
236 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
238 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
239 sa forme mathématique la plus naturelle possible.
241 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
243 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
245 Permet de définir dans l'algorithme des paramètres requis et leurs
246 caractéristiques par défaut.
249 raise ValueError("A name is mandatory to define a required parameter.")
251 self.__required_parameters[name] = {
253 "typecast" : typecast,
259 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
261 def getRequiredParameters(self, noDetails=True):
263 Renvoie la liste des noms de paramètres requis ou directement le
264 dictionnaire des paramètres requis.
267 ks = self.__required_parameters.keys()
271 return self.__required_parameters
273 def setParameterValue(self, name=None, value=None):
275 Renvoie la valeur d'un paramètre requis de manière contrôlée
277 default = self.__required_parameters[name]["default"]
278 typecast = self.__required_parameters[name]["typecast"]
279 minval = self.__required_parameters[name]["minval"]
280 maxval = self.__required_parameters[name]["maxval"]
281 listval = self.__required_parameters[name]["listval"]
283 if value is None and default is None:
285 elif value is None and default is not None:
286 if typecast is None: __val = default
287 else: __val = typecast( default )
289 if typecast is None: __val = value
290 else: __val = typecast( value )
292 if minval is not None and __val < minval:
293 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
294 if maxval is not None and __val > maxval:
295 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
296 if listval is not None:
297 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
300 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
301 elif __val not in listval:
302 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
305 def setParameters(self, fromDico={}):
307 Permet de stocker les paramètres reçus dans le dictionnaire interne.
309 self._parameters.update( fromDico )
310 for k in self.__required_parameters.keys():
311 if k in fromDico.keys():
312 self._parameters[k] = self.setParameterValue(k,fromDico[k])
314 self._parameters[k] = self.setParameterValue(k)
315 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
317 # ==============================================================================
320 Classe générale d'interface de type diagnostic
322 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
323 même temps que l'une des classes de persistance, comme "OneScalar" par
326 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
327 méthode "_formula" pour écrire explicitement et proprement la formule pour
328 l'écriture mathématique du calcul du diagnostic (méthode interne non
329 publique), et "calculate" pour activer la précédente tout en ayant vérifié
330 et préparé les données, et pour stocker les résultats à chaque pas (méthode
331 externe d'activation).
333 def __init__(self, name = "", parameters = {}):
334 self.name = str(name)
335 self.parameters = dict( parameters )
337 def _formula(self, *args):
339 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
340 mathématique la plus naturelle possible.
342 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
344 def calculate(self, *args):
346 Active la formule de calcul avec les arguments correctement rangés
348 raise NotImplementedError("Diagnostic activation method has not been implemented!")
350 # ==============================================================================
353 Classe générale d'interface de type covariance
356 name = "GenericCovariance",
358 asEyeByScalar = None,
359 asEyeByVector = None,
362 Permet de définir une covariance :
363 - asCovariance : entrée des données, comme une matrice compatible avec
364 le constructeur de numpy.matrix
365 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
366 multiplicatif d'une matrice de corrélation identité, aucune matrice
367 n'étant donc explicitement à donner
368 - asEyeByVector : entrée des données comme un seul vecteur de variance,
369 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
370 n'étant donc explicitement à donner
372 self.__name = str(name)
375 self.__is_scalar = False
376 self.__is_vector = False
377 self.__is_matrix = False
378 if asEyeByScalar is not None:
379 self.__is_scalar = True
380 self.__B = numpy.abs( float(asEyeByScalar) )
383 elif asEyeByVector is not None:
384 self.__is_vector = True
385 self.__B = numpy.abs( numpy.array( numpy.ravel( asEyeByVector ), float ) )
386 self.shape = (self.__B.size,self.__B.size)
387 self.size = self.__B.size**2
388 elif asCovariance is not None:
389 self.__is_matrix = True
390 self.__B = numpy.matrix( asCovariance, float )
391 self.shape = self.__B.shape
392 self.size = self.__B.size
395 # 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)
399 def __validate(self):
400 if self.ismatrix() and min(self.shape) != max(self.shape):
401 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))
402 if self.isscalar() and self.__B <= 0:
403 raise ValueError("The %s covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__B_scalar))
404 if self.isvector() and (self.__B <= 0).any():
405 raise ValueError("The %s covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
406 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
408 L = numpy.linalg.cholesky( self.__B )
410 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
413 return self.__is_scalar
416 return self.__is_vector
419 return self.__is_matrix
423 return Covariance(self.__name+"I", asCovariance = self.__B.I )
424 elif self.isvector():
425 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__B )
426 elif self.isscalar():
427 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__B )
433 return Covariance(self.__name+"T", asCovariance = self.__B.T )
434 elif self.isvector():
435 return Covariance(self.__name+"T", asEyeByVector = self.__B )
436 elif self.isscalar():
437 return Covariance(self.__name+"T", asEyeByScalar = self.__B )
441 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__B) )
442 elif self.isvector():
443 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__B ) )
444 elif self.isscalar():
445 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__B ) )
449 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__B).I )
450 elif self.isvector():
451 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__B ) )
452 elif self.isscalar():
453 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__B ) )
455 def diag(self, msize=None):
457 return numpy.diag(self.__B)
458 elif self.isvector():
460 elif self.isscalar():
462 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,))
464 return self.__B * numpy.ones(int(msize))
466 def asfullmatrix(self, msize=None):
469 elif self.isvector():
470 return numpy.matrix( numpy.diag(self.__B), float )
471 elif self.isscalar():
473 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,))
475 return numpy.matrix( self.__B * numpy.eye(int(msize)), float )
477 def trace(self, msize=None):
479 return numpy.trace(self.__B)
480 elif self.isvector():
481 return float(numpy.sum(self.__B))
482 elif self.isscalar():
484 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,))
486 return self.__B * int(msize)
489 return repr(self.__B)
494 def __add__(self, other):
496 return self.__B + numpy.asmatrix(other)
497 elif self.isvector() or self.isscalar():
498 _A = numpy.asarray(other)
499 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__B
500 return numpy.asmatrix(_A)
502 def __radd__(self, other):
503 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
505 def __sub__(self, other):
507 return self.__B - numpy.asmatrix(other)
508 elif self.isvector() or self.isscalar():
509 _A = numpy.asarray(other)
510 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__B - _A.reshape(_A.size)[::_A.shape[1]+1]
511 return numpy.asmatrix(_A)
513 def __rsub__(self, other):
514 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
519 def __mul__(self, other):
520 if self.ismatrix() and isinstance(other,numpy.matrix):
521 return self.__B * other
522 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
523 or isinstance(other,list) \
524 or isinstance(other,tuple)):
525 if numpy.ravel(other).size == self.shape[1]: # Vecteur
526 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
527 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
528 return self.__B * numpy.asmatrix(other)
530 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
531 elif self.isvector() and (isinstance(other,numpy.matrix) \
532 or isinstance(other,numpy.ndarray) \
533 or isinstance(other,list) \
534 or isinstance(other,tuple)):
535 if numpy.ravel(other).size == self.shape[1]: # Vecteur
536 return numpy.asmatrix(self.__B * numpy.ravel(other)).T
537 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
538 return numpy.asmatrix((self.__B * (numpy.asarray(other).transpose())).transpose())
540 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
541 elif self.isscalar() and isinstance(other,numpy.matrix):
542 return self.__B * other
543 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
544 or isinstance(other,list) \
545 or isinstance(other,tuple)):
546 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
547 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
549 return self.__B * numpy.asmatrix(other)
551 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
553 def __rmul__(self, other):
554 if self.ismatrix() and isinstance(other,numpy.matrix):
555 return other * self.__B
556 elif self.isvector() and isinstance(other,numpy.matrix):
557 if numpy.ravel(other).size == self.shape[0]: # Vecteur
558 return numpy.asmatrix(numpy.ravel(other) * self.__B)
559 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
560 return numpy.asmatrix(numpy.array(other) * self.__B)
562 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
563 elif self.isscalar() and isinstance(other,numpy.matrix):
564 return other * self.__B
566 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
571 # ==============================================================================
572 if __name__ == "__main__":
573 print '\n AUTODIAGNOSTIC \n'