1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2017 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, HValue = None):
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
153 if HValue is not None:
154 HxValue = numpy.asmatrix( numpy.ravel( HValue ) ).T
156 Operator.CM.storeValueInX(xValue,HxValue)
159 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
161 __alreadyCalculated = False
163 if __alreadyCalculated:
164 self.__addOneCacheCall()
167 if self.__Matrix is not None:
168 self.__addOneMatrixCall()
169 HxValue = self.__Matrix * xValue
171 self.__addOneMethodCall()
172 HxValue = self.__Method( xValue )
174 Operator.CM.storeValueInX(xValue,HxValue)
178 def appliedControledFormTo(self, (xValue, uValue) ):
180 Permet de restituer le résultat de l'application de l'opérateur à une
181 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
182 argument devant a priori être du bon type. Si la uValue est None,
183 on suppose que l'opérateur ne s'applique qu'à xValue.
185 - xValue : argument X adapté pour appliquer l'opérateur
186 - uValue : argument U adapté pour appliquer l'opérateur
188 if self.__Matrix is not None:
189 self.__addOneMatrixCall()
190 return self.__Matrix * xValue
191 elif uValue is not None:
192 self.__addOneMethodCall()
193 return self.__Method( (xValue, uValue) )
195 self.__addOneMethodCall()
196 return self.__Method( xValue )
198 def appliedInXTo(self, (xNominal, xValue) ):
200 Permet de restituer le résultat de l'application de l'opérateur à un
201 argument xValue, sachant que l'opérateur est valable en xNominal.
202 Cette méthode se contente d'appliquer, son argument devant a priori
203 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
204 alors il est valable en tout point nominal et il n'est pas nécessaire
206 Arguments : une liste contenant
207 - xNominal : argument permettant de donner le point où l'opérateur
208 est construit pour etre ensuite appliqué
209 - xValue : argument adapté pour appliquer l'opérateur
211 if self.__Matrix is not None:
212 self.__addOneMatrixCall()
213 return self.__Matrix * xValue
215 self.__addOneMethodCall()
216 return self.__Method( (xNominal, xValue) )
218 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
220 Permet de renvoyer l'opérateur sous la forme d'une matrice
222 if self.__Matrix is not None:
223 self.__addOneMatrixCall()
225 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
226 self.__addOneMethodCall()
227 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
229 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
233 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
234 la forme d'une matrice
236 if self.__Matrix is not None:
237 return self.__Matrix.shape
239 raise ValueError("Matrix form of the operator is not available, nor the shape")
241 def nbcalls(self, which=None):
243 Renvoie les nombres d'évaluations de l'opérateur
246 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
247 self.__NbCallsAsMatrix,
248 self.__NbCallsAsMethod,
249 self.__NbCallsOfCached,
250 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
251 Operator.NbCallsAsMatrix,
252 Operator.NbCallsAsMethod,
253 Operator.NbCallsOfCached,
255 if which is None: return __nbcalls
256 else: return __nbcalls[which]
258 def __addOneMatrixCall(self):
259 "Comptabilise un appel"
260 self.__NbCallsAsMatrix += 1 # Decompte local
261 Operator.NbCallsAsMatrix += 1 # Decompte global
263 def __addOneMethodCall(self):
264 "Comptabilise un appel"
265 self.__NbCallsAsMethod += 1 # Decompte local
266 Operator.NbCallsAsMethod += 1 # Decompte global
268 def __addOneCacheCall(self):
269 "Comptabilise un appel"
270 self.__NbCallsOfCached += 1 # Decompte local
271 Operator.NbCallsOfCached += 1 # Decompte global
273 # ==============================================================================
274 class Algorithm(object):
276 Classe générale d'interface de type algorithme
278 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
279 d'assimilation, en fournissant un container (dictionnaire) de variables
280 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
282 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
284 def __init__(self, name):
286 L'initialisation présente permet de fabriquer des variables de stockage
287 disponibles de manière générique dans les algorithmes élémentaires. Ces
288 variables de stockage sont ensuite conservées dans un dictionnaire
289 interne à l'objet, mais auquel on accède par la méthode "get".
291 Les variables prévues sont :
292 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
293 - CostFunctionJb : partie ébauche ou background de la fonction-cout
294 - CostFunctionJo : partie observations de la fonction-cout
295 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
296 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
297 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
298 - CurrentState : état courant lors d'itérations
299 - Analysis : l'analyse Xa
300 - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
301 - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
302 - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
303 - Innovation : l'innovation : d = Y - H(X)
304 - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
305 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
306 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
307 - MahalanobisConsistency : indicateur de consistance des covariances
308 - OMA : Observation moins Analysis : Y - Xa
309 - OMB : Observation moins Background : Y - Xb
310 - AMB : Analysis moins Background : Xa - Xb
311 - APosterioriCovariance : matrice A
312 - APosterioriVariances : variances de la matrice A
313 - APosterioriStandardDeviations : écart-types de la matrice A
314 - APosterioriCorrelations : correlations de la matrice A
315 - Residu : dans le cas des algorithmes de vérification
316 On peut rajouter des variables à stocker dans l'initialisation de
317 l'algorithme élémentaire qui va hériter de cette classe
319 logging.debug("%s Initialisation", str(name))
320 self._m = PlatformInfo.SystemUsage()
322 self._name = str( name )
323 self._parameters = {"StoreSupplementaryCalculations":[]}
324 self.__required_parameters = {}
325 self.StoredVariables = {}
327 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
328 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
329 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
330 self.StoredVariables["CostFunctionJAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJAtCurrentOptimum")
331 self.StoredVariables["CostFunctionJbAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJbAtCurrentOptimum")
332 self.StoredVariables["CostFunctionJoAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJoAtCurrentOptimum")
333 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
334 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
335 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
336 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
337 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
338 self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum")
339 self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
340 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
341 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
342 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
343 self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
344 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
345 self.StoredVariables["InnovationAtCurrentState"] = Persistence.OneVector(name = "InnovationAtCurrentState")
346 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
347 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
348 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
349 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
350 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
351 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
352 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
353 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
354 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
355 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
356 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
357 self.StoredVariables["Residu"] = Persistence.OneScalar(name = "Residu")
361 logging.debug("%s Lancement", self._name)
362 logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
365 def _post_run(self,_oH=None):
367 if ("StoreSupplementaryCalculations" in self._parameters) and \
368 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
369 for _A in self.StoredVariables["APosterioriCovariance"]:
370 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
371 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
372 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
373 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
374 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
375 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
376 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
377 self.StoredVariables["APosterioriCorrelations"].store( _C )
379 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))
380 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))
381 logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
382 logging.debug("%s Terminé", self._name)
385 def get(self, key=None):
387 Renvoie l'une des variables stockées identifiée par la clé, ou le
388 dictionnaire de l'ensemble des variables disponibles en l'absence de
389 clé. Ce sont directement les variables sous forme objet qui sont
390 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
391 des classes de persistance.
394 return self.StoredVariables[key]
396 return self.StoredVariables
398 def __contains__(self, key=None):
399 "D.__contains__(k) -> True if D has a key k, else False"
400 return key in self.StoredVariables
403 "D.keys() -> list of D's keys"
404 return self.StoredVariables.keys()
406 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
408 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
409 sa forme mathématique la plus naturelle possible.
411 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
413 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
415 Permet de définir dans l'algorithme des paramètres requis et leurs
416 caractéristiques par défaut.
419 raise ValueError("A name is mandatory to define a required parameter.")
421 self.__required_parameters[name] = {
423 "typecast" : typecast,
429 logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
431 def getRequiredParameters(self, noDetails=True):
433 Renvoie la liste des noms de paramètres requis ou directement le
434 dictionnaire des paramètres requis.
437 ks = self.__required_parameters.keys()
441 return self.__required_parameters
443 def setParameterValue(self, name=None, value=None):
445 Renvoie la valeur d'un paramètre requis de manière contrôlée
447 default = self.__required_parameters[name]["default"]
448 typecast = self.__required_parameters[name]["typecast"]
449 minval = self.__required_parameters[name]["minval"]
450 maxval = self.__required_parameters[name]["maxval"]
451 listval = self.__required_parameters[name]["listval"]
453 if value is None and default is None:
455 elif value is None and default is not None:
456 if typecast is None: __val = default
457 else: __val = typecast( default )
459 if typecast is None: __val = value
460 else: __val = typecast( value )
462 if minval is not None and (numpy.array(__val, float) < minval).any():
463 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
464 if maxval is not None and (numpy.array(__val, float) > maxval).any():
465 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
466 if listval is not None:
467 if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
470 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
471 elif __val not in listval:
472 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
475 def setParameters(self, fromDico={}):
477 Permet de stocker les paramètres reçus dans le dictionnaire interne.
479 self._parameters.update( fromDico )
480 for k in self.__required_parameters.keys():
481 if k in fromDico.keys():
482 self._parameters[k] = self.setParameterValue(k,fromDico[k])
484 self._parameters[k] = self.setParameterValue(k)
485 logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
487 # ==============================================================================
488 class Diagnostic(object):
490 Classe générale d'interface de type diagnostic
492 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
493 même temps que l'une des classes de persistance, comme "OneScalar" par
496 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
497 méthode "_formula" pour écrire explicitement et proprement la formule pour
498 l'écriture mathématique du calcul du diagnostic (méthode interne non
499 publique), et "calculate" pour activer la précédente tout en ayant vérifié
500 et préparé les données, et pour stocker les résultats à chaque pas (méthode
501 externe d'activation).
503 def __init__(self, name = "", parameters = {}):
505 self.name = str(name)
506 self.parameters = dict( parameters )
508 def _formula(self, *args):
510 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
511 mathématique la plus naturelle possible.
513 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
515 def calculate(self, *args):
517 Active la formule de calcul avec les arguments correctement rangés
519 raise NotImplementedError("Diagnostic activation method has not been implemented!")
521 # ==============================================================================
522 class Covariance(object):
524 Classe générale d'interface de type covariance
527 name = "GenericCovariance",
529 asEyeByScalar = None,
530 asEyeByVector = None,
535 Permet de définir une covariance :
536 - asCovariance : entrée des données, comme une matrice compatible avec
537 le constructeur de numpy.matrix
538 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
539 multiplicatif d'une matrice de corrélation identité, aucune matrice
540 n'étant donc explicitement à donner
541 - asEyeByVector : entrée des données comme un seul vecteur de variance,
542 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
543 n'étant donc explicitement à donner
544 - asCovObject : entrée des données comme un objet python, qui a les
545 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
546 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
547 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
548 - toBeChecked : booléen indiquant si le caractère SDP de la matrice
549 pleine doit être vérifié
551 self.__name = str(name)
552 self.__check = bool(toBeChecked)
555 self.__is_scalar = False
556 self.__is_vector = False
557 self.__is_matrix = False
558 self.__is_object = False
559 if asEyeByScalar is not None:
560 self.__is_scalar = True
561 self.__C = numpy.abs( float(asEyeByScalar) )
564 elif asEyeByVector is not None:
565 self.__is_vector = True
566 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
567 self.shape = (self.__C.size,self.__C.size)
568 self.size = self.__C.size**2
569 elif asCovariance is not None:
570 self.__is_matrix = True
571 self.__C = numpy.matrix( asCovariance, float )
572 self.shape = self.__C.shape
573 self.size = self.__C.size
574 elif asCovObject is not None:
575 self.__is_object = True
576 self.__C = asCovObject
577 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
578 if not hasattr(self.__C,at):
579 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
580 if hasattr(self.__C,"shape"):
581 self.shape = self.__C.shape
584 if hasattr(self.__C,"size"):
585 self.size = self.__C.size
590 # 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)
594 def __validate(self):
596 if self.ismatrix() and min(self.shape) != max(self.shape):
597 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))
598 if self.isobject() and min(self.shape) != max(self.shape):
599 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))
600 if self.isscalar() and self.__C <= 0:
601 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
602 if self.isvector() and (self.__C <= 0).any():
603 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
604 if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
606 L = numpy.linalg.cholesky( self.__C )
608 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
609 if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
611 L = self.__C.cholesky()
613 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
616 "Vérification du type interne"
617 return self.__is_scalar
620 "Vérification du type interne"
621 return self.__is_vector
624 "Vérification du type interne"
625 return self.__is_matrix
628 "Vérification du type interne"
629 return self.__is_object
634 return Covariance(self.__name+"I", asCovariance = self.__C.I )
635 elif self.isvector():
636 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
637 elif self.isscalar():
638 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
639 elif self.isobject():
640 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
642 return None # Indispensable
647 return Covariance(self.__name+"T", asCovariance = self.__C.T )
648 elif self.isvector():
649 return Covariance(self.__name+"T", asEyeByVector = self.__C )
650 elif self.isscalar():
651 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
652 elif self.isobject():
653 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
656 "Décomposition de Cholesky"
658 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
659 elif self.isvector():
660 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
661 elif self.isscalar():
662 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
663 elif self.isobject() and hasattr(self.__C,"cholesky"):
664 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
667 "Inversion de la décomposition de Cholesky"
669 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
670 elif self.isvector():
671 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
672 elif self.isscalar():
673 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
674 elif self.isobject() and hasattr(self.__C,"choleskyI"):
675 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
677 def diag(self, msize=None):
678 "Diagonale de la matrice"
680 return numpy.diag(self.__C)
681 elif self.isvector():
683 elif self.isscalar():
685 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,))
687 return self.__C * numpy.ones(int(msize))
688 elif self.isobject():
689 return self.__C.diag()
691 def asfullmatrix(self, msize=None):
695 elif self.isvector():
696 return numpy.matrix( numpy.diag(self.__C), float )
697 elif self.isscalar():
699 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,))
701 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
702 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
703 return self.__C.asfullmatrix()
705 def trace(self, msize=None):
706 "Trace de la matrice"
708 return numpy.trace(self.__C)
709 elif self.isvector():
710 return float(numpy.sum(self.__C))
711 elif self.isscalar():
713 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,))
715 return self.__C * int(msize)
716 elif self.isobject():
717 return self.__C.trace()
720 "x.__repr__() <==> repr(x)"
721 return repr(self.__C)
724 "x.__str__() <==> str(x)"
727 def __add__(self, other):
728 "x.__add__(y) <==> x+y"
729 if self.ismatrix() or self.isobject():
730 return self.__C + numpy.asmatrix(other)
731 elif self.isvector() or self.isscalar():
732 _A = numpy.asarray(other)
733 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
734 return numpy.asmatrix(_A)
736 def __radd__(self, other):
737 "x.__radd__(y) <==> y+x"
738 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
740 def __sub__(self, other):
741 "x.__sub__(y) <==> x-y"
742 if self.ismatrix() or self.isobject():
743 return self.__C - numpy.asmatrix(other)
744 elif self.isvector() or self.isscalar():
745 _A = numpy.asarray(other)
746 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
747 return numpy.asmatrix(_A)
749 def __rsub__(self, other):
750 "x.__rsub__(y) <==> y-x"
751 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
754 "x.__neg__() <==> -x"
757 def __mul__(self, other):
758 "x.__mul__(y) <==> x*y"
759 if self.ismatrix() and isinstance(other,numpy.matrix):
760 return self.__C * other
761 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
762 or isinstance(other,list) \
763 or isinstance(other,tuple)):
764 if numpy.ravel(other).size == self.shape[1]: # Vecteur
765 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
766 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
767 return self.__C * numpy.asmatrix(other)
769 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
770 elif self.isvector() and (isinstance(other,numpy.matrix) \
771 or isinstance(other,numpy.ndarray) \
772 or isinstance(other,list) \
773 or isinstance(other,tuple)):
774 if numpy.ravel(other).size == self.shape[1]: # Vecteur
775 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
776 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
777 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
779 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
780 elif self.isscalar() and isinstance(other,numpy.matrix):
781 return self.__C * other
782 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
783 or isinstance(other,list) \
784 or isinstance(other,tuple)):
785 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
786 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
788 return self.__C * numpy.asmatrix(other)
789 elif self.isobject():
790 return self.__C.__mul__(other)
792 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
794 def __rmul__(self, other):
795 "x.__rmul__(y) <==> y*x"
796 if self.ismatrix() and isinstance(other,numpy.matrix):
797 return other * self.__C
798 elif self.isvector() and isinstance(other,numpy.matrix):
799 if numpy.ravel(other).size == self.shape[0]: # Vecteur
800 return numpy.asmatrix(numpy.ravel(other) * self.__C)
801 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
802 return numpy.asmatrix(numpy.array(other) * self.__C)
804 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
805 elif self.isscalar() and isinstance(other,numpy.matrix):
806 return other * self.__C
807 elif self.isobject():
808 return self.__C.__rmul__(other)
810 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
813 "x.__len__() <==> len(x)"
816 # ==============================================================================
817 def CostFunction3D(_x,
818 _Hm = None, # Pour simuler Hm(x) : HO["Direct"].appliedTo
819 _HmX = None, # Simulation déjà faite de Hm(x)
820 _arg = None, # Arguments supplementaires pour Hm, sous la forme d'un tuple
825 _SIV = False, # A résorber pour la 8.0
826 _SSC = [], # self._parameters["StoreSupplementaryCalculations"]
827 _nPS = 0, # nbPreviousSteps
828 _QM = "DA", # QualityMeasure
829 _SSV = {}, # Entrée et/ou sortie : self.StoredVariables
830 _fRt = False, # Restitue ou pas la sortie étendue
831 _sSc = True, # Stocke ou pas les SSC
834 Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
835 et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
836 DFO, QuantileRegression
842 for k in ["CostFunctionJ",
848 "SimulatedObservationAtCurrentOptimum",
849 "SimulatedObservationAtCurrentState",
853 if hasattr(_SSV[k],"store"):
854 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
856 _X = numpy.asmatrix(numpy.ravel( _x )).T
857 if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
858 _SSV["CurrentState"].append( _X )
864 raise ValueError("COSTFUNCTION3D Operator has to be defined.")
868 _HX = _Hm( _X, *_arg )
869 _HX = numpy.asmatrix(numpy.ravel( _HX )).T
871 if "SimulatedObservationAtCurrentState" in _SSC or \
872 "SimulatedObservationAtCurrentOptimum" in _SSC:
873 _SSV["SimulatedObservationAtCurrentState"].append( _HX )
875 if numpy.any(numpy.isnan(_HX)):
876 Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
878 _Y = numpy.asmatrix(numpy.ravel( _Y )).T
879 if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
880 if _BI is None or _RI is None:
881 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
882 _Xb = numpy.asmatrix(numpy.ravel( _Xb )).T
883 Jb = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
884 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
885 elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
887 raise ValueError("Observation error covariance matrix has to be properly defined!")
889 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
890 elif _QM in ["LeastSquares", "LS", "L2"]:
892 Jo = 0.5 * (_Y - _HX).T * (_Y - _HX)
893 elif _QM in ["AbsoluteValue", "L1"]:
895 Jo = numpy.sum( numpy.abs(_Y - _HX) )
896 elif _QM in ["MaximumError", "ME"]:
898 Jo = numpy.max( numpy.abs(_Y - _HX) )
899 elif _QM in ["QR", "Null"]:
903 raise ValueError("Unknown asked quality measure!")
905 J = float( Jb ) + float( Jo )
908 _SSV["CostFunctionJb"].append( Jb )
909 _SSV["CostFunctionJo"].append( Jo )
910 _SSV["CostFunctionJ" ].append( J )
912 if "IndexOfOptimum" in _SSC or \
913 "CurrentOptimum" in _SSC or \
914 "SimulatedObservationAtCurrentOptimum" in _SSC:
915 IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
916 if "IndexOfOptimum" in _SSC:
917 _SSV["IndexOfOptimum"].append( IndexMin )
918 if "CurrentOptimum" in _SSC:
919 _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
920 if "SimulatedObservationAtCurrentOptimum" in _SSC:
921 _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
926 if _QM in ["QR"]: # Pour le QuantileRegression
931 # ==============================================================================
932 if __name__ == "__main__":
933 print '\n AUTODIAGNOSTIC \n'