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["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
300 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
301 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
302 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
303 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
304 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
305 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
306 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
307 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
308 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
309 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
310 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
311 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
312 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
313 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
316 logging.debug("%s Lancement"%self._name)
317 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
320 def _post_run(self,_oH=None):
321 if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
322 for _A in self.StoredVariables["APosterioriCovariance"]:
323 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
324 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
325 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
326 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
327 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
328 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
329 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
330 self.StoredVariables["APosterioriCorrelations"].store( _C )
332 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)))
333 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)))
334 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
335 logging.debug("%s Terminé"%self._name)
338 def get(self, key=None):
340 Renvoie l'une des variables stockées identifiée par la clé, ou le
341 dictionnaire de l'ensemble des variables disponibles en l'absence de
342 clé. Ce sont directement les variables sous forme objet qui sont
343 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
344 des classes de persistance.
347 return self.StoredVariables[key]
349 return self.StoredVariables
351 def has_key(self, key=None):
353 Vérifie si l'une des variables stockées est identifiée par la clé.
355 return self.StoredVariables.has_key(key)
359 Renvoie la liste des clés de variables stockées.
361 return self.StoredVariables.keys()
363 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
365 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
366 sa forme mathématique la plus naturelle possible.
368 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
370 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
372 Permet de définir dans l'algorithme des paramètres requis et leurs
373 caractéristiques par défaut.
376 raise ValueError("A name is mandatory to define a required parameter.")
378 self.__required_parameters[name] = {
380 "typecast" : typecast,
386 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
388 def getRequiredParameters(self, noDetails=True):
390 Renvoie la liste des noms de paramètres requis ou directement le
391 dictionnaire des paramètres requis.
394 ks = self.__required_parameters.keys()
398 return self.__required_parameters
400 def setParameterValue(self, name=None, value=None):
402 Renvoie la valeur d'un paramètre requis de manière contrôlée
404 default = self.__required_parameters[name]["default"]
405 typecast = self.__required_parameters[name]["typecast"]
406 minval = self.__required_parameters[name]["minval"]
407 maxval = self.__required_parameters[name]["maxval"]
408 listval = self.__required_parameters[name]["listval"]
410 if value is None and default is None:
412 elif value is None and default is not None:
413 if typecast is None: __val = default
414 else: __val = typecast( default )
416 if typecast is None: __val = value
417 else: __val = typecast( value )
419 if minval is not None and (numpy.array(__val) < minval).any():
420 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
421 if maxval is not None and (numpy.array(__val) > maxval).any():
422 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
423 if listval is not None:
424 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
427 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
428 elif __val not in listval:
429 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
432 def setParameters(self, fromDico={}):
434 Permet de stocker les paramètres reçus dans le dictionnaire interne.
436 self._parameters.update( fromDico )
437 for k in self.__required_parameters.keys():
438 if k in fromDico.keys():
439 self._parameters[k] = self.setParameterValue(k,fromDico[k])
441 self._parameters[k] = self.setParameterValue(k)
442 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
444 # ==============================================================================
447 Classe générale d'interface de type diagnostic
449 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
450 même temps que l'une des classes de persistance, comme "OneScalar" par
453 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
454 méthode "_formula" pour écrire explicitement et proprement la formule pour
455 l'écriture mathématique du calcul du diagnostic (méthode interne non
456 publique), et "calculate" pour activer la précédente tout en ayant vérifié
457 et préparé les données, et pour stocker les résultats à chaque pas (méthode
458 externe d'activation).
460 def __init__(self, name = "", parameters = {}):
461 self.name = str(name)
462 self.parameters = dict( parameters )
464 def _formula(self, *args):
466 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
467 mathématique la plus naturelle possible.
469 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
471 def calculate(self, *args):
473 Active la formule de calcul avec les arguments correctement rangés
475 raise NotImplementedError("Diagnostic activation method has not been implemented!")
477 # ==============================================================================
480 Classe générale d'interface de type covariance
483 name = "GenericCovariance",
485 asEyeByScalar = None,
486 asEyeByVector = None,
490 Permet de définir une covariance :
491 - asCovariance : entrée des données, comme une matrice compatible avec
492 le constructeur de numpy.matrix
493 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
494 multiplicatif d'une matrice de corrélation identité, aucune matrice
495 n'étant donc explicitement à donner
496 - asEyeByVector : entrée des données comme un seul vecteur de variance,
497 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
498 n'étant donc explicitement à donner
499 - asCovObject : entrée des données comme un objet python, qui a les
500 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
501 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
502 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
504 self.__name = str(name)
507 self.__is_scalar = False
508 self.__is_vector = False
509 self.__is_matrix = False
510 self.__is_object = False
511 if asEyeByScalar is not None:
512 self.__is_scalar = True
513 self.__C = numpy.abs( float(asEyeByScalar) )
516 elif asEyeByVector is not None:
517 self.__is_vector = True
518 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
519 self.shape = (self.__C.size,self.__C.size)
520 self.size = self.__C.size**2
521 elif asCovariance is not None:
522 self.__is_matrix = True
523 self.__C = numpy.matrix( asCovariance, float )
524 self.shape = self.__C.shape
525 self.size = self.__C.size
526 elif asCovObject is not None:
527 self.__is_object = True
528 self.__C = asCovObject
529 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
530 if not hasattr(self.__C,at):
531 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
532 if hasattr(self.__C,"shape"):
533 self.shape = self.__C.shape
536 if hasattr(self.__C,"size"):
537 self.size = self.__C.size
542 # 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)
546 def __validate(self):
547 if self.ismatrix() and min(self.shape) != max(self.shape):
548 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))
549 if self.isobject() and min(self.shape) != max(self.shape):
550 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))
551 if self.isscalar() and self.__C <= 0:
552 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
553 if self.isvector() and (self.__C <= 0).any():
554 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
555 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
557 L = numpy.linalg.cholesky( self.__C )
559 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
562 return self.__is_scalar
565 return self.__is_vector
568 return self.__is_matrix
571 return self.__is_object
575 return Covariance(self.__name+"I", asCovariance = self.__C.I )
576 elif self.isvector():
577 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
578 elif self.isscalar():
579 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
580 elif self.isobject():
581 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
583 return None # Indispensable
587 return Covariance(self.__name+"T", asCovariance = self.__C.T )
588 elif self.isvector():
589 return Covariance(self.__name+"T", asEyeByVector = self.__C )
590 elif self.isscalar():
591 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
592 elif self.isobject():
593 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
597 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
598 elif self.isvector():
599 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
600 elif self.isscalar():
601 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
602 elif self.isobject() and hasattr(self.__C,"cholesky"):
603 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
607 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
608 elif self.isvector():
609 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
610 elif self.isscalar():
611 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
612 elif self.isobject() and hasattr(self.__C,"choleskyI"):
613 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
615 def diag(self, msize=None):
617 return numpy.diag(self.__C)
618 elif self.isvector():
620 elif self.isscalar():
622 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,))
624 return self.__C * numpy.ones(int(msize))
625 elif self.isobject():
626 return self.__C.diag()
628 def asfullmatrix(self, msize=None):
631 elif self.isvector():
632 return numpy.matrix( numpy.diag(self.__C), float )
633 elif self.isscalar():
635 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,))
637 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
638 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
639 return self.__C.asfullmatrix()
641 def trace(self, msize=None):
643 return numpy.trace(self.__C)
644 elif self.isvector():
645 return float(numpy.sum(self.__C))
646 elif self.isscalar():
648 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,))
650 return self.__C * int(msize)
651 elif self.isobject():
652 return self.__C.trace()
655 return repr(self.__C)
660 def __add__(self, other):
661 if self.ismatrix() or self.isobject():
662 return self.__C + numpy.asmatrix(other)
663 elif self.isvector() or self.isscalar():
664 _A = numpy.asarray(other)
665 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
666 return numpy.asmatrix(_A)
668 def __radd__(self, other):
669 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
671 def __sub__(self, other):
672 if self.ismatrix() or self.isobject():
673 return self.__C - numpy.asmatrix(other)
674 elif self.isvector() or self.isscalar():
675 _A = numpy.asarray(other)
676 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
677 return numpy.asmatrix(_A)
679 def __rsub__(self, other):
680 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
685 def __mul__(self, other):
686 if self.ismatrix() and isinstance(other,numpy.matrix):
687 return self.__C * other
688 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
689 or isinstance(other,list) \
690 or isinstance(other,tuple)):
691 if numpy.ravel(other).size == self.shape[1]: # Vecteur
692 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
693 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
694 return self.__C * numpy.asmatrix(other)
696 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
697 elif self.isvector() and (isinstance(other,numpy.matrix) \
698 or isinstance(other,numpy.ndarray) \
699 or isinstance(other,list) \
700 or isinstance(other,tuple)):
701 if numpy.ravel(other).size == self.shape[1]: # Vecteur
702 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
703 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
704 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
706 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
707 elif self.isscalar() and isinstance(other,numpy.matrix):
708 return self.__C * other
709 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
710 or isinstance(other,list) \
711 or isinstance(other,tuple)):
712 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
713 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
715 return self.__C * numpy.asmatrix(other)
716 elif self.isobject():
717 return self.__C.__mul__(other)
719 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
721 def __rmul__(self, other):
722 if self.ismatrix() and isinstance(other,numpy.matrix):
723 return other * self.__C
724 elif self.isvector() and isinstance(other,numpy.matrix):
725 if numpy.ravel(other).size == self.shape[0]: # Vecteur
726 return numpy.asmatrix(numpy.ravel(other) * self.__C)
727 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
728 return numpy.asmatrix(numpy.array(other) * self.__C)
730 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
731 elif self.isscalar() and isinstance(other,numpy.matrix):
732 return other * self.__C
733 elif self.isobject():
734 return self.__C.__rmul__(other)
736 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
741 # ==============================================================================
742 if __name__ == "__main__":
743 print '\n AUTODIAGNOSTIC \n'