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 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
290 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
291 - MahalanobisConsistency : indicateur de consistance des covariances
292 - OMA : Observation moins Analysis : Y - Xa
293 - OMB : Observation moins Background : Y - Xb
294 - AMB : Analysis moins Background : Xa - Xb
295 - APosterioriCovariance : matrice A
296 - APosterioriVariances : variances de la matrice A
297 - APosterioriStandardDeviations : écart-types de la matrice A
298 - APosterioriCorrelations : correlations de la matrice A
299 On peut rajouter des variables à stocker dans l'initialisation de
300 l'algorithme élémentaire qui va hériter de cette classe
302 logging.debug("%s Initialisation"%str(name))
303 self._m = PlatformInfo.SystemUsage()
305 self._name = str( name )
306 self._parameters = {"StoreSupplementaryCalculations":[]}
307 self.__required_parameters = {}
308 self.StoredVariables = {}
310 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
311 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
312 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
313 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
314 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
315 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
316 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
317 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
318 self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum")
319 self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
320 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
321 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
322 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
323 self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
324 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
325 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
326 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
327 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
328 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
329 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
330 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
331 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
332 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
333 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
334 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
335 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
338 logging.debug("%s Lancement"%self._name)
339 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
342 def _post_run(self,_oH=None):
343 if self._parameters.has_key("StoreSupplementaryCalculations") and \
344 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
345 for _A in self.StoredVariables["APosterioriCovariance"]:
346 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
347 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
348 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
349 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
350 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
351 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
352 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
353 self.StoredVariables["APosterioriCorrelations"].store( _C )
355 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)))
356 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)))
357 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
358 logging.debug("%s Terminé"%self._name)
361 def get(self, key=None):
363 Renvoie l'une des variables stockées identifiée par la clé, ou le
364 dictionnaire de l'ensemble des variables disponibles en l'absence de
365 clé. Ce sont directement les variables sous forme objet qui sont
366 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
367 des classes de persistance.
370 return self.StoredVariables[key]
372 return self.StoredVariables
374 def has_key(self, key=None):
376 Vérifie si l'une des variables stockées est identifiée par la clé.
378 return self.StoredVariables.has_key(key)
382 Renvoie la liste des clés de variables stockées.
384 return self.StoredVariables.keys()
386 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
388 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
389 sa forme mathématique la plus naturelle possible.
391 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
393 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
395 Permet de définir dans l'algorithme des paramètres requis et leurs
396 caractéristiques par défaut.
399 raise ValueError("A name is mandatory to define a required parameter.")
401 self.__required_parameters[name] = {
403 "typecast" : typecast,
409 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
411 def getRequiredParameters(self, noDetails=True):
413 Renvoie la liste des noms de paramètres requis ou directement le
414 dictionnaire des paramètres requis.
417 ks = self.__required_parameters.keys()
421 return self.__required_parameters
423 def setParameterValue(self, name=None, value=None):
425 Renvoie la valeur d'un paramètre requis de manière contrôlée
427 default = self.__required_parameters[name]["default"]
428 typecast = self.__required_parameters[name]["typecast"]
429 minval = self.__required_parameters[name]["minval"]
430 maxval = self.__required_parameters[name]["maxval"]
431 listval = self.__required_parameters[name]["listval"]
433 if value is None and default is None:
435 elif value is None and default is not None:
436 if typecast is None: __val = default
437 else: __val = typecast( default )
439 if typecast is None: __val = value
440 else: __val = typecast( value )
442 if minval is not None and (numpy.array(__val, float) < minval).any():
443 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
444 if maxval is not None and (numpy.array(__val, float) > maxval).any():
445 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
446 if listval is not None:
447 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
450 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
451 elif __val not in listval:
452 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
455 def setParameters(self, fromDico={}):
457 Permet de stocker les paramètres reçus dans le dictionnaire interne.
459 self._parameters.update( fromDico )
460 for k in self.__required_parameters.keys():
461 if k in fromDico.keys():
462 self._parameters[k] = self.setParameterValue(k,fromDico[k])
464 self._parameters[k] = self.setParameterValue(k)
465 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
467 # ==============================================================================
470 Classe générale d'interface de type diagnostic
472 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
473 même temps que l'une des classes de persistance, comme "OneScalar" par
476 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
477 méthode "_formula" pour écrire explicitement et proprement la formule pour
478 l'écriture mathématique du calcul du diagnostic (méthode interne non
479 publique), et "calculate" pour activer la précédente tout en ayant vérifié
480 et préparé les données, et pour stocker les résultats à chaque pas (méthode
481 externe d'activation).
483 def __init__(self, name = "", parameters = {}):
484 self.name = str(name)
485 self.parameters = dict( parameters )
487 def _formula(self, *args):
489 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
490 mathématique la plus naturelle possible.
492 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
494 def calculate(self, *args):
496 Active la formule de calcul avec les arguments correctement rangés
498 raise NotImplementedError("Diagnostic activation method has not been implemented!")
500 # ==============================================================================
503 Classe générale d'interface de type covariance
506 name = "GenericCovariance",
508 asEyeByScalar = None,
509 asEyeByVector = None,
513 Permet de définir une covariance :
514 - asCovariance : entrée des données, comme une matrice compatible avec
515 le constructeur de numpy.matrix
516 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
517 multiplicatif d'une matrice de corrélation identité, aucune matrice
518 n'étant donc explicitement à donner
519 - asEyeByVector : entrée des données comme un seul vecteur de variance,
520 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
521 n'étant donc explicitement à donner
522 - asCovObject : entrée des données comme un objet python, qui a les
523 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
524 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
525 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
527 self.__name = str(name)
530 self.__is_scalar = False
531 self.__is_vector = False
532 self.__is_matrix = False
533 self.__is_object = False
534 if asEyeByScalar is not None:
535 self.__is_scalar = True
536 self.__C = numpy.abs( float(asEyeByScalar) )
539 elif asEyeByVector is not None:
540 self.__is_vector = True
541 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
542 self.shape = (self.__C.size,self.__C.size)
543 self.size = self.__C.size**2
544 elif asCovariance is not None:
545 self.__is_matrix = True
546 self.__C = numpy.matrix( asCovariance, float )
547 self.shape = self.__C.shape
548 self.size = self.__C.size
549 elif asCovObject is not None:
550 self.__is_object = True
551 self.__C = asCovObject
552 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
553 if not hasattr(self.__C,at):
554 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
555 if hasattr(self.__C,"shape"):
556 self.shape = self.__C.shape
559 if hasattr(self.__C,"size"):
560 self.size = self.__C.size
565 # 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)
569 def __validate(self):
570 if self.ismatrix() and min(self.shape) != max(self.shape):
571 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))
572 if self.isobject() and min(self.shape) != max(self.shape):
573 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))
574 if self.isscalar() and self.__C <= 0:
575 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
576 if self.isvector() and (self.__C <= 0).any():
577 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
578 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
580 L = numpy.linalg.cholesky( self.__C )
582 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
585 return self.__is_scalar
588 return self.__is_vector
591 return self.__is_matrix
594 return self.__is_object
598 return Covariance(self.__name+"I", asCovariance = self.__C.I )
599 elif self.isvector():
600 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
601 elif self.isscalar():
602 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
603 elif self.isobject():
604 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
606 return None # Indispensable
610 return Covariance(self.__name+"T", asCovariance = self.__C.T )
611 elif self.isvector():
612 return Covariance(self.__name+"T", asEyeByVector = self.__C )
613 elif self.isscalar():
614 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
615 elif self.isobject():
616 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
620 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
621 elif self.isvector():
622 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
623 elif self.isscalar():
624 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
625 elif self.isobject() and hasattr(self.__C,"cholesky"):
626 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
630 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
631 elif self.isvector():
632 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
633 elif self.isscalar():
634 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
635 elif self.isobject() and hasattr(self.__C,"choleskyI"):
636 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
638 def diag(self, msize=None):
640 return numpy.diag(self.__C)
641 elif self.isvector():
643 elif self.isscalar():
645 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,))
647 return self.__C * numpy.ones(int(msize))
648 elif self.isobject():
649 return self.__C.diag()
651 def asfullmatrix(self, msize=None):
654 elif self.isvector():
655 return numpy.matrix( numpy.diag(self.__C), float )
656 elif self.isscalar():
658 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,))
660 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
661 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
662 return self.__C.asfullmatrix()
664 def trace(self, msize=None):
666 return numpy.trace(self.__C)
667 elif self.isvector():
668 return float(numpy.sum(self.__C))
669 elif self.isscalar():
671 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,))
673 return self.__C * int(msize)
674 elif self.isobject():
675 return self.__C.trace()
678 return repr(self.__C)
683 def __add__(self, other):
684 if self.ismatrix() or self.isobject():
685 return self.__C + numpy.asmatrix(other)
686 elif self.isvector() or self.isscalar():
687 _A = numpy.asarray(other)
688 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
689 return numpy.asmatrix(_A)
691 def __radd__(self, other):
692 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
694 def __sub__(self, other):
695 if self.ismatrix() or self.isobject():
696 return self.__C - numpy.asmatrix(other)
697 elif self.isvector() or self.isscalar():
698 _A = numpy.asarray(other)
699 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
700 return numpy.asmatrix(_A)
702 def __rsub__(self, other):
703 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
708 def __mul__(self, other):
709 if self.ismatrix() and isinstance(other,numpy.matrix):
710 return self.__C * other
711 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
712 or isinstance(other,list) \
713 or isinstance(other,tuple)):
714 if numpy.ravel(other).size == self.shape[1]: # Vecteur
715 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
716 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
717 return self.__C * numpy.asmatrix(other)
719 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
720 elif self.isvector() and (isinstance(other,numpy.matrix) \
721 or isinstance(other,numpy.ndarray) \
722 or isinstance(other,list) \
723 or isinstance(other,tuple)):
724 if numpy.ravel(other).size == self.shape[1]: # Vecteur
725 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
726 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
727 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
729 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
730 elif self.isscalar() and isinstance(other,numpy.matrix):
731 return self.__C * other
732 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
733 or isinstance(other,list) \
734 or isinstance(other,tuple)):
735 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
736 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
738 return self.__C * numpy.asmatrix(other)
739 elif self.isobject():
740 return self.__C.__mul__(other)
742 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
744 def __rmul__(self, other):
745 if self.ismatrix() and isinstance(other,numpy.matrix):
746 return other * self.__C
747 elif self.isvector() and isinstance(other,numpy.matrix):
748 if numpy.ravel(other).size == self.shape[0]: # Vecteur
749 return numpy.asmatrix(numpy.ravel(other) * self.__C)
750 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
751 return numpy.asmatrix(numpy.array(other) * self.__C)
753 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
754 elif self.isscalar() and isinstance(other,numpy.matrix):
755 return other * self.__C
756 elif self.isobject():
757 return self.__C.__rmul__(other)
759 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
764 # ==============================================================================
765 if __name__ == "__main__":
766 print '\n AUTODIAGNOSTIC \n'