1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2015 EDF R&D
5 # This library is free software; you can redistribute it and/or
6 # modify it under the terms of the GNU Lesser General Public
7 # License as published by the Free Software Foundation; either
8 # version 2.1 of the License.
10 # This library is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 # Lesser General Public License for more details.
15 # You should have received a copy of the GNU Lesser General Public
16 # License along with this library; if not, write to the Free Software
17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 Définit les outils généraux élémentaires.
26 Ce module est destiné à etre appelée par AssimilationStudy pour constituer
27 les objets élémentaires de l'algorithme.
29 __author__ = "Jean-Philippe ARGAUD"
37 # ==============================================================================
40 Classe générale de gestion d'un cache de calculs
43 toleranceInRedundancy = 1.e-18,
44 lenghtOfRedundancy = -1,
47 Les caractéristiques de tolérance peuvent être modifées à la création.
49 self.__tolerBP = float(toleranceInRedundancy)
50 self.__lenghtOR = int(lenghtOfRedundancy)
51 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="" ):
61 for i in xrange(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
62 if xValue.size != self.__listOPCV[i][0].size:
63 # 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))
65 if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
67 __HxV = self.__listOPCV[i][1]
68 # logging.debug("CM Cas%s déja calculé, portant le numéro %i"%(info,i))
72 def storeValueInX(self, xValue, HxValue ):
73 if self.__lenghtOR < 0:
74 self.__lenghtOR = 2 * xValue.size + 2
75 self.__initlnOR = self.__lenghtOR
76 while len(self.__listOPCV) > self.__lenghtOR:
77 # logging.debug("CM Réduction de la liste des cas à %i éléments par suppression du premier"%self.__lenghtOR)
78 self.__listOPCV.pop(0)
79 self.__listOPCV.append( (
80 copy.copy(numpy.ravel(xValue)),
82 numpy.linalg.norm(xValue),
86 self.__initlnOR = self.__lenghtOR
90 self.__lenghtOR = self.__initlnOR
92 # ==============================================================================
95 Classe générale d'interface de type opérateur
102 def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
104 On construit un objet de ce type en fournissant à l'aide de l'un des
105 deux mots-clé, soit une fonction python, soit une matrice.
107 - fromMethod : argument de type fonction Python
108 - fromMatrix : argument adapté au constructeur numpy.matrix
109 - avoidingRedundancy : évite ou pas les calculs redondants
111 self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
112 self.__AvoidRC = bool( avoidingRedundancy )
113 if fromMethod is not None:
114 self.__Method = fromMethod
116 self.__Type = "Method"
117 elif fromMatrix is not None:
119 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
120 self.__Type = "Matrix"
126 def disableAvoidingRedundancy(self):
127 Operator.CM.disable()
129 def enableAvoidingRedundancy(self):
133 Operator.CM.disable()
138 def appliedTo(self, xValue):
140 Permet de restituer le résultat de l'application de l'opérateur à un
141 argument xValue. Cette méthode se contente d'appliquer, son argument
142 devant a priori être du bon type.
144 - xValue : argument adapté pour appliquer l'opérateur
147 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
149 __alreadyCalculated = False
151 if __alreadyCalculated:
152 self.__addOneCacheCall()
155 if self.__Matrix is not None:
156 self.__addOneMatrixCall()
157 HxValue = self.__Matrix * xValue
159 self.__addOneMethodCall()
160 HxValue = self.__Method( xValue )
162 Operator.CM.storeValueInX(xValue,HxValue)
166 def appliedControledFormTo(self, (xValue, uValue) ):
168 Permet de restituer le résultat de l'application de l'opérateur à une
169 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
170 argument devant a priori être du bon type. Si la uValue est None,
171 on suppose que l'opérateur ne s'applique qu'à xValue.
173 - xValue : argument X adapté pour appliquer l'opérateur
174 - uValue : argument U adapté pour appliquer l'opérateur
176 if self.__Matrix is not None:
177 self.__addOneMatrixCall()
178 return self.__Matrix * xValue
179 elif uValue is not None:
180 self.__addOneMethodCall()
181 return self.__Method( (xValue, uValue) )
183 self.__addOneMethodCall()
184 return self.__Method( xValue )
186 def appliedInXTo(self, (xNominal, xValue) ):
188 Permet de restituer le résultat de l'application de l'opérateur à un
189 argument xValue, sachant que l'opérateur est valable en xNominal.
190 Cette méthode se contente d'appliquer, son argument devant a priori
191 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
192 alors il est valable en tout point nominal et il n'est pas nécessaire
194 Arguments : une liste contenant
195 - xNominal : argument permettant de donner le point où l'opérateur
196 est construit pour etre ensuite appliqué
197 - xValue : argument adapté pour appliquer l'opérateur
199 if self.__Matrix is not None:
200 self.__addOneMatrixCall()
201 return self.__Matrix * xValue
203 self.__addOneMethodCall()
204 return self.__Method( (xNominal, xValue) )
206 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
208 Permet de renvoyer l'opérateur sous la forme d'une matrice
210 if self.__Matrix is not None:
211 self.__addOneMatrixCall()
213 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
214 self.__addOneMethodCall()
215 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
217 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
221 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
222 la forme d'une matrice
224 if self.__Matrix is not None:
225 return self.__Matrix.shape
227 raise ValueError("Matrix form of the operator is not available, nor the shape")
229 def nbcalls(self, which=None):
231 Renvoie les nombres d'évaluations de l'opérateur
234 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
235 self.__NbCallsAsMatrix,
236 self.__NbCallsAsMethod,
237 self.__NbCallsOfCached,
238 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
239 Operator.NbCallsAsMatrix,
240 Operator.NbCallsAsMethod,
241 Operator.NbCallsOfCached,
243 if which is None: return __nbcalls
244 else: return __nbcalls[which]
246 def __addOneMatrixCall(self):
247 self.__NbCallsAsMatrix += 1 # Decompte local
248 Operator.NbCallsAsMatrix += 1 # Decompte global
250 def __addOneMethodCall(self):
251 self.__NbCallsAsMethod += 1 # Decompte local
252 Operator.NbCallsAsMethod += 1 # Decompte global
254 def __addOneCacheCall(self):
255 self.__NbCallsOfCached += 1 # Decompte local
256 Operator.NbCallsOfCached += 1 # Decompte global
258 # ==============================================================================
261 Classe générale d'interface de type algorithme
263 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
264 d'assimilation, en fournissant un container (dictionnaire) de variables
265 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
267 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
269 def __init__(self, name):
271 L'initialisation présente permet de fabriquer des variables de stockage
272 disponibles de manière générique dans les algorithmes élémentaires. Ces
273 variables de stockage sont ensuite conservées dans un dictionnaire
274 interne à l'objet, mais auquel on accède par la méthode "get".
276 Les variables prévues sont :
277 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
278 - CostFunctionJb : partie ébauche ou background de la fonction-cout
279 - CostFunctionJo : partie observations de la fonction-cout
280 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
281 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
282 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
283 - CurrentState : état courant lors d'itérations
284 - Analysis : l'analyse Xa
285 - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
286 - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
287 - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
288 - Innovation : l'innovation : d = Y - H(X)
289 - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
290 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
291 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
292 - MahalanobisConsistency : indicateur de consistance des covariances
293 - OMA : Observation moins Analysis : Y - Xa
294 - OMB : Observation moins Background : Y - Xb
295 - AMB : Analysis moins Background : Xa - Xb
296 - APosterioriCovariance : matrice A
297 - APosterioriVariances : variances de la matrice A
298 - APosterioriStandardDeviations : écart-types de la matrice A
299 - APosterioriCorrelations : correlations de la matrice A
300 On peut rajouter des variables à stocker dans l'initialisation de
301 l'algorithme élémentaire qui va hériter de cette classe
303 logging.debug("%s Initialisation"%str(name))
304 self._m = PlatformInfo.SystemUsage()
306 self._name = str( name )
307 self._parameters = {"StoreSupplementaryCalculations":[]}
308 self.__required_parameters = {}
309 self.StoredVariables = {}
311 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
312 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
313 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
314 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
315 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
316 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
317 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
318 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
319 self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum")
320 self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
321 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
322 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
323 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
324 self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
325 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
326 self.StoredVariables["InnovationAtCurrentState"] = Persistence.OneVector(name = "InnovationAtCurrentState")
327 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
328 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
329 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
330 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
331 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
332 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
333 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
334 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
335 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
336 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
337 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
340 logging.debug("%s Lancement"%self._name)
341 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
344 def _post_run(self,_oH=None):
345 if ("StoreSupplementaryCalculations" in self._parameters) and \
346 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
347 for _A in self.StoredVariables["APosterioriCovariance"]:
348 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
349 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
350 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
351 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
352 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
353 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
354 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
355 self.StoredVariables["APosterioriCorrelations"].store( _C )
357 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)))
358 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)))
359 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
360 logging.debug("%s Terminé"%self._name)
363 def get(self, key=None):
365 Renvoie l'une des variables stockées identifiée par la clé, ou le
366 dictionnaire de l'ensemble des variables disponibles en l'absence de
367 clé. Ce sont directement les variables sous forme objet qui sont
368 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
369 des classes de persistance.
372 return self.StoredVariables[key]
374 return self.StoredVariables
376 def __contains__(self, key=None):
378 Vérifie si l'une des variables stockées est identifiée par la clé.
380 return (key in self.StoredVariables)
384 Renvoie la liste des clés de variables stockées.
386 return self.StoredVariables.keys()
388 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
390 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
391 sa forme mathématique la plus naturelle possible.
393 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
395 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
397 Permet de définir dans l'algorithme des paramètres requis et leurs
398 caractéristiques par défaut.
401 raise ValueError("A name is mandatory to define a required parameter.")
403 self.__required_parameters[name] = {
405 "typecast" : typecast,
411 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
413 def getRequiredParameters(self, noDetails=True):
415 Renvoie la liste des noms de paramètres requis ou directement le
416 dictionnaire des paramètres requis.
419 ks = self.__required_parameters.keys()
423 return self.__required_parameters
425 def setParameterValue(self, name=None, value=None):
427 Renvoie la valeur d'un paramètre requis de manière contrôlée
429 default = self.__required_parameters[name]["default"]
430 typecast = self.__required_parameters[name]["typecast"]
431 minval = self.__required_parameters[name]["minval"]
432 maxval = self.__required_parameters[name]["maxval"]
433 listval = self.__required_parameters[name]["listval"]
435 if value is None and default is None:
437 elif value is None and default is not None:
438 if typecast is None: __val = default
439 else: __val = typecast( default )
441 if typecast is None: __val = value
442 else: __val = typecast( value )
444 if minval is not None and (numpy.array(__val, float) < minval).any():
445 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
446 if maxval is not None and (numpy.array(__val, float) > maxval).any():
447 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
448 if listval is not None:
449 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
452 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
453 elif __val not in listval:
454 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
457 def setParameters(self, fromDico={}):
459 Permet de stocker les paramètres reçus dans le dictionnaire interne.
461 self._parameters.update( fromDico )
462 for k in self.__required_parameters.keys():
463 if k in fromDico.keys():
464 self._parameters[k] = self.setParameterValue(k,fromDico[k])
466 self._parameters[k] = self.setParameterValue(k)
467 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
469 # ==============================================================================
472 Classe générale d'interface de type diagnostic
474 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
475 même temps que l'une des classes de persistance, comme "OneScalar" par
478 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
479 méthode "_formula" pour écrire explicitement et proprement la formule pour
480 l'écriture mathématique du calcul du diagnostic (méthode interne non
481 publique), et "calculate" pour activer la précédente tout en ayant vérifié
482 et préparé les données, et pour stocker les résultats à chaque pas (méthode
483 externe d'activation).
485 def __init__(self, name = "", parameters = {}):
486 self.name = str(name)
487 self.parameters = dict( parameters )
489 def _formula(self, *args):
491 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
492 mathématique la plus naturelle possible.
494 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
496 def calculate(self, *args):
498 Active la formule de calcul avec les arguments correctement rangés
500 raise NotImplementedError("Diagnostic activation method has not been implemented!")
502 # ==============================================================================
505 Classe générale d'interface de type covariance
508 name = "GenericCovariance",
510 asEyeByScalar = None,
511 asEyeByVector = None,
515 Permet de définir une covariance :
516 - asCovariance : entrée des données, comme une matrice compatible avec
517 le constructeur de numpy.matrix
518 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
519 multiplicatif d'une matrice de corrélation identité, aucune matrice
520 n'étant donc explicitement à donner
521 - asEyeByVector : entrée des données comme un seul vecteur de variance,
522 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
523 n'étant donc explicitement à donner
524 - asCovObject : entrée des données comme un objet python, qui a les
525 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
526 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
527 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
529 self.__name = str(name)
532 self.__is_scalar = False
533 self.__is_vector = False
534 self.__is_matrix = False
535 self.__is_object = False
536 if asEyeByScalar is not None:
537 self.__is_scalar = True
538 self.__C = numpy.abs( float(asEyeByScalar) )
541 elif asEyeByVector is not None:
542 self.__is_vector = True
543 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
544 self.shape = (self.__C.size,self.__C.size)
545 self.size = self.__C.size**2
546 elif asCovariance is not None:
547 self.__is_matrix = True
548 self.__C = numpy.matrix( asCovariance, float )
549 self.shape = self.__C.shape
550 self.size = self.__C.size
551 elif asCovObject is not None:
552 self.__is_object = True
553 self.__C = asCovObject
554 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
555 if not hasattr(self.__C,at):
556 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
557 if hasattr(self.__C,"shape"):
558 self.shape = self.__C.shape
561 if hasattr(self.__C,"size"):
562 self.size = self.__C.size
567 # 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)
571 def __validate(self):
572 if self.ismatrix() and min(self.shape) != max(self.shape):
573 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))
574 if self.isobject() and min(self.shape) != max(self.shape):
575 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))
576 if self.isscalar() and self.__C <= 0:
577 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
578 if self.isvector() and (self.__C <= 0).any():
579 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
580 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
582 L = numpy.linalg.cholesky( self.__C )
584 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
587 return self.__is_scalar
590 return self.__is_vector
593 return self.__is_matrix
596 return self.__is_object
600 return Covariance(self.__name+"I", asCovariance = self.__C.I )
601 elif self.isvector():
602 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
603 elif self.isscalar():
604 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
605 elif self.isobject():
606 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
608 return None # Indispensable
612 return Covariance(self.__name+"T", asCovariance = self.__C.T )
613 elif self.isvector():
614 return Covariance(self.__name+"T", asEyeByVector = self.__C )
615 elif self.isscalar():
616 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
617 elif self.isobject():
618 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
622 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
623 elif self.isvector():
624 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
625 elif self.isscalar():
626 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
627 elif self.isobject() and hasattr(self.__C,"cholesky"):
628 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
632 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
633 elif self.isvector():
634 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
635 elif self.isscalar():
636 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
637 elif self.isobject() and hasattr(self.__C,"choleskyI"):
638 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
640 def diag(self, msize=None):
642 return numpy.diag(self.__C)
643 elif self.isvector():
645 elif self.isscalar():
647 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,))
649 return self.__C * numpy.ones(int(msize))
650 elif self.isobject():
651 return self.__C.diag()
653 def asfullmatrix(self, msize=None):
656 elif self.isvector():
657 return numpy.matrix( numpy.diag(self.__C), float )
658 elif self.isscalar():
660 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,))
662 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
663 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
664 return self.__C.asfullmatrix()
666 def trace(self, msize=None):
668 return numpy.trace(self.__C)
669 elif self.isvector():
670 return float(numpy.sum(self.__C))
671 elif self.isscalar():
673 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,))
675 return self.__C * int(msize)
676 elif self.isobject():
677 return self.__C.trace()
680 return repr(self.__C)
685 def __add__(self, other):
686 if self.ismatrix() or self.isobject():
687 return self.__C + numpy.asmatrix(other)
688 elif self.isvector() or self.isscalar():
689 _A = numpy.asarray(other)
690 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
691 return numpy.asmatrix(_A)
693 def __radd__(self, other):
694 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
696 def __sub__(self, other):
697 if self.ismatrix() or self.isobject():
698 return self.__C - numpy.asmatrix(other)
699 elif self.isvector() or self.isscalar():
700 _A = numpy.asarray(other)
701 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
702 return numpy.asmatrix(_A)
704 def __rsub__(self, other):
705 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
710 def __mul__(self, other):
711 if self.ismatrix() and isinstance(other,numpy.matrix):
712 return self.__C * other
713 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
714 or isinstance(other,list) \
715 or isinstance(other,tuple)):
716 if numpy.ravel(other).size == self.shape[1]: # Vecteur
717 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
718 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
719 return self.__C * numpy.asmatrix(other)
721 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
722 elif self.isvector() and (isinstance(other,numpy.matrix) \
723 or isinstance(other,numpy.ndarray) \
724 or isinstance(other,list) \
725 or isinstance(other,tuple)):
726 if numpy.ravel(other).size == self.shape[1]: # Vecteur
727 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
728 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
729 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
731 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
732 elif self.isscalar() and isinstance(other,numpy.matrix):
733 return self.__C * other
734 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
735 or isinstance(other,list) \
736 or isinstance(other,tuple)):
737 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
738 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
740 return self.__C * numpy.asmatrix(other)
741 elif self.isobject():
742 return self.__C.__mul__(other)
744 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
746 def __rmul__(self, other):
747 if self.ismatrix() and isinstance(other,numpy.matrix):
748 return other * self.__C
749 elif self.isvector() and isinstance(other,numpy.matrix):
750 if numpy.ravel(other).size == self.shape[0]: # Vecteur
751 return numpy.asmatrix(numpy.ravel(other) * self.__C)
752 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
753 return numpy.asmatrix(numpy.array(other) * self.__C)
755 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
756 elif self.isscalar() and isinstance(other,numpy.matrix):
757 return other * self.__C
758 elif self.isobject():
759 return self.__C.__rmul__(other)
761 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
766 # ==============================================================================
769 _Hm = None, # Pour simuler Hm(x) : HO["Direct"].appliedTo
770 _HmX = None, # Simulation déjà faite de Hm(x)
771 _arg = None, # Arguments supplementaires pour Hm, sous la forme d'un tuple
776 _SIV = False, # A résorber pour la 8.0
777 _SSC = [], # self._parameters["StoreSupplementaryCalculations"]
778 _nPS = 0, # nbPreviousSteps
779 _QM = "DA", # QualityMeasure
780 _SSV = {}, # Entrée et/ou sortie : self.StoredVariables
781 _fRt = False, # Restitue ou pas la sortie étendue
782 _sSc = True, # Stocke ou pas les SSC
785 Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
786 et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
787 DFO, QuantileRegression
793 for k in ["CostFunctionJ",
799 "SimulatedObservationAtCurrentOptimum",
800 "SimulatedObservationAtCurrentState",
804 if hasattr(_SSV[k],"store"):
805 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
807 _X = numpy.asmatrix(numpy.ravel( _x )).T
808 if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
809 _SSV["CurrentState"].append( _X )
815 raise ValueError("%s Operator has to be defined."%(self.__name,))
819 _HX = _Hm( _X, *_arg )
820 _HX = numpy.asmatrix(numpy.ravel( _HX )).T
822 if "SimulatedObservationAtCurrentState" in _SSC or \
823 "SimulatedObservationAtCurrentOptimum" in _SSC:
824 _SSV["SimulatedObservationAtCurrentState"].append( _HX )
826 if numpy.any(numpy.isnan(_HX)):
827 Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
829 _Y = numpy.asmatrix(numpy.ravel( _Y )).T
830 if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
831 if _BI is None or _RI is None:
832 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
833 _Xb = numpy.asmatrix(numpy.ravel( _Xb )).T
834 Jb = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
835 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
836 elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
838 raise ValueError("Observation error covariance matrix has to be properly defined!")
840 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
841 elif _QM in ["LeastSquares", "LS", "L2"]:
843 Jo = 0.5 * (_Y - _HX).T * (_Y - _HX)
844 elif _QM in ["AbsoluteValue", "L1"]:
846 Jo = numpy.sum( numpy.abs(_Y - _HX) )
847 elif _QM in ["MaximumError", "ME"]:
849 Jo = numpy.max( numpy.abs(_Y - _HX) )
850 elif _QM in ["QR", "Null"]:
854 raise ValueError("Unknown asked quality measure!")
856 J = float( Jb ) + float( Jo )
859 _SSV["CostFunctionJb"].append( Jb )
860 _SSV["CostFunctionJo"].append( Jo )
861 _SSV["CostFunctionJ" ].append( J )
863 if "IndexOfOptimum" in _SSC or \
864 "CurrentOptimum" in _SSC or \
865 "SimulatedObservationAtCurrentOptimum" in _SSC:
866 IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
867 if "IndexOfOptimum" in _SSC:
868 _SSV["IndexOfOptimum"].append( IndexMin )
869 if "CurrentOptimum" in _SSC:
870 _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
871 if "SimulatedObservationAtCurrentOptimum" in _SSC:
872 _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
877 if _QM in ["QR"]: # Pour le QuantileRegression
882 # ==============================================================================
883 if __name__ == "__main__":
884 print '\n AUTODIAGNOSTIC \n'