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 %.0f 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 %.0f 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 if numpy.matrix(asEyeByScalar).size != 1:
586 raise ValueError(' The diagonal multiplier given to define a sparse matrix is not a unique scalar value.\n Its actual measured size is %i. Please check your scalar input.'%numpy.matrix(asEyeByScalar).size)
587 self.__is_scalar = True
588 self.__C = numpy.abs( float(asEyeByScalar) )
591 elif asEyeByVector is not None:
592 self.__is_vector = True
593 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
594 self.shape = (self.__C.size,self.__C.size)
595 self.size = self.__C.size**2
596 elif asCovariance is not None:
597 self.__is_matrix = True
598 self.__C = numpy.matrix( asCovariance, float )
599 self.shape = self.__C.shape
600 self.size = self.__C.size
601 elif asCovObject is not None:
602 self.__is_object = True
603 self.__C = asCovObject
604 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
605 if not hasattr(self.__C,at):
606 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
607 if hasattr(self.__C,"shape"):
608 self.shape = self.__C.shape
611 if hasattr(self.__C,"size"):
612 self.size = self.__C.size
617 # 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)
621 def __validate(self):
623 if self.ismatrix() and min(self.shape) != max(self.shape):
624 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))
625 if self.isobject() and min(self.shape) != max(self.shape):
626 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))
627 if self.isscalar() and self.__C <= 0:
628 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
629 if self.isvector() and (self.__C <= 0).any():
630 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
631 if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
633 L = numpy.linalg.cholesky( self.__C )
635 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
636 if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
638 L = self.__C.cholesky()
640 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
643 "Vérification du type interne"
644 return self.__is_scalar
647 "Vérification du type interne"
648 return self.__is_vector
651 "Vérification du type interne"
652 return self.__is_matrix
655 "Vérification du type interne"
656 return self.__is_object
661 return Covariance(self.__name+"I", asCovariance = self.__C.I )
662 elif self.isvector():
663 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
664 elif self.isscalar():
665 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
666 elif self.isobject():
667 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
669 return None # Indispensable
674 return Covariance(self.__name+"T", asCovariance = self.__C.T )
675 elif self.isvector():
676 return Covariance(self.__name+"T", asEyeByVector = self.__C )
677 elif self.isscalar():
678 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
679 elif self.isobject():
680 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
683 "Décomposition de Cholesky"
685 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
686 elif self.isvector():
687 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
688 elif self.isscalar():
689 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
690 elif self.isobject() and hasattr(self.__C,"cholesky"):
691 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
694 "Inversion de la décomposition de Cholesky"
696 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
697 elif self.isvector():
698 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
699 elif self.isscalar():
700 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
701 elif self.isobject() and hasattr(self.__C,"choleskyI"):
702 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
704 def diag(self, msize=None):
705 "Diagonale de la matrice"
707 return numpy.diag(self.__C)
708 elif self.isvector():
710 elif self.isscalar():
712 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,))
714 return self.__C * numpy.ones(int(msize))
715 elif self.isobject():
716 return self.__C.diag()
718 def asfullmatrix(self, msize=None):
722 elif self.isvector():
723 return numpy.matrix( numpy.diag(self.__C), float )
724 elif self.isscalar():
726 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,))
728 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
729 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
730 return self.__C.asfullmatrix()
732 def trace(self, msize=None):
733 "Trace de la matrice"
735 return numpy.trace(self.__C)
736 elif self.isvector():
737 return float(numpy.sum(self.__C))
738 elif self.isscalar():
740 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,))
742 return self.__C * int(msize)
743 elif self.isobject():
744 return self.__C.trace()
747 "x.__repr__() <==> repr(x)"
748 return repr(self.__C)
751 "x.__str__() <==> str(x)"
754 def __add__(self, other):
755 "x.__add__(y) <==> x+y"
756 if self.ismatrix() or self.isobject():
757 return self.__C + numpy.asmatrix(other)
758 elif self.isvector() or self.isscalar():
759 _A = numpy.asarray(other)
760 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
761 return numpy.asmatrix(_A)
763 def __radd__(self, other):
764 "x.__radd__(y) <==> y+x"
765 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
767 def __sub__(self, other):
768 "x.__sub__(y) <==> x-y"
769 if self.ismatrix() or self.isobject():
770 return self.__C - numpy.asmatrix(other)
771 elif self.isvector() or self.isscalar():
772 _A = numpy.asarray(other)
773 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
774 return numpy.asmatrix(_A)
776 def __rsub__(self, other):
777 "x.__rsub__(y) <==> y-x"
778 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
781 "x.__neg__() <==> -x"
784 def __mul__(self, other):
785 "x.__mul__(y) <==> x*y"
786 if self.ismatrix() and isinstance(other,numpy.matrix):
787 return self.__C * other
788 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
789 or isinstance(other,list) \
790 or isinstance(other,tuple)):
791 if numpy.ravel(other).size == self.shape[1]: # Vecteur
792 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
793 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
794 return self.__C * numpy.asmatrix(other)
796 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
797 elif self.isvector() and (isinstance(other,numpy.matrix) \
798 or isinstance(other,numpy.ndarray) \
799 or isinstance(other,list) \
800 or isinstance(other,tuple)):
801 if numpy.ravel(other).size == self.shape[1]: # Vecteur
802 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
803 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
804 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
806 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
807 elif self.isscalar() and isinstance(other,numpy.matrix):
808 return self.__C * other
809 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
810 or isinstance(other,list) \
811 or isinstance(other,tuple)):
812 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
813 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
815 return self.__C * numpy.asmatrix(other)
816 elif self.isobject():
817 return self.__C.__mul__(other)
819 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
821 def __rmul__(self, other):
822 "x.__rmul__(y) <==> y*x"
823 if self.ismatrix() and isinstance(other,numpy.matrix):
824 return other * self.__C
825 elif self.isvector() and isinstance(other,numpy.matrix):
826 if numpy.ravel(other).size == self.shape[0]: # Vecteur
827 return numpy.asmatrix(numpy.ravel(other) * self.__C)
828 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
829 return numpy.asmatrix(numpy.array(other) * self.__C)
831 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
832 elif self.isscalar() and isinstance(other,numpy.matrix):
833 return other * self.__C
834 elif self.isobject():
835 return self.__C.__rmul__(other)
837 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
840 "x.__len__() <==> len(x)"
843 # ==============================================================================
844 def CostFunction3D(_x,
845 _Hm = None, # Pour simuler Hm(x) : HO["Direct"].appliedTo
846 _HmX = None, # Simulation déjà faite de Hm(x)
847 _arg = None, # Arguments supplementaires pour Hm, sous la forme d'un tuple
852 _SIV = False, # A résorber pour la 8.0
853 _SSC = [], # self._parameters["StoreSupplementaryCalculations"]
854 _nPS = 0, # nbPreviousSteps
855 _QM = "DA", # QualityMeasure
856 _SSV = {}, # Entrée et/ou sortie : self.StoredVariables
857 _fRt = False, # Restitue ou pas la sortie étendue
858 _sSc = True, # Stocke ou pas les SSC
861 Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
862 et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
863 DFO, QuantileRegression
869 for k in ["CostFunctionJ",
875 "SimulatedObservationAtCurrentOptimum",
876 "SimulatedObservationAtCurrentState",
880 if hasattr(_SSV[k],"store"):
881 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
883 _X = numpy.asmatrix(numpy.ravel( _x )).T
884 if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
885 _SSV["CurrentState"].append( _X )
891 raise ValueError("COSTFUNCTION3D Operator has to be defined.")
895 _HX = _Hm( _X, *_arg )
896 _HX = numpy.asmatrix(numpy.ravel( _HX )).T
898 if "SimulatedObservationAtCurrentState" in _SSC or \
899 "SimulatedObservationAtCurrentOptimum" in _SSC:
900 _SSV["SimulatedObservationAtCurrentState"].append( _HX )
902 if numpy.any(numpy.isnan(_HX)):
903 Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
905 _Y = numpy.asmatrix(numpy.ravel( _Y )).T
906 if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
907 if _BI is None or _RI is None:
908 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
909 _Xb = numpy.asmatrix(numpy.ravel( _Xb )).T
910 Jb = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
911 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
912 elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
914 raise ValueError("Observation error covariance matrix has to be properly defined!")
916 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
917 elif _QM in ["LeastSquares", "LS", "L2"]:
919 Jo = 0.5 * (_Y - _HX).T * (_Y - _HX)
920 elif _QM in ["AbsoluteValue", "L1"]:
922 Jo = numpy.sum( numpy.abs(_Y - _HX) )
923 elif _QM in ["MaximumError", "ME"]:
925 Jo = numpy.max( numpy.abs(_Y - _HX) )
926 elif _QM in ["QR", "Null"]:
930 raise ValueError("Unknown asked quality measure!")
932 J = float( Jb ) + float( Jo )
935 _SSV["CostFunctionJb"].append( Jb )
936 _SSV["CostFunctionJo"].append( Jo )
937 _SSV["CostFunctionJ" ].append( J )
939 if "IndexOfOptimum" in _SSC or \
940 "CurrentOptimum" in _SSC or \
941 "SimulatedObservationAtCurrentOptimum" in _SSC:
942 IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
943 if "IndexOfOptimum" in _SSC:
944 _SSV["IndexOfOptimum"].append( IndexMin )
945 if "CurrentOptimum" in _SSC:
946 _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
947 if "SimulatedObservationAtCurrentOptimum" in _SSC:
948 _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
953 if _QM in ["QR"]: # Pour le QuantileRegression
958 # ==============================================================================
959 if __name__ == "__main__":
960 print('\n AUTODIAGNOSTIC \n')