1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2015 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"
37 # ==============================================================================
40 Classe générale de gestion d'un cache de calculs
43 toleranceInRedundancy = 1.e-18,
44 lenghtOfRedundancy = -1,
47 Les caractéristiques de tolérance peuvent être modifées à la création.
49 self.__tolerBP = float(toleranceInRedundancy)
50 self.__lenghtOR = int(lenghtOfRedundancy)
51 self.__initlnOR = self.__lenghtOR
55 self.__listOPCV = [] # Operator Previous Calculated Points, Results, Point Norms
56 # logging.debug("CM Tolerance de determination des doublons : %.2e"%self.__tolerBP)
58 def wasCalculatedIn(self, xValue ): #, info="" ):
61 for i in xrange(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
62 if xValue.size != self.__listOPCV[i][0].size:
63 # logging.debug("CM Différence de la taille %s de X et de celle %s du point %i déjà calculé"%(xValue.shape,i,self.__listOPCP[i].shape))
65 if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
67 __HxV = self.__listOPCV[i][1]
68 # logging.debug("CM Cas%s déja calculé, portant le numéro %i"%(info,i))
72 def storeValueInX(self, xValue, HxValue ):
73 if self.__lenghtOR < 0:
74 self.__lenghtOR = 2 * xValue.size + 2
75 self.__initlnOR = self.__lenghtOR
76 while len(self.__listOPCV) > self.__lenghtOR:
77 # logging.debug("CM Réduction de la liste des cas à %i éléments par suppression du premier"%self.__lenghtOR)
78 self.__listOPCV.pop(0)
79 self.__listOPCV.append( (
80 copy.copy(numpy.ravel(xValue)),
82 numpy.linalg.norm(xValue),
86 self.__initlnOR = self.__lenghtOR
90 self.__lenghtOR = self.__initlnOR
92 # ==============================================================================
95 Classe générale d'interface de type opérateur
102 def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
104 On construit un objet de ce type en fournissant à l'aide de l'un des
105 deux mots-clé, soit une fonction python, soit une matrice.
107 - fromMethod : argument de type fonction Python
108 - fromMatrix : argument adapté au constructeur numpy.matrix
109 - avoidingRedundancy : évite ou pas les calculs redondants
111 self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
112 self.__AvoidRC = bool( avoidingRedundancy )
113 if fromMethod is not None:
114 self.__Method = fromMethod
116 self.__Type = "Method"
117 elif fromMatrix is not None:
119 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
120 self.__Type = "Matrix"
126 def disableAvoidingRedundancy(self):
127 Operator.CM.disable()
129 def enableAvoidingRedundancy(self):
133 Operator.CM.disable()
138 def appliedTo(self, xValue):
140 Permet de restituer le résultat de l'application de l'opérateur à un
141 argument xValue. Cette méthode se contente d'appliquer, son argument
142 devant a priori être du bon type.
144 - xValue : argument adapté pour appliquer l'opérateur
147 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
149 __alreadyCalculated = False
151 if __alreadyCalculated:
152 self.__addOneCacheCall()
155 if self.__Matrix is not None:
156 self.__addOneMatrixCall()
157 HxValue = self.__Matrix * xValue
159 self.__addOneMethodCall()
160 HxValue = self.__Method( xValue )
162 Operator.CM.storeValueInX(xValue,HxValue)
166 def appliedControledFormTo(self, (xValue, uValue) ):
168 Permet de restituer le résultat de l'application de l'opérateur à une
169 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
170 argument devant a priori être du bon type. Si la uValue est None,
171 on suppose que l'opérateur ne s'applique qu'à xValue.
173 - xValue : argument X adapté pour appliquer l'opérateur
174 - uValue : argument U adapté pour appliquer l'opérateur
176 if self.__Matrix is not None:
177 self.__addOneMatrixCall()
178 return self.__Matrix * xValue
179 elif uValue is not None:
180 self.__addOneMethodCall()
181 return self.__Method( (xValue, uValue) )
183 self.__addOneMethodCall()
184 return self.__Method( xValue )
186 def appliedInXTo(self, (xNominal, xValue) ):
188 Permet de restituer le résultat de l'application de l'opérateur à un
189 argument xValue, sachant que l'opérateur est valable en xNominal.
190 Cette méthode se contente d'appliquer, son argument devant a priori
191 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
192 alors il est valable en tout point nominal et il n'est pas nécessaire
194 Arguments : une liste contenant
195 - xNominal : argument permettant de donner le point où l'opérateur
196 est construit pour etre ensuite appliqué
197 - xValue : argument adapté pour appliquer l'opérateur
199 if self.__Matrix is not None:
200 self.__addOneMatrixCall()
201 return self.__Matrix * xValue
203 self.__addOneMethodCall()
204 return self.__Method( (xNominal, xValue) )
206 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
208 Permet de renvoyer l'opérateur sous la forme d'une matrice
210 if self.__Matrix is not None:
211 self.__addOneMatrixCall()
213 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
214 self.__addOneMethodCall()
215 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
217 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
221 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
222 la forme d'une matrice
224 if self.__Matrix is not None:
225 return self.__Matrix.shape
227 raise ValueError("Matrix form of the operator is not available, nor the shape")
229 def nbcalls(self, which=None):
231 Renvoie les nombres d'évaluations de l'opérateur
234 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
235 self.__NbCallsAsMatrix,
236 self.__NbCallsAsMethod,
237 self.__NbCallsOfCached,
238 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
239 Operator.NbCallsAsMatrix,
240 Operator.NbCallsAsMethod,
241 Operator.NbCallsOfCached,
243 if which is None: return __nbcalls
244 else: return __nbcalls[which]
246 def __addOneMatrixCall(self):
247 self.__NbCallsAsMatrix += 1 # Decompte local
248 Operator.NbCallsAsMatrix += 1 # Decompte global
250 def __addOneMethodCall(self):
251 self.__NbCallsAsMethod += 1 # Decompte local
252 Operator.NbCallsAsMethod += 1 # Decompte global
254 def __addOneCacheCall(self):
255 self.__NbCallsOfCached += 1 # Decompte local
256 Operator.NbCallsOfCached += 1 # Decompte global
258 # ==============================================================================
261 Classe générale d'interface de type algorithme
263 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
264 d'assimilation, en fournissant un container (dictionnaire) de variables
265 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
267 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
269 def __init__(self, name):
271 L'initialisation présente permet de fabriquer des variables de stockage
272 disponibles de manière générique dans les algorithmes élémentaires. Ces
273 variables de stockage sont ensuite conservées dans un dictionnaire
274 interne à l'objet, mais auquel on accède par la méthode "get".
276 Les variables prévues sont :
277 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
278 - CostFunctionJb : partie ébauche ou background de la fonction-cout
279 - CostFunctionJo : partie observations de la fonction-cout
280 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
281 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
282 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
283 - CurrentState : état courant lors d'itérations
284 - Analysis : l'analyse Xa
285 - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
286 - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
287 - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
288 - Innovation : l'innovation : d = Y - H(X)
289 - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
290 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
291 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
292 - MahalanobisConsistency : indicateur de consistance des covariances
293 - OMA : Observation moins Analysis : Y - Xa
294 - OMB : Observation moins Background : Y - Xb
295 - AMB : Analysis moins Background : Xa - Xb
296 - APosterioriCovariance : matrice A
297 - APosterioriVariances : variances de la matrice A
298 - APosterioriStandardDeviations : écart-types de la matrice A
299 - APosterioriCorrelations : correlations de la matrice A
300 On peut rajouter des variables à stocker dans l'initialisation de
301 l'algorithme élémentaire qui va hériter de cette classe
303 logging.debug("%s Initialisation"%str(name))
304 self._m = PlatformInfo.SystemUsage()
306 self._name = str( name )
307 self._parameters = {"StoreSupplementaryCalculations":[]}
308 self.__required_parameters = {}
309 self.StoredVariables = {}
311 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
312 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
313 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
314 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
315 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
316 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
317 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
318 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
319 self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum")
320 self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
321 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
322 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
323 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
324 self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
325 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
326 self.StoredVariables["InnovationAtCurrentState"] = Persistence.OneVector(name = "InnovationAtCurrentState")
327 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
328 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
329 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
330 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
331 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
332 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
333 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
334 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
335 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
336 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
337 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
340 logging.debug("%s Lancement"%self._name)
341 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
344 def _post_run(self,_oH=None):
345 if ("StoreSupplementaryCalculations" in self._parameters) and \
346 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
347 for _A in self.StoredVariables["APosterioriCovariance"]:
348 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
349 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
350 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
351 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
352 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
353 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
354 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
355 self.StoredVariables["APosterioriCorrelations"].store( _C )
357 logging.debug("%s Nombre d'évaluation(s) de l'opérateur d'observation direct/tangent/adjoint.: %i/%i/%i"%(self._name, _oH["Direct"].nbcalls(0),_oH["Tangent"].nbcalls(0),_oH["Adjoint"].nbcalls(0)))
358 logging.debug("%s Nombre d'appels au cache d'opérateur d'observation direct/tangent/adjoint..: %i/%i/%i"%(self._name, _oH["Direct"].nbcalls(3),_oH["Tangent"].nbcalls(3),_oH["Adjoint"].nbcalls(3)))
359 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
360 logging.debug("%s Terminé"%self._name)
363 def get(self, key=None):
365 Renvoie l'une des variables stockées identifiée par la clé, ou le
366 dictionnaire de l'ensemble des variables disponibles en l'absence de
367 clé. Ce sont directement les variables sous forme objet qui sont
368 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
369 des classes de persistance.
372 return self.StoredVariables[key]
374 return self.StoredVariables
376 def __contains__(self, key=None):
378 Vérifie si l'une des variables stockées est identifiée par la clé.
380 return (key in self.StoredVariables)
384 Renvoie la liste des clés de variables stockées.
386 return self.StoredVariables.keys()
388 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
390 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
391 sa forme mathématique la plus naturelle possible.
393 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
395 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
397 Permet de définir dans l'algorithme des paramètres requis et leurs
398 caractéristiques par défaut.
401 raise ValueError("A name is mandatory to define a required parameter.")
403 self.__required_parameters[name] = {
405 "typecast" : typecast,
411 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
413 def getRequiredParameters(self, noDetails=True):
415 Renvoie la liste des noms de paramètres requis ou directement le
416 dictionnaire des paramètres requis.
419 ks = self.__required_parameters.keys()
423 return self.__required_parameters
425 def setParameterValue(self, name=None, value=None):
427 Renvoie la valeur d'un paramètre requis de manière contrôlée
429 default = self.__required_parameters[name]["default"]
430 typecast = self.__required_parameters[name]["typecast"]
431 minval = self.__required_parameters[name]["minval"]
432 maxval = self.__required_parameters[name]["maxval"]
433 listval = self.__required_parameters[name]["listval"]
435 if value is None and default is None:
437 elif value is None and default is not None:
438 if typecast is None: __val = default
439 else: __val = typecast( default )
441 if typecast is None: __val = value
442 else: __val = typecast( value )
444 if minval is not None and (numpy.array(__val, float) < minval).any():
445 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
446 if maxval is not None and (numpy.array(__val, float) > maxval).any():
447 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
448 if listval is not None:
449 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
452 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
453 elif __val not in listval:
454 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
457 def setParameters(self, fromDico={}):
459 Permet de stocker les paramètres reçus dans le dictionnaire interne.
461 self._parameters.update( fromDico )
462 for k in self.__required_parameters.keys():
463 if k in fromDico.keys():
464 self._parameters[k] = self.setParameterValue(k,fromDico[k])
466 self._parameters[k] = self.setParameterValue(k)
467 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
469 # ==============================================================================
472 Classe générale d'interface de type diagnostic
474 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
475 même temps que l'une des classes de persistance, comme "OneScalar" par
478 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
479 méthode "_formula" pour écrire explicitement et proprement la formule pour
480 l'écriture mathématique du calcul du diagnostic (méthode interne non
481 publique), et "calculate" pour activer la précédente tout en ayant vérifié
482 et préparé les données, et pour stocker les résultats à chaque pas (méthode
483 externe d'activation).
485 def __init__(self, name = "", parameters = {}):
486 self.name = str(name)
487 self.parameters = dict( parameters )
489 def _formula(self, *args):
491 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
492 mathématique la plus naturelle possible.
494 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
496 def calculate(self, *args):
498 Active la formule de calcul avec les arguments correctement rangés
500 raise NotImplementedError("Diagnostic activation method has not been implemented!")
502 # ==============================================================================
505 Classe générale d'interface de type covariance
508 name = "GenericCovariance",
510 asEyeByScalar = None,
511 asEyeByVector = None,
516 Permet de définir une covariance :
517 - asCovariance : entrée des données, comme une matrice compatible avec
518 le constructeur de numpy.matrix
519 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
520 multiplicatif d'une matrice de corrélation identité, aucune matrice
521 n'étant donc explicitement à donner
522 - asEyeByVector : entrée des données comme un seul vecteur de variance,
523 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
524 n'étant donc explicitement à donner
525 - asCovObject : entrée des données comme un objet python, qui a les
526 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
527 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
528 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
529 - toBeChecked : booléen indiquant si le caractère SDP de la matrice
530 pleine doit être vérifié
532 self.__name = str(name)
533 self.__check = bool(toBeChecked)
536 self.__is_scalar = False
537 self.__is_vector = False
538 self.__is_matrix = False
539 self.__is_object = False
540 if asEyeByScalar is not None:
541 self.__is_scalar = True
542 self.__C = numpy.abs( float(asEyeByScalar) )
545 elif asEyeByVector is not None:
546 self.__is_vector = True
547 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
548 self.shape = (self.__C.size,self.__C.size)
549 self.size = self.__C.size**2
550 elif asCovariance is not None:
551 self.__is_matrix = True
552 self.__C = numpy.matrix( asCovariance, float )
553 self.shape = self.__C.shape
554 self.size = self.__C.size
555 elif asCovObject is not None:
556 self.__is_object = True
557 self.__C = asCovObject
558 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
559 if not hasattr(self.__C,at):
560 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
561 if hasattr(self.__C,"shape"):
562 self.shape = self.__C.shape
565 if hasattr(self.__C,"size"):
566 self.size = self.__C.size
571 # 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)
575 def __validate(self):
576 if self.ismatrix() and min(self.shape) != max(self.shape):
577 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))
578 if self.isobject() and min(self.shape) != max(self.shape):
579 raise ValueError("The matrix given for \"%s\" is not a square one, its shape is %s. Please check your object input."%(self.__name,self.shape))
580 if self.isscalar() and self.__C <= 0:
581 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
582 if self.isvector() and (self.__C <= 0).any():
583 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
584 if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
586 L = numpy.linalg.cholesky( self.__C )
588 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
589 if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
591 L = self.__C.cholesky()
593 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
596 return self.__is_scalar
599 return self.__is_vector
602 return self.__is_matrix
605 return self.__is_object
609 return Covariance(self.__name+"I", asCovariance = self.__C.I )
610 elif self.isvector():
611 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
612 elif self.isscalar():
613 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
614 elif self.isobject():
615 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
617 return None # Indispensable
621 return Covariance(self.__name+"T", asCovariance = self.__C.T )
622 elif self.isvector():
623 return Covariance(self.__name+"T", asEyeByVector = self.__C )
624 elif self.isscalar():
625 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
626 elif self.isobject():
627 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
631 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
632 elif self.isvector():
633 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
634 elif self.isscalar():
635 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
636 elif self.isobject() and hasattr(self.__C,"cholesky"):
637 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
641 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
642 elif self.isvector():
643 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
644 elif self.isscalar():
645 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
646 elif self.isobject() and hasattr(self.__C,"choleskyI"):
647 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
649 def diag(self, msize=None):
651 return numpy.diag(self.__C)
652 elif self.isvector():
654 elif self.isscalar():
656 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,))
658 return self.__C * numpy.ones(int(msize))
659 elif self.isobject():
660 return self.__C.diag()
662 def asfullmatrix(self, msize=None):
665 elif self.isvector():
666 return numpy.matrix( numpy.diag(self.__C), float )
667 elif self.isscalar():
669 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,))
671 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
672 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
673 return self.__C.asfullmatrix()
675 def trace(self, msize=None):
677 return numpy.trace(self.__C)
678 elif self.isvector():
679 return float(numpy.sum(self.__C))
680 elif self.isscalar():
682 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,))
684 return self.__C * int(msize)
685 elif self.isobject():
686 return self.__C.trace()
689 return repr(self.__C)
694 def __add__(self, other):
695 if self.ismatrix() or self.isobject():
696 return self.__C + numpy.asmatrix(other)
697 elif self.isvector() or self.isscalar():
698 _A = numpy.asarray(other)
699 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
700 return numpy.asmatrix(_A)
702 def __radd__(self, other):
703 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
705 def __sub__(self, other):
706 if self.ismatrix() or self.isobject():
707 return self.__C - numpy.asmatrix(other)
708 elif self.isvector() or self.isscalar():
709 _A = numpy.asarray(other)
710 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
711 return numpy.asmatrix(_A)
713 def __rsub__(self, other):
714 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
719 def __mul__(self, other):
720 if self.ismatrix() and isinstance(other,numpy.matrix):
721 return self.__C * other
722 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
723 or isinstance(other,list) \
724 or isinstance(other,tuple)):
725 if numpy.ravel(other).size == self.shape[1]: # Vecteur
726 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
727 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
728 return self.__C * numpy.asmatrix(other)
730 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
731 elif self.isvector() and (isinstance(other,numpy.matrix) \
732 or isinstance(other,numpy.ndarray) \
733 or isinstance(other,list) \
734 or isinstance(other,tuple)):
735 if numpy.ravel(other).size == self.shape[1]: # Vecteur
736 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
737 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
738 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
740 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
741 elif self.isscalar() and isinstance(other,numpy.matrix):
742 return self.__C * other
743 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
744 or isinstance(other,list) \
745 or isinstance(other,tuple)):
746 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
747 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
749 return self.__C * numpy.asmatrix(other)
750 elif self.isobject():
751 return self.__C.__mul__(other)
753 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
755 def __rmul__(self, other):
756 if self.ismatrix() and isinstance(other,numpy.matrix):
757 return other * self.__C
758 elif self.isvector() and isinstance(other,numpy.matrix):
759 if numpy.ravel(other).size == self.shape[0]: # Vecteur
760 return numpy.asmatrix(numpy.ravel(other) * self.__C)
761 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
762 return numpy.asmatrix(numpy.array(other) * self.__C)
764 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
765 elif self.isscalar() and isinstance(other,numpy.matrix):
766 return other * self.__C
767 elif self.isobject():
768 return self.__C.__rmul__(other)
770 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
775 # ==============================================================================
776 class TemplateStorage(object):
778 Classe générale de stockage de type dictionnaire étendu
781 def __init__( self, language = "fr_FR" ):
782 self.__preferedLanguage = "fr"
785 def store( self, name = None, content = None, fr_FR = "", en_EN = "" ):
786 if name is None or content is None:
787 raise ValueError("To be consistent, the storage of a template must provide a name and a content.")
788 self.__values[str(name)] = {
789 'content': str(content),
790 'fr_FR' : str(fr_FR),
791 'en_EN' : str(en_EN),
795 __keys = self.__values.keys()
800 return len(self.__values)
802 def __getitem__(self, name=None ):
803 return self.__values[name]['content']
805 def getdoc(self, name = None, lang = "fr_FR"):
806 if lang not in self.__values[name]: lang = "fr_FR"
807 return self.__values[name][lang]
809 # ==============================================================================
812 _Hm = None, # Pour simuler Hm(x) : HO["Direct"].appliedTo
813 _HmX = None, # Simulation déjà faite de Hm(x)
814 _arg = None, # Arguments supplementaires pour Hm, sous la forme d'un tuple
819 _SIV = False, # A résorber pour la 8.0
820 _SSC = [], # self._parameters["StoreSupplementaryCalculations"]
821 _nPS = 0, # nbPreviousSteps
822 _QM = "DA", # QualityMeasure
823 _SSV = {}, # Entrée et/ou sortie : self.StoredVariables
824 _fRt = False, # Restitue ou pas la sortie étendue
825 _sSc = True, # Stocke ou pas les SSC
828 Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
829 et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
830 DFO, QuantileRegression
836 for k in ["CostFunctionJ",
842 "SimulatedObservationAtCurrentOptimum",
843 "SimulatedObservationAtCurrentState",
847 if hasattr(_SSV[k],"store"):
848 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
850 _X = numpy.asmatrix(numpy.ravel( _x )).T
851 if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
852 _SSV["CurrentState"].append( _X )
858 raise ValueError("%s Operator has to be defined."%(self.__name,))
862 _HX = _Hm( _X, *_arg )
863 _HX = numpy.asmatrix(numpy.ravel( _HX )).T
865 if "SimulatedObservationAtCurrentState" in _SSC or \
866 "SimulatedObservationAtCurrentOptimum" in _SSC:
867 _SSV["SimulatedObservationAtCurrentState"].append( _HX )
869 if numpy.any(numpy.isnan(_HX)):
870 Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
872 _Y = numpy.asmatrix(numpy.ravel( _Y )).T
873 if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
874 if _BI is None or _RI is None:
875 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
876 _Xb = numpy.asmatrix(numpy.ravel( _Xb )).T
877 Jb = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
878 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
879 elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
881 raise ValueError("Observation error covariance matrix has to be properly defined!")
883 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
884 elif _QM in ["LeastSquares", "LS", "L2"]:
886 Jo = 0.5 * (_Y - _HX).T * (_Y - _HX)
887 elif _QM in ["AbsoluteValue", "L1"]:
889 Jo = numpy.sum( numpy.abs(_Y - _HX) )
890 elif _QM in ["MaximumError", "ME"]:
892 Jo = numpy.max( numpy.abs(_Y - _HX) )
893 elif _QM in ["QR", "Null"]:
897 raise ValueError("Unknown asked quality measure!")
899 J = float( Jb ) + float( Jo )
902 _SSV["CostFunctionJb"].append( Jb )
903 _SSV["CostFunctionJo"].append( Jo )
904 _SSV["CostFunctionJ" ].append( J )
906 if "IndexOfOptimum" in _SSC or \
907 "CurrentOptimum" in _SSC or \
908 "SimulatedObservationAtCurrentOptimum" in _SSC:
909 IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
910 if "IndexOfOptimum" in _SSC:
911 _SSV["IndexOfOptimum"].append( IndexMin )
912 if "CurrentOptimum" in _SSC:
913 _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
914 if "SimulatedObservationAtCurrentOptimum" in _SSC:
915 _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
920 if _QM in ["QR"]: # Pour le QuantileRegression
925 # ==============================================================================
926 if __name__ == "__main__":
927 print '\n AUTODIAGNOSTIC \n'