1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2016 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["CostFunctionJAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJAtCurrentOptimum")
326 self.StoredVariables["CostFunctionJbAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJbAtCurrentOptimum")
327 self.StoredVariables["CostFunctionJoAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJoAtCurrentOptimum")
328 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
329 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
330 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
331 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
332 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
333 self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum")
334 self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
335 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
336 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
337 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
338 self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
339 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
340 self.StoredVariables["InnovationAtCurrentState"] = Persistence.OneVector(name = "InnovationAtCurrentState")
341 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
342 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
343 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
344 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
345 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
346 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
347 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
348 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
349 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
350 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
351 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
352 self.StoredVariables["Residu"] = Persistence.OneScalar(name = "Residu")
356 logging.debug("%s Lancement", self._name)
357 logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
360 def _post_run(self,_oH=None):
362 if ("StoreSupplementaryCalculations" in self._parameters) and \
363 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
364 for _A in self.StoredVariables["APosterioriCovariance"]:
365 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
366 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
367 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
368 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
369 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
370 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
371 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
372 self.StoredVariables["APosterioriCorrelations"].store( _C )
374 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))
375 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))
376 logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
377 logging.debug("%s Terminé", self._name)
380 def get(self, key=None):
382 Renvoie l'une des variables stockées identifiée par la clé, ou le
383 dictionnaire de l'ensemble des variables disponibles en l'absence de
384 clé. Ce sont directement les variables sous forme objet qui sont
385 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
386 des classes de persistance.
389 return self.StoredVariables[key]
391 return self.StoredVariables
393 def __contains__(self, key=None):
394 "D.__contains__(k) -> True if D has a key k, else False"
395 return key in self.StoredVariables
398 "D.keys() -> list of D's keys"
399 return self.StoredVariables.keys()
401 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
403 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
404 sa forme mathématique la plus naturelle possible.
406 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
408 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
410 Permet de définir dans l'algorithme des paramètres requis et leurs
411 caractéristiques par défaut.
414 raise ValueError("A name is mandatory to define a required parameter.")
416 self.__required_parameters[name] = {
418 "typecast" : typecast,
424 logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
426 def getRequiredParameters(self, noDetails=True):
428 Renvoie la liste des noms de paramètres requis ou directement le
429 dictionnaire des paramètres requis.
432 ks = self.__required_parameters.keys()
436 return self.__required_parameters
438 def setParameterValue(self, name=None, value=None):
440 Renvoie la valeur d'un paramètre requis de manière contrôlée
442 default = self.__required_parameters[name]["default"]
443 typecast = self.__required_parameters[name]["typecast"]
444 minval = self.__required_parameters[name]["minval"]
445 maxval = self.__required_parameters[name]["maxval"]
446 listval = self.__required_parameters[name]["listval"]
448 if value is None and default is None:
450 elif value is None and default is not None:
451 if typecast is None: __val = default
452 else: __val = typecast( default )
454 if typecast is None: __val = value
455 else: __val = typecast( value )
457 if minval is not None and (numpy.array(__val, float) < minval).any():
458 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
459 if maxval is not None and (numpy.array(__val, float) > maxval).any():
460 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
461 if listval is not None:
462 if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
465 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
466 elif __val not in listval:
467 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
470 def setParameters(self, fromDico={}):
472 Permet de stocker les paramètres reçus dans le dictionnaire interne.
474 self._parameters.update( fromDico )
475 for k in self.__required_parameters.keys():
476 if k in fromDico.keys():
477 self._parameters[k] = self.setParameterValue(k,fromDico[k])
479 self._parameters[k] = self.setParameterValue(k)
480 logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
482 # ==============================================================================
483 class Diagnostic(object):
485 Classe générale d'interface de type diagnostic
487 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
488 même temps que l'une des classes de persistance, comme "OneScalar" par
491 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
492 méthode "_formula" pour écrire explicitement et proprement la formule pour
493 l'écriture mathématique du calcul du diagnostic (méthode interne non
494 publique), et "calculate" pour activer la précédente tout en ayant vérifié
495 et préparé les données, et pour stocker les résultats à chaque pas (méthode
496 externe d'activation).
498 def __init__(self, name = "", parameters = {}):
500 self.name = str(name)
501 self.parameters = dict( parameters )
503 def _formula(self, *args):
505 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
506 mathématique la plus naturelle possible.
508 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
510 def calculate(self, *args):
512 Active la formule de calcul avec les arguments correctement rangés
514 raise NotImplementedError("Diagnostic activation method has not been implemented!")
516 # ==============================================================================
517 class Covariance(object):
519 Classe générale d'interface de type covariance
522 name = "GenericCovariance",
524 asEyeByScalar = None,
525 asEyeByVector = None,
530 Permet de définir une covariance :
531 - asCovariance : entrée des données, comme une matrice compatible avec
532 le constructeur de numpy.matrix
533 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
534 multiplicatif d'une matrice de corrélation identité, aucune matrice
535 n'étant donc explicitement à donner
536 - asEyeByVector : entrée des données comme un seul vecteur de variance,
537 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
538 n'étant donc explicitement à donner
539 - asCovObject : entrée des données comme un objet python, qui a les
540 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
541 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
542 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
543 - toBeChecked : booléen indiquant si le caractère SDP de la matrice
544 pleine doit être vérifié
546 self.__name = str(name)
547 self.__check = bool(toBeChecked)
550 self.__is_scalar = False
551 self.__is_vector = False
552 self.__is_matrix = False
553 self.__is_object = False
554 if asEyeByScalar is not None:
555 self.__is_scalar = True
556 self.__C = numpy.abs( float(asEyeByScalar) )
559 elif asEyeByVector is not None:
560 self.__is_vector = True
561 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
562 self.shape = (self.__C.size,self.__C.size)
563 self.size = self.__C.size**2
564 elif asCovariance is not None:
565 self.__is_matrix = True
566 self.__C = numpy.matrix( asCovariance, float )
567 self.shape = self.__C.shape
568 self.size = self.__C.size
569 elif asCovObject is not None:
570 self.__is_object = True
571 self.__C = asCovObject
572 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
573 if not hasattr(self.__C,at):
574 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
575 if hasattr(self.__C,"shape"):
576 self.shape = self.__C.shape
579 if hasattr(self.__C,"size"):
580 self.size = self.__C.size
585 # 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)
589 def __validate(self):
591 if self.ismatrix() and min(self.shape) != max(self.shape):
592 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))
593 if self.isobject() and min(self.shape) != max(self.shape):
594 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))
595 if self.isscalar() and self.__C <= 0:
596 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
597 if self.isvector() and (self.__C <= 0).any():
598 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
599 if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
601 L = numpy.linalg.cholesky( self.__C )
603 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
604 if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
606 L = self.__C.cholesky()
608 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
611 "Vérification du type interne"
612 return self.__is_scalar
615 "Vérification du type interne"
616 return self.__is_vector
619 "Vérification du type interne"
620 return self.__is_matrix
623 "Vérification du type interne"
624 return self.__is_object
629 return Covariance(self.__name+"I", asCovariance = self.__C.I )
630 elif self.isvector():
631 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
632 elif self.isscalar():
633 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
634 elif self.isobject():
635 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
637 return None # Indispensable
642 return Covariance(self.__name+"T", asCovariance = self.__C.T )
643 elif self.isvector():
644 return Covariance(self.__name+"T", asEyeByVector = self.__C )
645 elif self.isscalar():
646 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
647 elif self.isobject():
648 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
651 "Décomposition de Cholesky"
653 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
654 elif self.isvector():
655 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
656 elif self.isscalar():
657 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
658 elif self.isobject() and hasattr(self.__C,"cholesky"):
659 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
662 "Inversion de la décomposition de Cholesky"
664 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
665 elif self.isvector():
666 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
667 elif self.isscalar():
668 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
669 elif self.isobject() and hasattr(self.__C,"choleskyI"):
670 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
672 def diag(self, msize=None):
673 "Diagonale de la matrice"
675 return numpy.diag(self.__C)
676 elif self.isvector():
678 elif self.isscalar():
680 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,))
682 return self.__C * numpy.ones(int(msize))
683 elif self.isobject():
684 return self.__C.diag()
686 def asfullmatrix(self, msize=None):
690 elif self.isvector():
691 return numpy.matrix( numpy.diag(self.__C), float )
692 elif self.isscalar():
694 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,))
696 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
697 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
698 return self.__C.asfullmatrix()
700 def trace(self, msize=None):
701 "Trace de la matrice"
703 return numpy.trace(self.__C)
704 elif self.isvector():
705 return float(numpy.sum(self.__C))
706 elif self.isscalar():
708 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,))
710 return self.__C * int(msize)
711 elif self.isobject():
712 return self.__C.trace()
715 "x.__repr__() <==> repr(x)"
716 return repr(self.__C)
719 "x.__str__() <==> str(x)"
722 def __add__(self, other):
723 "x.__add__(y) <==> x+y"
724 if self.ismatrix() or self.isobject():
725 return self.__C + numpy.asmatrix(other)
726 elif self.isvector() or self.isscalar():
727 _A = numpy.asarray(other)
728 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
729 return numpy.asmatrix(_A)
731 def __radd__(self, other):
732 "x.__radd__(y) <==> y+x"
733 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
735 def __sub__(self, other):
736 "x.__sub__(y) <==> x-y"
737 if self.ismatrix() or self.isobject():
738 return self.__C - numpy.asmatrix(other)
739 elif self.isvector() or self.isscalar():
740 _A = numpy.asarray(other)
741 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
742 return numpy.asmatrix(_A)
744 def __rsub__(self, other):
745 "x.__rsub__(y) <==> y-x"
746 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
749 "x.__neg__() <==> -x"
752 def __mul__(self, other):
753 "x.__mul__(y) <==> x*y"
754 if self.ismatrix() and isinstance(other,numpy.matrix):
755 return self.__C * other
756 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
757 or isinstance(other,list) \
758 or isinstance(other,tuple)):
759 if numpy.ravel(other).size == self.shape[1]: # Vecteur
760 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
761 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
762 return self.__C * numpy.asmatrix(other)
764 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
765 elif self.isvector() and (isinstance(other,numpy.matrix) \
766 or isinstance(other,numpy.ndarray) \
767 or isinstance(other,list) \
768 or isinstance(other,tuple)):
769 if numpy.ravel(other).size == self.shape[1]: # Vecteur
770 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
771 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
772 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
774 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
775 elif self.isscalar() and isinstance(other,numpy.matrix):
776 return self.__C * other
777 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
778 or isinstance(other,list) \
779 or isinstance(other,tuple)):
780 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
781 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
783 return self.__C * numpy.asmatrix(other)
784 elif self.isobject():
785 return self.__C.__mul__(other)
787 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
789 def __rmul__(self, other):
790 "x.__rmul__(y) <==> y*x"
791 if self.ismatrix() and isinstance(other,numpy.matrix):
792 return other * self.__C
793 elif self.isvector() and isinstance(other,numpy.matrix):
794 if numpy.ravel(other).size == self.shape[0]: # Vecteur
795 return numpy.asmatrix(numpy.ravel(other) * self.__C)
796 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
797 return numpy.asmatrix(numpy.array(other) * self.__C)
799 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
800 elif self.isscalar() and isinstance(other,numpy.matrix):
801 return other * self.__C
802 elif self.isobject():
803 return self.__C.__rmul__(other)
805 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
808 "x.__len__() <==> len(x)"
811 # ==============================================================================
812 def CostFunction3D(_x,
813 _Hm = None, # Pour simuler Hm(x) : HO["Direct"].appliedTo
814 _HmX = None, # Simulation déjà faite de Hm(x)
815 _arg = None, # Arguments supplementaires pour Hm, sous la forme d'un tuple
820 _SIV = False, # A résorber pour la 8.0
821 _SSC = [], # self._parameters["StoreSupplementaryCalculations"]
822 _nPS = 0, # nbPreviousSteps
823 _QM = "DA", # QualityMeasure
824 _SSV = {}, # Entrée et/ou sortie : self.StoredVariables
825 _fRt = False, # Restitue ou pas la sortie étendue
826 _sSc = True, # Stocke ou pas les SSC
829 Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
830 et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
831 DFO, QuantileRegression
837 for k in ["CostFunctionJ",
843 "SimulatedObservationAtCurrentOptimum",
844 "SimulatedObservationAtCurrentState",
848 if hasattr(_SSV[k],"store"):
849 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
851 _X = numpy.asmatrix(numpy.ravel( _x )).T
852 if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
853 _SSV["CurrentState"].append( _X )
859 raise ValueError("COSTFUNCTION3D Operator has to be defined.")
863 _HX = _Hm( _X, *_arg )
864 _HX = numpy.asmatrix(numpy.ravel( _HX )).T
866 if "SimulatedObservationAtCurrentState" in _SSC or \
867 "SimulatedObservationAtCurrentOptimum" in _SSC:
868 _SSV["SimulatedObservationAtCurrentState"].append( _HX )
870 if numpy.any(numpy.isnan(_HX)):
871 Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
873 _Y = numpy.asmatrix(numpy.ravel( _Y )).T
874 if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
875 if _BI is None or _RI is None:
876 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
877 _Xb = numpy.asmatrix(numpy.ravel( _Xb )).T
878 Jb = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
879 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
880 elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
882 raise ValueError("Observation error covariance matrix has to be properly defined!")
884 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
885 elif _QM in ["LeastSquares", "LS", "L2"]:
887 Jo = 0.5 * (_Y - _HX).T * (_Y - _HX)
888 elif _QM in ["AbsoluteValue", "L1"]:
890 Jo = numpy.sum( numpy.abs(_Y - _HX) )
891 elif _QM in ["MaximumError", "ME"]:
893 Jo = numpy.max( numpy.abs(_Y - _HX) )
894 elif _QM in ["QR", "Null"]:
898 raise ValueError("Unknown asked quality measure!")
900 J = float( Jb ) + float( Jo )
903 _SSV["CostFunctionJb"].append( Jb )
904 _SSV["CostFunctionJo"].append( Jo )
905 _SSV["CostFunctionJ" ].append( J )
907 if "IndexOfOptimum" in _SSC or \
908 "CurrentOptimum" in _SSC or \
909 "SimulatedObservationAtCurrentOptimum" in _SSC:
910 IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
911 if "IndexOfOptimum" in _SSC:
912 _SSV["IndexOfOptimum"].append( IndexMin )
913 if "CurrentOptimum" in _SSC:
914 _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
915 if "SimulatedObservationAtCurrentOptimum" in _SSC:
916 _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
921 if _QM in ["QR"]: # Pour le QuantileRegression
926 # ==============================================================================
927 if __name__ == "__main__":
928 print '\n AUTODIAGNOSTIC \n'