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)
54 self.__listOPCV = [] # Operator Previous Calculated Points, Results, Point Norms
55 # logging.debug("CM Tolerance de determination des doublons : %.2e"%self.__tolerBP)
57 def wasCalculatedIn(self, xValue ): #, info="" ):
60 for i in xrange(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
61 if xValue.size != self.__listOPCV[i][0].size:
62 # 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))
64 if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
66 __HxV = self.__listOPCV[i][1]
67 # logging.debug("CM Cas%s déja calculé, portant le numéro %i"%(info,i))
71 def storeValueInX(self, xValue, HxValue ):
72 if self.__lenghtOR < 0: self.__lenghtOR = 2 * xValue.size + 2
73 while len(self.__listOPCV) > self.__lenghtOR:
74 # logging.debug("CM Réduction de la liste des cas à %i éléments par suppression du premier"%self.__lenghtOR)
75 self.__listOPCV.pop(0)
76 self.__listOPCV.append( (
77 copy.copy(numpy.ravel(xValue)),
79 numpy.linalg.norm(xValue),
82 # ==============================================================================
85 Classe générale d'interface de type opérateur
92 def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
94 On construit un objet de ce type en fournissant à l'aide de l'un des
95 deux mots-clé, soit une fonction python, soit une matrice.
97 - fromMethod : argument de type fonction Python
98 - fromMatrix : argument adapté au constructeur numpy.matrix
99 - avoidingRedundancy : évite ou pas les calculs redondants
101 self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
102 self.__AvoidRC = bool( avoidingRedundancy )
103 if fromMethod is not None:
104 self.__Method = fromMethod
106 self.__Type = "Method"
107 elif fromMatrix is not None:
109 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
110 self.__Type = "Matrix"
119 def appliedTo(self, xValue):
121 Permet de restituer le résultat de l'application de l'opérateur à un
122 argument xValue. Cette méthode se contente d'appliquer, son argument
123 devant a priori être du bon type.
125 - xValue : argument adapté pour appliquer l'opérateur
128 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
130 __alreadyCalculated = False
132 if __alreadyCalculated:
133 self.__addOneCacheCall()
136 if self.__Matrix is not None:
137 self.__addOneMatrixCall()
138 HxValue = self.__Matrix * xValue
140 self.__addOneMethodCall()
141 HxValue = self.__Method( xValue )
143 Operator.CM.storeValueInX(xValue,HxValue)
147 def appliedControledFormTo(self, (xValue, uValue) ):
149 Permet de restituer le résultat de l'application de l'opérateur à une
150 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
151 argument devant a priori être du bon type. Si la uValue est None,
152 on suppose que l'opérateur ne s'applique qu'à xValue.
154 - xValue : argument X adapté pour appliquer l'opérateur
155 - uValue : argument U adapté pour appliquer l'opérateur
157 if self.__Matrix is not None:
158 self.__addOneMatrixCall()
159 return self.__Matrix * xValue
160 elif uValue is not None:
161 self.__addOneMethodCall()
162 return self.__Method( (xValue, uValue) )
164 self.__addOneMethodCall()
165 return self.__Method( xValue )
167 def appliedInXTo(self, (xNominal, xValue) ):
169 Permet de restituer le résultat de l'application de l'opérateur à un
170 argument xValue, sachant que l'opérateur est valable en xNominal.
171 Cette méthode se contente d'appliquer, son argument devant a priori
172 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
173 alors il est valable en tout point nominal et il n'est pas nécessaire
175 Arguments : une liste contenant
176 - xNominal : argument permettant de donner le point où l'opérateur
177 est construit pour etre ensuite appliqué
178 - xValue : argument adapté pour appliquer l'opérateur
180 if self.__Matrix is not None:
181 self.__addOneMatrixCall()
182 return self.__Matrix * xValue
184 self.__addOneMethodCall()
185 return self.__Method( (xNominal, xValue) )
187 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
189 Permet de renvoyer l'opérateur sous la forme d'une matrice
191 if self.__Matrix is not None:
192 self.__addOneMatrixCall()
194 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
195 self.__addOneMethodCall()
196 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
198 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
202 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
203 la forme d'une matrice
205 if self.__Matrix is not None:
206 return self.__Matrix.shape
208 raise ValueError("Matrix form of the operator is not available, nor the shape")
210 def nbcalls(self, which=None):
212 Renvoie les nombres d'évaluations de l'opérateur
215 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
216 self.__NbCallsAsMatrix,
217 self.__NbCallsAsMethod,
218 self.__NbCallsOfCached,
219 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
220 Operator.NbCallsAsMatrix,
221 Operator.NbCallsAsMethod,
222 Operator.NbCallsOfCached,
224 if which is None: return __nbcalls
225 else: return __nbcalls[which]
227 def __addOneMatrixCall(self):
228 self.__NbCallsAsMatrix += 1 # Decompte local
229 Operator.NbCallsAsMatrix += 1 # Decompte global
231 def __addOneMethodCall(self):
232 self.__NbCallsAsMethod += 1 # Decompte local
233 Operator.NbCallsAsMethod += 1 # Decompte global
235 def __addOneCacheCall(self):
236 self.__NbCallsOfCached += 1 # Decompte local
237 Operator.NbCallsOfCached += 1 # Decompte global
239 # ==============================================================================
242 Classe générale d'interface de type algorithme
244 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
245 d'assimilation, en fournissant un container (dictionnaire) de variables
246 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
248 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
250 def __init__(self, name):
252 L'initialisation présente permet de fabriquer des variables de stockage
253 disponibles de manière générique dans les algorithmes élémentaires. Ces
254 variables de stockage sont ensuite conservées dans un dictionnaire
255 interne à l'objet, mais auquel on accède par la méthode "get".
257 Les variables prévues sont :
258 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
259 - CostFunctionJb : partie ébauche ou background de la fonction-cout
260 - CostFunctionJo : partie observations de la fonction-cout
261 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
262 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
263 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
264 - CurrentState : état courant lors d'itérations
265 - Analysis : l'analyse Xa
266 - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
267 - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
268 - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
269 - Innovation : l'innovation : d = Y - H(X)
270 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
271 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
272 - MahalanobisConsistency : indicateur de consistance des covariances
273 - OMA : Observation moins Analysis : Y - Xa
274 - OMB : Observation moins Background : Y - Xb
275 - AMB : Analysis moins Background : Xa - Xb
276 - APosterioriCovariance : matrice A
277 - APosterioriVariances : variances de la matrice A
278 - APosterioriStandardDeviations : écart-types de la matrice A
279 - APosterioriCorrelations : correlations de la matrice A
280 On peut rajouter des variables à stocker dans l'initialisation de
281 l'algorithme élémentaire qui va hériter de cette classe
283 logging.debug("%s Initialisation"%str(name))
284 self._m = PlatformInfo.SystemUsage()
286 self._name = str( name )
287 self._parameters = {"StoreSupplementaryCalculations":[]}
288 self.__required_parameters = {}
289 self.StoredVariables = {}
291 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
292 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
293 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
294 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
295 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
296 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
297 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
298 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
299 self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum")
300 self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
301 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
302 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
303 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
304 self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
305 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
306 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
307 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
308 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
309 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
310 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
311 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
312 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
313 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
314 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
315 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
316 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
319 logging.debug("%s Lancement"%self._name)
320 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
323 def _post_run(self,_oH=None):
324 if self._parameters.has_key("StoreSupplementaryCalculations") and \
325 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
326 for _A in self.StoredVariables["APosterioriCovariance"]:
327 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
328 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
329 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
330 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
331 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
332 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
333 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
334 self.StoredVariables["APosterioriCorrelations"].store( _C )
336 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)))
337 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)))
338 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
339 logging.debug("%s Terminé"%self._name)
342 def get(self, key=None):
344 Renvoie l'une des variables stockées identifiée par la clé, ou le
345 dictionnaire de l'ensemble des variables disponibles en l'absence de
346 clé. Ce sont directement les variables sous forme objet qui sont
347 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
348 des classes de persistance.
351 return self.StoredVariables[key]
353 return self.StoredVariables
355 def has_key(self, key=None):
357 Vérifie si l'une des variables stockées est identifiée par la clé.
359 return self.StoredVariables.has_key(key)
363 Renvoie la liste des clés de variables stockées.
365 return self.StoredVariables.keys()
367 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
369 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
370 sa forme mathématique la plus naturelle possible.
372 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
374 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
376 Permet de définir dans l'algorithme des paramètres requis et leurs
377 caractéristiques par défaut.
380 raise ValueError("A name is mandatory to define a required parameter.")
382 self.__required_parameters[name] = {
384 "typecast" : typecast,
390 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
392 def getRequiredParameters(self, noDetails=True):
394 Renvoie la liste des noms de paramètres requis ou directement le
395 dictionnaire des paramètres requis.
398 ks = self.__required_parameters.keys()
402 return self.__required_parameters
404 def setParameterValue(self, name=None, value=None):
406 Renvoie la valeur d'un paramètre requis de manière contrôlée
408 default = self.__required_parameters[name]["default"]
409 typecast = self.__required_parameters[name]["typecast"]
410 minval = self.__required_parameters[name]["minval"]
411 maxval = self.__required_parameters[name]["maxval"]
412 listval = self.__required_parameters[name]["listval"]
414 if value is None and default is None:
416 elif value is None and default is not None:
417 if typecast is None: __val = default
418 else: __val = typecast( default )
420 if typecast is None: __val = value
421 else: __val = typecast( value )
423 if minval is not None and (numpy.array(__val, float) < minval).any():
424 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
425 if maxval is not None and (numpy.array(__val, float) > maxval).any():
426 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
427 if listval is not None:
428 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
431 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
432 elif __val not in listval:
433 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
436 def setParameters(self, fromDico={}):
438 Permet de stocker les paramètres reçus dans le dictionnaire interne.
440 self._parameters.update( fromDico )
441 for k in self.__required_parameters.keys():
442 if k in fromDico.keys():
443 self._parameters[k] = self.setParameterValue(k,fromDico[k])
445 self._parameters[k] = self.setParameterValue(k)
446 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
448 # ==============================================================================
451 Classe générale d'interface de type diagnostic
453 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
454 même temps que l'une des classes de persistance, comme "OneScalar" par
457 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
458 méthode "_formula" pour écrire explicitement et proprement la formule pour
459 l'écriture mathématique du calcul du diagnostic (méthode interne non
460 publique), et "calculate" pour activer la précédente tout en ayant vérifié
461 et préparé les données, et pour stocker les résultats à chaque pas (méthode
462 externe d'activation).
464 def __init__(self, name = "", parameters = {}):
465 self.name = str(name)
466 self.parameters = dict( parameters )
468 def _formula(self, *args):
470 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
471 mathématique la plus naturelle possible.
473 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
475 def calculate(self, *args):
477 Active la formule de calcul avec les arguments correctement rangés
479 raise NotImplementedError("Diagnostic activation method has not been implemented!")
481 # ==============================================================================
484 Classe générale d'interface de type covariance
487 name = "GenericCovariance",
489 asEyeByScalar = None,
490 asEyeByVector = None,
494 Permet de définir une covariance :
495 - asCovariance : entrée des données, comme une matrice compatible avec
496 le constructeur de numpy.matrix
497 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
498 multiplicatif d'une matrice de corrélation identité, aucune matrice
499 n'étant donc explicitement à donner
500 - asEyeByVector : entrée des données comme un seul vecteur de variance,
501 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
502 n'étant donc explicitement à donner
503 - asCovObject : entrée des données comme un objet python, qui a les
504 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
505 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
506 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
508 self.__name = str(name)
511 self.__is_scalar = False
512 self.__is_vector = False
513 self.__is_matrix = False
514 self.__is_object = False
515 if asEyeByScalar is not None:
516 self.__is_scalar = True
517 self.__C = numpy.abs( float(asEyeByScalar) )
520 elif asEyeByVector is not None:
521 self.__is_vector = True
522 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
523 self.shape = (self.__C.size,self.__C.size)
524 self.size = self.__C.size**2
525 elif asCovariance is not None:
526 self.__is_matrix = True
527 self.__C = numpy.matrix( asCovariance, float )
528 self.shape = self.__C.shape
529 self.size = self.__C.size
530 elif asCovObject is not None:
531 self.__is_object = True
532 self.__C = asCovObject
533 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
534 if not hasattr(self.__C,at):
535 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
536 if hasattr(self.__C,"shape"):
537 self.shape = self.__C.shape
540 if hasattr(self.__C,"size"):
541 self.size = self.__C.size
546 # 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)
550 def __validate(self):
551 if self.ismatrix() and min(self.shape) != max(self.shape):
552 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))
553 if self.isobject() and min(self.shape) != max(self.shape):
554 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))
555 if self.isscalar() and self.__C <= 0:
556 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
557 if self.isvector() and (self.__C <= 0).any():
558 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
559 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
561 L = numpy.linalg.cholesky( self.__C )
563 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
566 return self.__is_scalar
569 return self.__is_vector
572 return self.__is_matrix
575 return self.__is_object
579 return Covariance(self.__name+"I", asCovariance = self.__C.I )
580 elif self.isvector():
581 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
582 elif self.isscalar():
583 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
584 elif self.isobject():
585 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
587 return None # Indispensable
591 return Covariance(self.__name+"T", asCovariance = self.__C.T )
592 elif self.isvector():
593 return Covariance(self.__name+"T", asEyeByVector = self.__C )
594 elif self.isscalar():
595 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
596 elif self.isobject():
597 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
601 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
602 elif self.isvector():
603 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
604 elif self.isscalar():
605 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
606 elif self.isobject() and hasattr(self.__C,"cholesky"):
607 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
611 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
612 elif self.isvector():
613 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
614 elif self.isscalar():
615 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
616 elif self.isobject() and hasattr(self.__C,"choleskyI"):
617 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
619 def diag(self, msize=None):
621 return numpy.diag(self.__C)
622 elif self.isvector():
624 elif self.isscalar():
626 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,))
628 return self.__C * numpy.ones(int(msize))
629 elif self.isobject():
630 return self.__C.diag()
632 def asfullmatrix(self, msize=None):
635 elif self.isvector():
636 return numpy.matrix( numpy.diag(self.__C), float )
637 elif self.isscalar():
639 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,))
641 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
642 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
643 return self.__C.asfullmatrix()
645 def trace(self, msize=None):
647 return numpy.trace(self.__C)
648 elif self.isvector():
649 return float(numpy.sum(self.__C))
650 elif self.isscalar():
652 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,))
654 return self.__C * int(msize)
655 elif self.isobject():
656 return self.__C.trace()
659 return repr(self.__C)
664 def __add__(self, other):
665 if self.ismatrix() or self.isobject():
666 return self.__C + numpy.asmatrix(other)
667 elif self.isvector() or self.isscalar():
668 _A = numpy.asarray(other)
669 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
670 return numpy.asmatrix(_A)
672 def __radd__(self, other):
673 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
675 def __sub__(self, other):
676 if self.ismatrix() or self.isobject():
677 return self.__C - numpy.asmatrix(other)
678 elif self.isvector() or self.isscalar():
679 _A = numpy.asarray(other)
680 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
681 return numpy.asmatrix(_A)
683 def __rsub__(self, other):
684 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
689 def __mul__(self, other):
690 if self.ismatrix() and isinstance(other,numpy.matrix):
691 return self.__C * other
692 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
693 or isinstance(other,list) \
694 or isinstance(other,tuple)):
695 if numpy.ravel(other).size == self.shape[1]: # Vecteur
696 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
697 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
698 return self.__C * numpy.asmatrix(other)
700 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
701 elif self.isvector() and (isinstance(other,numpy.matrix) \
702 or isinstance(other,numpy.ndarray) \
703 or isinstance(other,list) \
704 or isinstance(other,tuple)):
705 if numpy.ravel(other).size == self.shape[1]: # Vecteur
706 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
707 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
708 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
710 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
711 elif self.isscalar() and isinstance(other,numpy.matrix):
712 return self.__C * other
713 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
714 or isinstance(other,list) \
715 or isinstance(other,tuple)):
716 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
717 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
719 return self.__C * numpy.asmatrix(other)
720 elif self.isobject():
721 return self.__C.__mul__(other)
723 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
725 def __rmul__(self, other):
726 if self.ismatrix() and isinstance(other,numpy.matrix):
727 return other * self.__C
728 elif self.isvector() and isinstance(other,numpy.matrix):
729 if numpy.ravel(other).size == self.shape[0]: # Vecteur
730 return numpy.asmatrix(numpy.ravel(other) * self.__C)
731 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
732 return numpy.asmatrix(numpy.array(other) * self.__C)
734 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
735 elif self.isscalar() and isinstance(other,numpy.matrix):
736 return other * self.__C
737 elif self.isobject():
738 return self.__C.__rmul__(other)
740 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
745 # ==============================================================================
746 if __name__ == "__main__":
747 print '\n AUTODIAGNOSTIC \n'