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")
359 def _pre_run(self, Parameters ):
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"))
364 # Mise a jour de self._parameters avec Parameters
365 self.__setParameters(Parameters)
367 # Corrections et complements
368 if self._parameters.has_key("Bounds") and (type(self._parameters["Bounds"]) is type([]) or type(self._parameters["Bounds"]) is type(())) and (len(self._parameters["Bounds"]) > 0):
369 logging.debug("%s Prise en compte des bornes effectuee"%(self._name,))
371 self._parameters["Bounds"] = None
373 if logging.getLogger().level < logging.WARNING:
374 self._parameters["optiprint"], self._parameters["optdisp"] = 1, 1
375 if PlatformInfo.has_scipy:
376 import scipy.optimize
377 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_ALL
379 self._parameters["optmessages"] = 15
381 self._parameters["optiprint"], self._parameters["optdisp"] = -1, 0
382 if PlatformInfo.has_scipy:
383 import scipy.optimize
384 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_NONE
386 self._parameters["optmessages"] = 15
390 def _post_run(self,_oH=None):
392 if ("StoreSupplementaryCalculations" in self._parameters) and \
393 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
394 for _A in self.StoredVariables["APosterioriCovariance"]:
395 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
396 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
397 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
398 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
399 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
400 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
401 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
402 self.StoredVariables["APosterioriCorrelations"].store( _C )
404 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))
405 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))
406 logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
407 logging.debug("%s Terminé", self._name)
410 def get(self, key=None):
412 Renvoie l'une des variables stockées identifiée par la clé, ou le
413 dictionnaire de l'ensemble des variables disponibles en l'absence de
414 clé. Ce sont directement les variables sous forme objet qui sont
415 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
416 des classes de persistance.
419 return self.StoredVariables[key]
421 return self.StoredVariables
423 def __contains__(self, key=None):
424 "D.__contains__(k) -> True if D has a key k, else False"
425 return key in self.StoredVariables
428 "D.keys() -> list of D's keys"
429 return self.StoredVariables.keys()
431 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
433 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
434 sa forme mathématique la plus naturelle possible.
436 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
438 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
440 Permet de définir dans l'algorithme des paramètres requis et leurs
441 caractéristiques par défaut.
444 raise ValueError("A name is mandatory to define a required parameter.")
446 self.__required_parameters[name] = {
448 "typecast" : typecast,
454 logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
456 def getRequiredParameters(self, noDetails=True):
458 Renvoie la liste des noms de paramètres requis ou directement le
459 dictionnaire des paramètres requis.
462 ks = self.__required_parameters.keys()
466 return self.__required_parameters
468 def setParameterValue(self, name=None, value=None):
470 Renvoie la valeur d'un paramètre requis de manière contrôlée
472 default = self.__required_parameters[name]["default"]
473 typecast = self.__required_parameters[name]["typecast"]
474 minval = self.__required_parameters[name]["minval"]
475 maxval = self.__required_parameters[name]["maxval"]
476 listval = self.__required_parameters[name]["listval"]
478 if value is None and default is None:
480 elif value is None and default is not None:
481 if typecast is None: __val = default
482 else: __val = typecast( default )
484 if typecast is None: __val = value
485 else: __val = typecast( value )
487 if minval is not None and (numpy.array(__val, float) < minval).any():
488 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
489 if maxval is not None and (numpy.array(__val, float) > maxval).any():
490 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
491 if listval is not None:
492 if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
495 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
496 elif __val not in listval:
497 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
500 def __setParameters(self, fromDico={}):
502 Permet de stocker les paramètres reçus dans le dictionnaire interne.
504 self._parameters.update( fromDico )
505 for k in self.__required_parameters.keys():
506 if k in fromDico.keys():
507 self._parameters[k] = self.setParameterValue(k,fromDico[k])
509 self._parameters[k] = self.setParameterValue(k)
510 logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
512 # ==============================================================================
513 class Diagnostic(object):
515 Classe générale d'interface de type diagnostic
517 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
518 même temps que l'une des classes de persistance, comme "OneScalar" par
521 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
522 méthode "_formula" pour écrire explicitement et proprement la formule pour
523 l'écriture mathématique du calcul du diagnostic (méthode interne non
524 publique), et "calculate" pour activer la précédente tout en ayant vérifié
525 et préparé les données, et pour stocker les résultats à chaque pas (méthode
526 externe d'activation).
528 def __init__(self, name = "", parameters = {}):
530 self.name = str(name)
531 self.parameters = dict( parameters )
533 def _formula(self, *args):
535 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
536 mathématique la plus naturelle possible.
538 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
540 def calculate(self, *args):
542 Active la formule de calcul avec les arguments correctement rangés
544 raise NotImplementedError("Diagnostic activation method has not been implemented!")
546 # ==============================================================================
547 class Covariance(object):
549 Classe générale d'interface de type covariance
552 name = "GenericCovariance",
554 asEyeByScalar = None,
555 asEyeByVector = None,
560 Permet de définir une covariance :
561 - asCovariance : entrée des données, comme une matrice compatible avec
562 le constructeur de numpy.matrix
563 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
564 multiplicatif d'une matrice de corrélation identité, aucune matrice
565 n'étant donc explicitement à donner
566 - asEyeByVector : entrée des données comme un seul vecteur de variance,
567 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
568 n'étant donc explicitement à donner
569 - asCovObject : entrée des données comme un objet python, qui a les
570 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
571 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
572 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
573 - toBeChecked : booléen indiquant si le caractère SDP de la matrice
574 pleine doit être vérifié
576 self.__name = str(name)
577 self.__check = bool(toBeChecked)
580 self.__is_scalar = False
581 self.__is_vector = False
582 self.__is_matrix = False
583 self.__is_object = False
584 if asEyeByScalar is not None:
585 self.__is_scalar = True
586 self.__C = numpy.abs( float(asEyeByScalar) )
589 elif asEyeByVector is not None:
590 self.__is_vector = True
591 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
592 self.shape = (self.__C.size,self.__C.size)
593 self.size = self.__C.size**2
594 elif asCovariance is not None:
595 self.__is_matrix = True
596 self.__C = numpy.matrix( asCovariance, float )
597 self.shape = self.__C.shape
598 self.size = self.__C.size
599 elif asCovObject is not None:
600 self.__is_object = True
601 self.__C = asCovObject
602 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
603 if not hasattr(self.__C,at):
604 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
605 if hasattr(self.__C,"shape"):
606 self.shape = self.__C.shape
609 if hasattr(self.__C,"size"):
610 self.size = self.__C.size
615 # 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)
619 def __validate(self):
621 if self.ismatrix() and min(self.shape) != max(self.shape):
622 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))
623 if self.isobject() and min(self.shape) != max(self.shape):
624 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))
625 if self.isscalar() and self.__C <= 0:
626 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
627 if self.isvector() and (self.__C <= 0).any():
628 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
629 if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
631 L = numpy.linalg.cholesky( self.__C )
633 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
634 if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
636 L = self.__C.cholesky()
638 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
641 "Vérification du type interne"
642 return self.__is_scalar
645 "Vérification du type interne"
646 return self.__is_vector
649 "Vérification du type interne"
650 return self.__is_matrix
653 "Vérification du type interne"
654 return self.__is_object
659 return Covariance(self.__name+"I", asCovariance = self.__C.I )
660 elif self.isvector():
661 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
662 elif self.isscalar():
663 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
664 elif self.isobject():
665 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
667 return None # Indispensable
672 return Covariance(self.__name+"T", asCovariance = self.__C.T )
673 elif self.isvector():
674 return Covariance(self.__name+"T", asEyeByVector = self.__C )
675 elif self.isscalar():
676 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
677 elif self.isobject():
678 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
681 "Décomposition de Cholesky"
683 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
684 elif self.isvector():
685 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
686 elif self.isscalar():
687 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
688 elif self.isobject() and hasattr(self.__C,"cholesky"):
689 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
692 "Inversion de la décomposition de Cholesky"
694 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
695 elif self.isvector():
696 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
697 elif self.isscalar():
698 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
699 elif self.isobject() and hasattr(self.__C,"choleskyI"):
700 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
702 def diag(self, msize=None):
703 "Diagonale de la matrice"
705 return numpy.diag(self.__C)
706 elif self.isvector():
708 elif self.isscalar():
710 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,))
712 return self.__C * numpy.ones(int(msize))
713 elif self.isobject():
714 return self.__C.diag()
716 def asfullmatrix(self, msize=None):
720 elif self.isvector():
721 return numpy.matrix( numpy.diag(self.__C), float )
722 elif self.isscalar():
724 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,))
726 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
727 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
728 return self.__C.asfullmatrix()
730 def trace(self, msize=None):
731 "Trace de la matrice"
733 return numpy.trace(self.__C)
734 elif self.isvector():
735 return float(numpy.sum(self.__C))
736 elif self.isscalar():
738 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,))
740 return self.__C * int(msize)
741 elif self.isobject():
742 return self.__C.trace()
745 "x.__repr__() <==> repr(x)"
746 return repr(self.__C)
749 "x.__str__() <==> str(x)"
752 def __add__(self, other):
753 "x.__add__(y) <==> x+y"
754 if self.ismatrix() or self.isobject():
755 return self.__C + numpy.asmatrix(other)
756 elif self.isvector() or self.isscalar():
757 _A = numpy.asarray(other)
758 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
759 return numpy.asmatrix(_A)
761 def __radd__(self, other):
762 "x.__radd__(y) <==> y+x"
763 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
765 def __sub__(self, other):
766 "x.__sub__(y) <==> x-y"
767 if self.ismatrix() or self.isobject():
768 return self.__C - numpy.asmatrix(other)
769 elif self.isvector() or self.isscalar():
770 _A = numpy.asarray(other)
771 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
772 return numpy.asmatrix(_A)
774 def __rsub__(self, other):
775 "x.__rsub__(y) <==> y-x"
776 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
779 "x.__neg__() <==> -x"
782 def __mul__(self, other):
783 "x.__mul__(y) <==> x*y"
784 if self.ismatrix() and isinstance(other,numpy.matrix):
785 return self.__C * other
786 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
787 or isinstance(other,list) \
788 or isinstance(other,tuple)):
789 if numpy.ravel(other).size == self.shape[1]: # Vecteur
790 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
791 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
792 return self.__C * numpy.asmatrix(other)
794 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
795 elif self.isvector() and (isinstance(other,numpy.matrix) \
796 or isinstance(other,numpy.ndarray) \
797 or isinstance(other,list) \
798 or isinstance(other,tuple)):
799 if numpy.ravel(other).size == self.shape[1]: # Vecteur
800 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
801 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
802 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
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 self.__C * other
807 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
808 or isinstance(other,list) \
809 or isinstance(other,tuple)):
810 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
811 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
813 return self.__C * numpy.asmatrix(other)
814 elif self.isobject():
815 return self.__C.__mul__(other)
817 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
819 def __rmul__(self, other):
820 "x.__rmul__(y) <==> y*x"
821 if self.ismatrix() and isinstance(other,numpy.matrix):
822 return other * self.__C
823 elif self.isvector() and isinstance(other,numpy.matrix):
824 if numpy.ravel(other).size == self.shape[0]: # Vecteur
825 return numpy.asmatrix(numpy.ravel(other) * self.__C)
826 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
827 return numpy.asmatrix(numpy.array(other) * self.__C)
829 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
830 elif self.isscalar() and isinstance(other,numpy.matrix):
831 return other * self.__C
832 elif self.isobject():
833 return self.__C.__rmul__(other)
835 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
838 "x.__len__() <==> len(x)"
841 # ==============================================================================
842 def CostFunction3D(_x,
843 _Hm = None, # Pour simuler Hm(x) : HO["Direct"].appliedTo
844 _HmX = None, # Simulation déjà faite de Hm(x)
845 _arg = None, # Arguments supplementaires pour Hm, sous la forme d'un tuple
850 _SIV = False, # A résorber pour la 8.0
851 _SSC = [], # self._parameters["StoreSupplementaryCalculations"]
852 _nPS = 0, # nbPreviousSteps
853 _QM = "DA", # QualityMeasure
854 _SSV = {}, # Entrée et/ou sortie : self.StoredVariables
855 _fRt = False, # Restitue ou pas la sortie étendue
856 _sSc = True, # Stocke ou pas les SSC
859 Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
860 et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
861 DFO, QuantileRegression
867 for k in ["CostFunctionJ",
873 "SimulatedObservationAtCurrentOptimum",
874 "SimulatedObservationAtCurrentState",
878 if hasattr(_SSV[k],"store"):
879 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
881 _X = numpy.asmatrix(numpy.ravel( _x )).T
882 if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
883 _SSV["CurrentState"].append( _X )
889 raise ValueError("COSTFUNCTION3D Operator has to be defined.")
893 _HX = _Hm( _X, *_arg )
894 _HX = numpy.asmatrix(numpy.ravel( _HX )).T
896 if "SimulatedObservationAtCurrentState" in _SSC or \
897 "SimulatedObservationAtCurrentOptimum" in _SSC:
898 _SSV["SimulatedObservationAtCurrentState"].append( _HX )
900 if numpy.any(numpy.isnan(_HX)):
901 Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
903 _Y = numpy.asmatrix(numpy.ravel( _Y )).T
904 if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
905 if _BI is None or _RI is None:
906 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
907 _Xb = numpy.asmatrix(numpy.ravel( _Xb )).T
908 Jb = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
909 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
910 elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
912 raise ValueError("Observation error covariance matrix has to be properly defined!")
914 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
915 elif _QM in ["LeastSquares", "LS", "L2"]:
917 Jo = 0.5 * (_Y - _HX).T * (_Y - _HX)
918 elif _QM in ["AbsoluteValue", "L1"]:
920 Jo = numpy.sum( numpy.abs(_Y - _HX) )
921 elif _QM in ["MaximumError", "ME"]:
923 Jo = numpy.max( numpy.abs(_Y - _HX) )
924 elif _QM in ["QR", "Null"]:
928 raise ValueError("Unknown asked quality measure!")
930 J = float( Jb ) + float( Jo )
933 _SSV["CostFunctionJb"].append( Jb )
934 _SSV["CostFunctionJo"].append( Jo )
935 _SSV["CostFunctionJ" ].append( J )
937 if "IndexOfOptimum" in _SSC or \
938 "CurrentOptimum" in _SSC or \
939 "SimulatedObservationAtCurrentOptimum" in _SSC:
940 IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
941 if "IndexOfOptimum" in _SSC:
942 _SSV["IndexOfOptimum"].append( IndexMin )
943 if "CurrentOptimum" in _SSC:
944 _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
945 if "SimulatedObservationAtCurrentOptimum" in _SSC:
946 _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
951 if _QM in ["QR"]: # Pour le QuantileRegression
956 # ==============================================================================
957 if __name__ == "__main__":
958 print '\n AUTODIAGNOSTIC \n'