1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2014 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
43 def __init__(self, fromMethod=None, fromMatrix=None):
45 On construit un objet de ce type en fournissant à l'aide de l'un des
46 deux mots-clé, soit une fonction python, soit matrice.
48 - fromMethod : argument de type fonction Python
49 - fromMatrix : argument adapté au constructeur numpy.matrix
51 self.__NbCallsAsMatrix, self.__NbCallsAsMethod = 0, 0
52 if fromMethod is not None:
53 self.__Method = fromMethod
55 self.__Type = "Method"
56 elif fromMatrix is not None:
58 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
59 self.__Type = "Matrix"
68 def appliedTo(self, xValue):
70 Permet de restituer le résultat de l'application de l'opérateur à un
71 argument xValue. Cette méthode se contente d'appliquer, son argument
72 devant a priori être du bon type.
74 - xValue : argument adapté pour appliquer l'opérateur
76 if self.__Matrix is not None:
77 self.__addOneMatrixCall()
78 return self.__Matrix * xValue
80 self.__addOneMethodCall()
81 return self.__Method( xValue )
83 def appliedControledFormTo(self, (xValue, uValue) ):
85 Permet de restituer le résultat de l'application de l'opérateur à une
86 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
87 argument devant a priori être du bon type. Si la uValue est None,
88 on suppose que l'opérateur ne s'applique qu'à xValue.
90 - xValue : argument X adapté pour appliquer l'opérateur
91 - uValue : argument U adapté pour appliquer l'opérateur
93 if self.__Matrix is not None:
94 self.__addOneMatrixCall()
95 return self.__Matrix * xValue
96 elif uValue is not None:
97 self.__addOneMethodCall()
98 return self.__Method( (xValue, uValue) )
100 self.__addOneMethodCall()
101 return self.__Method( xValue )
103 def appliedInXTo(self, (xNominal, xValue) ):
105 Permet de restituer le résultat de l'application de l'opérateur à un
106 argument xValue, sachant que l'opérateur est valable en xNominal.
107 Cette méthode se contente d'appliquer, son argument devant a priori
108 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
109 alors il est valable en tout point nominal et il n'est pas nécessaire
111 Arguments : une liste contenant
112 - xNominal : argument permettant de donner le point où l'opérateur
113 est construit pour etre ensuite appliqué
114 - xValue : argument adapté pour appliquer l'opérateur
116 if self.__Matrix is not None:
117 self.__addOneMatrixCall()
118 return self.__Matrix * xValue
120 self.__addOneMethodCall()
121 return self.__Method( (xNominal, xValue) )
123 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
125 Permet de renvoyer l'opérateur sous la forme d'une matrice
127 if self.__Matrix is not None:
128 self.__addOneMatrixCall()
130 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
131 self.__addOneMethodCall()
132 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
134 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
138 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
139 la forme d'une matrice
141 if self.__Matrix is not None:
142 return self.__Matrix.shape
144 raise ValueError("Matrix form of the operator is not available, nor the shape")
146 def nbcalls(self, which=None):
148 Renvoie les nombres d'évaluations de l'opérateur
151 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
152 self.__NbCallsAsMatrix,
153 self.__NbCallsAsMethod,
154 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
155 Operator.NbCallsAsMatrix,
156 Operator.NbCallsAsMethod,
158 if which is None: return __nbcalls
159 else: return __nbcalls[which]
161 def __addOneMatrixCall(self):
162 self.__NbCallsAsMatrix += 1 # Decompte local
163 Operator.NbCallsAsMatrix += 1 # Decompte global
165 def __addOneMethodCall(self):
166 self.__NbCallsAsMethod += 1 # Decompte local
167 Operator.NbCallsAsMethod += 1 # Decompte global
169 # ==============================================================================
172 Classe générale d'interface de type algorithme
174 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
175 d'assimilation, en fournissant un container (dictionnaire) de variables
176 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
178 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
180 def __init__(self, name):
182 L'initialisation présente permet de fabriquer des variables de stockage
183 disponibles de manière générique dans les algorithmes élémentaires. Ces
184 variables de stockage sont ensuite conservées dans un dictionnaire
185 interne à l'objet, mais auquel on accède par la méthode "get".
187 Les variables prévues sont :
188 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
189 - CostFunctionJb : partie ébauche ou background de la fonction-cout
190 - CostFunctionJo : partie observations de la fonction-cout
191 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
192 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
193 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
194 - CurrentState : état courant lors d'itérations
195 - Analysis : l'analyse
196 - Innovation : l'innovation : d = Y - H Xb
197 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
198 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
199 - MahalanobisConsistency : indicateur de consistance des covariances
200 - OMA : Observation moins Analysis : Y - Xa
201 - OMB : Observation moins Background : Y - Xb
202 - AMB : Analysis moins Background : Xa - Xb
203 - APosterioriCovariance : matrice A
204 On peut rajouter des variables à stocker dans l'initialisation de
205 l'algorithme élémentaire qui va hériter de cette classe
207 logging.debug("%s Initialisation"%str(name))
208 self._name = str( name )
209 self._parameters = {}
210 self.__required_parameters = {}
211 self.StoredVariables = {}
213 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
214 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
215 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
216 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
217 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
218 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
219 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
220 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
221 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
222 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
223 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
224 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
225 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
226 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
227 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
228 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
229 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
231 def get(self, key=None):
233 Renvoie l'une des variables stockées identifiée par la clé, ou le
234 dictionnaire de l'ensemble des variables disponibles en l'absence de
235 clé. Ce sont directement les variables sous forme objet qui sont
236 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
237 des classes de persistance.
240 return self.StoredVariables[key]
242 return self.StoredVariables
244 def has_key(self, key=None):
246 Vérifie si l'une des variables stockées est identifiée par la clé.
248 return self.StoredVariables.has_key(key)
252 Renvoie la liste des clés de variables stockées.
254 return self.StoredVariables.keys()
256 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
258 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
259 sa forme mathématique la plus naturelle possible.
261 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
263 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
265 Permet de définir dans l'algorithme des paramètres requis et leurs
266 caractéristiques par défaut.
269 raise ValueError("A name is mandatory to define a required parameter.")
271 self.__required_parameters[name] = {
273 "typecast" : typecast,
279 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
281 def getRequiredParameters(self, noDetails=True):
283 Renvoie la liste des noms de paramètres requis ou directement le
284 dictionnaire des paramètres requis.
287 ks = self.__required_parameters.keys()
291 return self.__required_parameters
293 def setParameterValue(self, name=None, value=None):
295 Renvoie la valeur d'un paramètre requis de manière contrôlée
297 default = self.__required_parameters[name]["default"]
298 typecast = self.__required_parameters[name]["typecast"]
299 minval = self.__required_parameters[name]["minval"]
300 maxval = self.__required_parameters[name]["maxval"]
301 listval = self.__required_parameters[name]["listval"]
303 if value is None and default is None:
305 elif value is None and default is not None:
306 if typecast is None: __val = default
307 else: __val = typecast( default )
309 if typecast is None: __val = value
310 else: __val = typecast( value )
312 if minval is not None and __val < minval:
313 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
314 if maxval is not None and __val > maxval:
315 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
316 if listval is not None:
317 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
320 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
321 elif __val not in listval:
322 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
325 def setParameters(self, fromDico={}):
327 Permet de stocker les paramètres reçus dans le dictionnaire interne.
329 self._parameters.update( fromDico )
330 for k in self.__required_parameters.keys():
331 if k in fromDico.keys():
332 self._parameters[k] = self.setParameterValue(k,fromDico[k])
334 self._parameters[k] = self.setParameterValue(k)
335 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
337 # ==============================================================================
340 Classe générale d'interface de type diagnostic
342 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
343 même temps que l'une des classes de persistance, comme "OneScalar" par
346 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
347 méthode "_formula" pour écrire explicitement et proprement la formule pour
348 l'écriture mathématique du calcul du diagnostic (méthode interne non
349 publique), et "calculate" pour activer la précédente tout en ayant vérifié
350 et préparé les données, et pour stocker les résultats à chaque pas (méthode
351 externe d'activation).
353 def __init__(self, name = "", parameters = {}):
354 self.name = str(name)
355 self.parameters = dict( parameters )
357 def _formula(self, *args):
359 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
360 mathématique la plus naturelle possible.
362 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
364 def calculate(self, *args):
366 Active la formule de calcul avec les arguments correctement rangés
368 raise NotImplementedError("Diagnostic activation method has not been implemented!")
370 # ==============================================================================
373 Classe générale d'interface de type covariance
376 name = "GenericCovariance",
378 asEyeByScalar = None,
379 asEyeByVector = None,
382 Permet de définir une covariance :
383 - asCovariance : entrée des données, comme une matrice compatible avec
384 le constructeur de numpy.matrix
385 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
386 multiplicatif d'une matrice de corrélation identité, aucune matrice
387 n'étant donc explicitement à donner
388 - asEyeByVector : entrée des données comme un seul vecteur de variance,
389 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
390 n'étant donc explicitement à donner
392 self.__name = str(name)
395 self.__is_scalar = False
396 self.__is_vector = False
397 self.__is_matrix = False
398 if asEyeByScalar is not None:
399 self.__is_scalar = True
400 self.__B = numpy.abs( float(asEyeByScalar) )
403 elif asEyeByVector is not None:
404 self.__is_vector = True
405 self.__B = numpy.abs( numpy.array( numpy.ravel( asEyeByVector ), float ) )
406 self.shape = (self.__B.size,self.__B.size)
407 self.size = self.__B.size**2
408 elif asCovariance is not None:
409 self.__is_matrix = True
410 self.__B = numpy.matrix( asCovariance, float )
411 self.shape = self.__B.shape
412 self.size = self.__B.size
415 # 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)
419 def __validate(self):
420 if self.ismatrix() and min(self.shape) != max(self.shape):
421 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))
422 if self.isscalar() and self.__B <= 0:
423 raise ValueError("The %s covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__B_scalar))
424 if self.isvector() and (self.__B <= 0).any():
425 raise ValueError("The %s covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
426 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
428 L = numpy.linalg.cholesky( self.__B )
430 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
433 return self.__is_scalar
436 return self.__is_vector
439 return self.__is_matrix
443 return Covariance(self.__name+"I", asCovariance = self.__B.I )
444 elif self.isvector():
445 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__B )
446 elif self.isscalar():
447 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__B )
453 return Covariance(self.__name+"T", asCovariance = self.__B.T )
454 elif self.isvector():
455 return Covariance(self.__name+"T", asEyeByVector = self.__B )
456 elif self.isscalar():
457 return Covariance(self.__name+"T", asEyeByScalar = self.__B )
461 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__B) )
462 elif self.isvector():
463 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__B ) )
464 elif self.isscalar():
465 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__B ) )
469 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__B).I )
470 elif self.isvector():
471 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__B ) )
472 elif self.isscalar():
473 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__B ) )
475 def diag(self, msize=None):
477 return numpy.diag(self.__B)
478 elif self.isvector():
480 elif self.isscalar():
482 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,))
484 return self.__B * numpy.ones(int(msize))
486 def asfullmatrix(self, msize=None):
489 elif self.isvector():
490 return numpy.matrix( numpy.diag(self.__B), float )
491 elif self.isscalar():
493 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,))
495 return numpy.matrix( self.__B * numpy.eye(int(msize)), float )
497 def trace(self, msize=None):
499 return numpy.trace(self.__B)
500 elif self.isvector():
501 return float(numpy.sum(self.__B))
502 elif self.isscalar():
504 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,))
506 return self.__B * int(msize)
509 return repr(self.__B)
514 def __add__(self, other):
516 return self.__B + numpy.asmatrix(other)
517 elif self.isvector() or self.isscalar():
518 _A = numpy.asarray(other)
519 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__B
520 return numpy.asmatrix(_A)
522 def __radd__(self, other):
523 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
525 def __sub__(self, other):
527 return self.__B - numpy.asmatrix(other)
528 elif self.isvector() or self.isscalar():
529 _A = numpy.asarray(other)
530 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__B - _A.reshape(_A.size)[::_A.shape[1]+1]
531 return numpy.asmatrix(_A)
533 def __rsub__(self, other):
534 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
539 def __mul__(self, other):
540 if self.ismatrix() and isinstance(other,numpy.matrix):
541 return self.__B * other
542 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
543 or isinstance(other,list) \
544 or isinstance(other,tuple)):
545 if numpy.ravel(other).size == self.shape[1]: # Vecteur
546 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
547 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
548 return self.__B * numpy.asmatrix(other)
550 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
551 elif self.isvector() and (isinstance(other,numpy.matrix) \
552 or isinstance(other,numpy.ndarray) \
553 or isinstance(other,list) \
554 or isinstance(other,tuple)):
555 if numpy.ravel(other).size == self.shape[1]: # Vecteur
556 return numpy.asmatrix(self.__B * numpy.ravel(other)).T
557 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
558 return numpy.asmatrix((self.__B * (numpy.asarray(other).transpose())).transpose())
560 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
561 elif self.isscalar() and isinstance(other,numpy.matrix):
562 return self.__B * other
563 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
564 or isinstance(other,list) \
565 or isinstance(other,tuple)):
566 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
567 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
569 return self.__B * numpy.asmatrix(other)
571 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
573 def __rmul__(self, other):
574 if self.ismatrix() and isinstance(other,numpy.matrix):
575 return other * self.__B
576 elif self.isvector() and isinstance(other,numpy.matrix):
577 if numpy.ravel(other).size == self.shape[0]: # Vecteur
578 return numpy.asmatrix(numpy.ravel(other) * self.__B)
579 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
580 return numpy.asmatrix(numpy.array(other) * self.__B)
582 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
583 elif self.isscalar() and isinstance(other,numpy.matrix):
584 return other * self.__B
586 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
591 # ==============================================================================
592 if __name__ == "__main__":
593 print '\n AUTODIAGNOSTIC \n'