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é à être appelée par AssimilationStudy.
28 __author__ = "Jean-Philippe ARGAUD"
36 # ==============================================================================
37 class CacheManager(object):
39 Classe générale de gestion d'un cache de calculs
42 toleranceInRedundancy = 1.e-18,
43 lenghtOfRedundancy = -1,
46 Les caractéristiques de tolérance peuvent être modifées à la création.
48 self.__tolerBP = float(toleranceInRedundancy)
49 self.__lenghtOR = int(lenghtOfRedundancy)
50 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="" ):
59 "Vérifie l'existence d'un calcul correspondant à la valeur"
62 for i in xrange(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
63 if xValue.size != self.__listOPCV[i][0].size:
64 # 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)
66 if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
68 __HxV = self.__listOPCV[i][1]
69 # logging.debug("CM Cas%s déja calculé, portant le numéro %i", info, i)
73 def storeValueInX(self, xValue, HxValue ):
74 "Stocke un calcul correspondant à la valeur"
75 if self.__lenghtOR < 0:
76 self.__lenghtOR = 2 * xValue.size + 2
77 self.__initlnOR = self.__lenghtOR
78 while len(self.__listOPCV) > self.__lenghtOR:
79 # logging.debug("CM Réduction de la liste des cas à %i éléments par suppression du premier", self.__lenghtOR)
80 self.__listOPCV.pop(0)
81 self.__listOPCV.append( (
82 copy.copy(numpy.ravel(xValue)),
84 numpy.linalg.norm(xValue),
89 self.__initlnOR = self.__lenghtOR
94 self.__lenghtOR = self.__initlnOR
96 # ==============================================================================
97 class Operator(object):
99 Classe générale d'interface de type opérateur
106 def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
108 On construit un objet de ce type en fournissant à l'aide de l'un des
109 deux mots-clé, soit une fonction python, soit une matrice.
111 - fromMethod : argument de type fonction Python
112 - fromMatrix : argument adapté au constructeur numpy.matrix
113 - avoidingRedundancy : évite ou pas les calculs redondants
115 self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
116 self.__AvoidRC = bool( avoidingRedundancy )
117 if fromMethod is not None:
118 self.__Method = fromMethod
120 self.__Type = "Method"
121 elif fromMatrix is not None:
123 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
124 self.__Type = "Matrix"
130 def disableAvoidingRedundancy(self):
132 Operator.CM.disable()
134 def enableAvoidingRedundancy(self):
139 Operator.CM.disable()
145 def appliedTo(self, xValue):
147 Permet de restituer le résultat de l'application de l'opérateur à un
148 argument xValue. Cette méthode se contente d'appliquer, son argument
149 devant a priori être du bon type.
151 - xValue : argument adapté pour appliquer l'opérateur
154 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
156 __alreadyCalculated = False
158 if __alreadyCalculated:
159 self.__addOneCacheCall()
162 if self.__Matrix is not None:
163 self.__addOneMatrixCall()
164 HxValue = self.__Matrix * xValue
166 self.__addOneMethodCall()
167 HxValue = self.__Method( xValue )
169 Operator.CM.storeValueInX(xValue,HxValue)
173 def appliedControledFormTo(self, (xValue, uValue) ):
175 Permet de restituer le résultat de l'application de l'opérateur à une
176 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
177 argument devant a priori être du bon type. Si la uValue est None,
178 on suppose que l'opérateur ne s'applique qu'à xValue.
180 - xValue : argument X adapté pour appliquer l'opérateur
181 - uValue : argument U adapté pour appliquer l'opérateur
183 if self.__Matrix is not None:
184 self.__addOneMatrixCall()
185 return self.__Matrix * xValue
186 elif uValue is not None:
187 self.__addOneMethodCall()
188 return self.__Method( (xValue, uValue) )
190 self.__addOneMethodCall()
191 return self.__Method( xValue )
193 def appliedInXTo(self, (xNominal, xValue) ):
195 Permet de restituer le résultat de l'application de l'opérateur à un
196 argument xValue, sachant que l'opérateur est valable en xNominal.
197 Cette méthode se contente d'appliquer, son argument devant a priori
198 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
199 alors il est valable en tout point nominal et il n'est pas nécessaire
201 Arguments : une liste contenant
202 - xNominal : argument permettant de donner le point où l'opérateur
203 est construit pour etre ensuite appliqué
204 - xValue : argument adapté pour appliquer l'opérateur
206 if self.__Matrix is not None:
207 self.__addOneMatrixCall()
208 return self.__Matrix * xValue
210 self.__addOneMethodCall()
211 return self.__Method( (xNominal, xValue) )
213 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
215 Permet de renvoyer l'opérateur sous la forme d'une matrice
217 if self.__Matrix is not None:
218 self.__addOneMatrixCall()
220 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
221 self.__addOneMethodCall()
222 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
224 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
228 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
229 la forme d'une matrice
231 if self.__Matrix is not None:
232 return self.__Matrix.shape
234 raise ValueError("Matrix form of the operator is not available, nor the shape")
236 def nbcalls(self, which=None):
238 Renvoie les nombres d'évaluations de l'opérateur
241 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
242 self.__NbCallsAsMatrix,
243 self.__NbCallsAsMethod,
244 self.__NbCallsOfCached,
245 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
246 Operator.NbCallsAsMatrix,
247 Operator.NbCallsAsMethod,
248 Operator.NbCallsOfCached,
250 if which is None: return __nbcalls
251 else: return __nbcalls[which]
253 def __addOneMatrixCall(self):
254 "Comptabilise un appel"
255 self.__NbCallsAsMatrix += 1 # Decompte local
256 Operator.NbCallsAsMatrix += 1 # Decompte global
258 def __addOneMethodCall(self):
259 "Comptabilise un appel"
260 self.__NbCallsAsMethod += 1 # Decompte local
261 Operator.NbCallsAsMethod += 1 # Decompte global
263 def __addOneCacheCall(self):
264 "Comptabilise un appel"
265 self.__NbCallsOfCached += 1 # Decompte local
266 Operator.NbCallsOfCached += 1 # Decompte global
268 # ==============================================================================
269 class Algorithm(object):
271 Classe générale d'interface de type algorithme
273 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
274 d'assimilation, en fournissant un container (dictionnaire) de variables
275 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
277 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
279 def __init__(self, name):
281 L'initialisation présente permet de fabriquer des variables de stockage
282 disponibles de manière générique dans les algorithmes élémentaires. Ces
283 variables de stockage sont ensuite conservées dans un dictionnaire
284 interne à l'objet, mais auquel on accède par la méthode "get".
286 Les variables prévues sont :
287 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
288 - CostFunctionJb : partie ébauche ou background de la fonction-cout
289 - CostFunctionJo : partie observations de la fonction-cout
290 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
291 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
292 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
293 - CurrentState : état courant lors d'itérations
294 - Analysis : l'analyse Xa
295 - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
296 - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
297 - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
298 - Innovation : l'innovation : d = Y - H(X)
299 - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
300 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
301 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
302 - MahalanobisConsistency : indicateur de consistance des covariances
303 - OMA : Observation moins Analysis : Y - Xa
304 - OMB : Observation moins Background : Y - Xb
305 - AMB : Analysis moins Background : Xa - Xb
306 - APosterioriCovariance : matrice A
307 - APosterioriVariances : variances de la matrice A
308 - APosterioriStandardDeviations : écart-types de la matrice A
309 - APosterioriCorrelations : correlations de la matrice A
310 - Residu : dans le cas des algorithmes de vérification
311 On peut rajouter des variables à stocker dans l'initialisation de
312 l'algorithme élémentaire qui va hériter de cette classe
314 logging.debug("%s Initialisation", str(name))
315 self._m = PlatformInfo.SystemUsage()
317 self._name = str( name )
318 self._parameters = {"StoreSupplementaryCalculations":[]}
319 self.__required_parameters = {}
320 self.StoredVariables = {}
322 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
323 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
324 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
325 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
326 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
327 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
328 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
329 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
330 self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum")
331 self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
332 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
333 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
334 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
335 self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
336 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
337 self.StoredVariables["InnovationAtCurrentState"] = Persistence.OneVector(name = "InnovationAtCurrentState")
338 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
339 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
340 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
341 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
342 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
343 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
344 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
345 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
346 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
347 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
348 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
349 self.StoredVariables["Residu"] = Persistence.OneScalar(name = "Residu")
353 logging.debug("%s Lancement", self._name)
354 logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
357 def _post_run(self,_oH=None):
359 if ("StoreSupplementaryCalculations" in self._parameters) and \
360 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
361 for _A in self.StoredVariables["APosterioriCovariance"]:
362 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
363 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
364 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
365 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
366 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
367 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
368 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
369 self.StoredVariables["APosterioriCorrelations"].store( _C )
371 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))
372 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))
373 logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
374 logging.debug("%s Terminé", self._name)
377 def get(self, key=None):
379 Renvoie l'une des variables stockées identifiée par la clé, ou le
380 dictionnaire de l'ensemble des variables disponibles en l'absence de
381 clé. Ce sont directement les variables sous forme objet qui sont
382 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
383 des classes de persistance.
386 return self.StoredVariables[key]
388 return self.StoredVariables
390 def __contains__(self, key=None):
391 "D.__contains__(k) -> True if D has a key k, else False"
392 return key in self.StoredVariables
395 "D.keys() -> list of D's keys"
396 return self.StoredVariables.keys()
398 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
400 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
401 sa forme mathématique la plus naturelle possible.
403 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
405 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
407 Permet de définir dans l'algorithme des paramètres requis et leurs
408 caractéristiques par défaut.
411 raise ValueError("A name is mandatory to define a required parameter.")
413 self.__required_parameters[name] = {
415 "typecast" : typecast,
421 logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
423 def getRequiredParameters(self, noDetails=True):
425 Renvoie la liste des noms de paramètres requis ou directement le
426 dictionnaire des paramètres requis.
429 ks = self.__required_parameters.keys()
433 return self.__required_parameters
435 def setParameterValue(self, name=None, value=None):
437 Renvoie la valeur d'un paramètre requis de manière contrôlée
439 default = self.__required_parameters[name]["default"]
440 typecast = self.__required_parameters[name]["typecast"]
441 minval = self.__required_parameters[name]["minval"]
442 maxval = self.__required_parameters[name]["maxval"]
443 listval = self.__required_parameters[name]["listval"]
445 if value is None and default is None:
447 elif value is None and default is not None:
448 if typecast is None: __val = default
449 else: __val = typecast( default )
451 if typecast is None: __val = value
452 else: __val = typecast( value )
454 if minval is not None and (numpy.array(__val, float) < minval).any():
455 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
456 if maxval is not None and (numpy.array(__val, float) > maxval).any():
457 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
458 if listval is not None:
459 if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
462 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
463 elif __val not in listval:
464 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
467 def setParameters(self, fromDico={}):
469 Permet de stocker les paramètres reçus dans le dictionnaire interne.
471 self._parameters.update( fromDico )
472 for k in self.__required_parameters.keys():
473 if k in fromDico.keys():
474 self._parameters[k] = self.setParameterValue(k,fromDico[k])
476 self._parameters[k] = self.setParameterValue(k)
477 logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
479 # ==============================================================================
480 class Diagnostic(object):
482 Classe générale d'interface de type diagnostic
484 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
485 même temps que l'une des classes de persistance, comme "OneScalar" par
488 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
489 méthode "_formula" pour écrire explicitement et proprement la formule pour
490 l'écriture mathématique du calcul du diagnostic (méthode interne non
491 publique), et "calculate" pour activer la précédente tout en ayant vérifié
492 et préparé les données, et pour stocker les résultats à chaque pas (méthode
493 externe d'activation).
495 def __init__(self, name = "", parameters = {}):
497 self.name = str(name)
498 self.parameters = dict( parameters )
500 def _formula(self, *args):
502 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
503 mathématique la plus naturelle possible.
505 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
507 def calculate(self, *args):
509 Active la formule de calcul avec les arguments correctement rangés
511 raise NotImplementedError("Diagnostic activation method has not been implemented!")
513 # ==============================================================================
514 class Covariance(object):
516 Classe générale d'interface de type covariance
519 name = "GenericCovariance",
521 asEyeByScalar = None,
522 asEyeByVector = None,
527 Permet de définir une covariance :
528 - asCovariance : entrée des données, comme une matrice compatible avec
529 le constructeur de numpy.matrix
530 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
531 multiplicatif d'une matrice de corrélation identité, aucune matrice
532 n'étant donc explicitement à donner
533 - asEyeByVector : entrée des données comme un seul vecteur de variance,
534 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
535 n'étant donc explicitement à donner
536 - asCovObject : entrée des données comme un objet python, qui a les
537 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
538 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
539 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
540 - toBeChecked : booléen indiquant si le caractère SDP de la matrice
541 pleine doit être vérifié
543 self.__name = str(name)
544 self.__check = bool(toBeChecked)
547 self.__is_scalar = False
548 self.__is_vector = False
549 self.__is_matrix = False
550 self.__is_object = False
551 if asEyeByScalar is not None:
552 self.__is_scalar = True
553 self.__C = numpy.abs( float(asEyeByScalar) )
556 elif asEyeByVector is not None:
557 self.__is_vector = True
558 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
559 self.shape = (self.__C.size,self.__C.size)
560 self.size = self.__C.size**2
561 elif asCovariance is not None:
562 self.__is_matrix = True
563 self.__C = numpy.matrix( asCovariance, float )
564 self.shape = self.__C.shape
565 self.size = self.__C.size
566 elif asCovObject is not None:
567 self.__is_object = True
568 self.__C = asCovObject
569 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
570 if not hasattr(self.__C,at):
571 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
572 if hasattr(self.__C,"shape"):
573 self.shape = self.__C.shape
576 if hasattr(self.__C,"size"):
577 self.size = self.__C.size
582 # 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)
586 def __validate(self):
588 if self.ismatrix() and min(self.shape) != max(self.shape):
589 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))
590 if self.isobject() and min(self.shape) != max(self.shape):
591 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))
592 if self.isscalar() and self.__C <= 0:
593 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
594 if self.isvector() and (self.__C <= 0).any():
595 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
596 if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
598 L = numpy.linalg.cholesky( self.__C )
600 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
601 if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
603 L = self.__C.cholesky()
605 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
608 "Vérification du type interne"
609 return self.__is_scalar
612 "Vérification du type interne"
613 return self.__is_vector
616 "Vérification du type interne"
617 return self.__is_matrix
620 "Vérification du type interne"
621 return self.__is_object
626 return Covariance(self.__name+"I", asCovariance = self.__C.I )
627 elif self.isvector():
628 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
629 elif self.isscalar():
630 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
631 elif self.isobject():
632 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
634 return None # Indispensable
639 return Covariance(self.__name+"T", asCovariance = self.__C.T )
640 elif self.isvector():
641 return Covariance(self.__name+"T", asEyeByVector = self.__C )
642 elif self.isscalar():
643 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
644 elif self.isobject():
645 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
648 "Décomposition de Cholesky"
650 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
651 elif self.isvector():
652 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
653 elif self.isscalar():
654 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
655 elif self.isobject() and hasattr(self.__C,"cholesky"):
656 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
659 "Inversion de la décomposition de Cholesky"
661 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
662 elif self.isvector():
663 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
664 elif self.isscalar():
665 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
666 elif self.isobject() and hasattr(self.__C,"choleskyI"):
667 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
669 def diag(self, msize=None):
670 "Diagonale de la matrice"
672 return numpy.diag(self.__C)
673 elif self.isvector():
675 elif self.isscalar():
677 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,))
679 return self.__C * numpy.ones(int(msize))
680 elif self.isobject():
681 return self.__C.diag()
683 def asfullmatrix(self, msize=None):
687 elif self.isvector():
688 return numpy.matrix( numpy.diag(self.__C), float )
689 elif self.isscalar():
691 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,))
693 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
694 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
695 return self.__C.asfullmatrix()
697 def trace(self, msize=None):
698 "Trace de la matrice"
700 return numpy.trace(self.__C)
701 elif self.isvector():
702 return float(numpy.sum(self.__C))
703 elif self.isscalar():
705 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,))
707 return self.__C * int(msize)
708 elif self.isobject():
709 return self.__C.trace()
712 "x.__repr__() <==> repr(x)"
713 return repr(self.__C)
716 "x.__str__() <==> str(x)"
719 def __add__(self, other):
720 "x.__add__(y) <==> x+y"
721 if self.ismatrix() or self.isobject():
722 return self.__C + numpy.asmatrix(other)
723 elif self.isvector() or self.isscalar():
724 _A = numpy.asarray(other)
725 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
726 return numpy.asmatrix(_A)
728 def __radd__(self, other):
729 "x.__radd__(y) <==> y+x"
730 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
732 def __sub__(self, other):
733 "x.__sub__(y) <==> x-y"
734 if self.ismatrix() or self.isobject():
735 return self.__C - numpy.asmatrix(other)
736 elif self.isvector() or self.isscalar():
737 _A = numpy.asarray(other)
738 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
739 return numpy.asmatrix(_A)
741 def __rsub__(self, other):
742 "x.__rsub__(y) <==> y-x"
743 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
746 "x.__neg__() <==> -x"
749 def __mul__(self, other):
750 "x.__mul__(y) <==> x*y"
751 if self.ismatrix() and isinstance(other,numpy.matrix):
752 return self.__C * other
753 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
754 or isinstance(other,list) \
755 or isinstance(other,tuple)):
756 if numpy.ravel(other).size == self.shape[1]: # Vecteur
757 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
758 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
759 return self.__C * numpy.asmatrix(other)
761 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
762 elif self.isvector() and (isinstance(other,numpy.matrix) \
763 or isinstance(other,numpy.ndarray) \
764 or isinstance(other,list) \
765 or isinstance(other,tuple)):
766 if numpy.ravel(other).size == self.shape[1]: # Vecteur
767 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
768 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
769 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
771 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
772 elif self.isscalar() and isinstance(other,numpy.matrix):
773 return self.__C * other
774 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
775 or isinstance(other,list) \
776 or isinstance(other,tuple)):
777 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
778 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
780 return self.__C * numpy.asmatrix(other)
781 elif self.isobject():
782 return self.__C.__mul__(other)
784 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
786 def __rmul__(self, other):
787 "x.__rmul__(y) <==> y*x"
788 if self.ismatrix() and isinstance(other,numpy.matrix):
789 return other * self.__C
790 elif self.isvector() and isinstance(other,numpy.matrix):
791 if numpy.ravel(other).size == self.shape[0]: # Vecteur
792 return numpy.asmatrix(numpy.ravel(other) * self.__C)
793 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
794 return numpy.asmatrix(numpy.array(other) * self.__C)
796 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
797 elif self.isscalar() and isinstance(other,numpy.matrix):
798 return other * self.__C
799 elif self.isobject():
800 return self.__C.__rmul__(other)
802 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
805 "x.__len__() <==> len(x)"
808 # ==============================================================================
809 def CostFunction3D(_x,
810 _Hm = None, # Pour simuler Hm(x) : HO["Direct"].appliedTo
811 _HmX = None, # Simulation déjà faite de Hm(x)
812 _arg = None, # Arguments supplementaires pour Hm, sous la forme d'un tuple
817 _SIV = False, # A résorber pour la 8.0
818 _SSC = [], # self._parameters["StoreSupplementaryCalculations"]
819 _nPS = 0, # nbPreviousSteps
820 _QM = "DA", # QualityMeasure
821 _SSV = {}, # Entrée et/ou sortie : self.StoredVariables
822 _fRt = False, # Restitue ou pas la sortie étendue
823 _sSc = True, # Stocke ou pas les SSC
826 Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
827 et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
828 DFO, QuantileRegression
834 for k in ["CostFunctionJ",
840 "SimulatedObservationAtCurrentOptimum",
841 "SimulatedObservationAtCurrentState",
845 if hasattr(_SSV[k],"store"):
846 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
848 _X = numpy.asmatrix(numpy.ravel( _x )).T
849 if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
850 _SSV["CurrentState"].append( _X )
856 raise ValueError("COSTFUNCTION3D Operator has to be defined.")
860 _HX = _Hm( _X, *_arg )
861 _HX = numpy.asmatrix(numpy.ravel( _HX )).T
863 if "SimulatedObservationAtCurrentState" in _SSC or \
864 "SimulatedObservationAtCurrentOptimum" in _SSC:
865 _SSV["SimulatedObservationAtCurrentState"].append( _HX )
867 if numpy.any(numpy.isnan(_HX)):
868 Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
870 _Y = numpy.asmatrix(numpy.ravel( _Y )).T
871 if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
872 if _BI is None or _RI is None:
873 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
874 _Xb = numpy.asmatrix(numpy.ravel( _Xb )).T
875 Jb = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
876 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
877 elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
879 raise ValueError("Observation error covariance matrix has to be properly defined!")
881 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
882 elif _QM in ["LeastSquares", "LS", "L2"]:
884 Jo = 0.5 * (_Y - _HX).T * (_Y - _HX)
885 elif _QM in ["AbsoluteValue", "L1"]:
887 Jo = numpy.sum( numpy.abs(_Y - _HX) )
888 elif _QM in ["MaximumError", "ME"]:
890 Jo = numpy.max( numpy.abs(_Y - _HX) )
891 elif _QM in ["QR", "Null"]:
895 raise ValueError("Unknown asked quality measure!")
897 J = float( Jb ) + float( Jo )
900 _SSV["CostFunctionJb"].append( Jb )
901 _SSV["CostFunctionJo"].append( Jo )
902 _SSV["CostFunctionJ" ].append( J )
904 if "IndexOfOptimum" in _SSC or \
905 "CurrentOptimum" in _SSC or \
906 "SimulatedObservationAtCurrentOptimum" in _SSC:
907 IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
908 if "IndexOfOptimum" in _SSC:
909 _SSV["IndexOfOptimum"].append( IndexMin )
910 if "CurrentOptimum" in _SSC:
911 _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
912 if "SimulatedObservationAtCurrentOptimum" in _SSC:
913 _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
918 if _QM in ["QR"]: # Pour le QuantileRegression
923 # ==============================================================================
924 if __name__ == "__main__":
925 print '\n AUTODIAGNOSTIC \n'