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 matrice.
97 - fromMethod : argument de type fonction Python
98 - fromMatrix : argument adapté au constructeur numpy.matrix
100 self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
101 self.__AvoidRC = bool( avoidingRedundancy )
102 if fromMethod is not None:
103 self.__Method = fromMethod
105 self.__Type = "Method"
106 elif fromMatrix is not None:
108 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
109 self.__Type = "Matrix"
118 def appliedTo(self, xValue):
120 Permet de restituer le résultat de l'application de l'opérateur à un
121 argument xValue. Cette méthode se contente d'appliquer, son argument
122 devant a priori être du bon type.
124 - xValue : argument adapté pour appliquer l'opérateur
127 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
129 __alreadyCalculated = False
131 if __alreadyCalculated:
132 self.__addOneCacheCall()
135 if self.__Matrix is not None:
136 self.__addOneMatrixCall()
137 HxValue = self.__Matrix * xValue
139 self.__addOneMethodCall()
140 HxValue = self.__Method( xValue )
142 Operator.CM.storeValueInX(xValue,HxValue)
146 def appliedControledFormTo(self, (xValue, uValue) ):
148 Permet de restituer le résultat de l'application de l'opérateur à une
149 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
150 argument devant a priori être du bon type. Si la uValue est None,
151 on suppose que l'opérateur ne s'applique qu'à xValue.
153 - xValue : argument X adapté pour appliquer l'opérateur
154 - uValue : argument U adapté pour appliquer l'opérateur
156 if self.__Matrix is not None:
157 self.__addOneMatrixCall()
158 return self.__Matrix * xValue
159 elif uValue is not None:
160 self.__addOneMethodCall()
161 return self.__Method( (xValue, uValue) )
163 self.__addOneMethodCall()
164 return self.__Method( xValue )
166 def appliedInXTo(self, (xNominal, xValue) ):
168 Permet de restituer le résultat de l'application de l'opérateur à un
169 argument xValue, sachant que l'opérateur est valable en xNominal.
170 Cette méthode se contente d'appliquer, son argument devant a priori
171 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
172 alors il est valable en tout point nominal et il n'est pas nécessaire
174 Arguments : une liste contenant
175 - xNominal : argument permettant de donner le point où l'opérateur
176 est construit pour etre ensuite appliqué
177 - xValue : argument adapté pour appliquer l'opérateur
179 if self.__Matrix is not None:
180 self.__addOneMatrixCall()
181 return self.__Matrix * xValue
183 self.__addOneMethodCall()
184 return self.__Method( (xNominal, xValue) )
186 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
188 Permet de renvoyer l'opérateur sous la forme d'une matrice
190 if self.__Matrix is not None:
191 self.__addOneMatrixCall()
193 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
194 self.__addOneMethodCall()
195 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
197 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
201 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
202 la forme d'une matrice
204 if self.__Matrix is not None:
205 return self.__Matrix.shape
207 raise ValueError("Matrix form of the operator is not available, nor the shape")
209 def nbcalls(self, which=None):
211 Renvoie les nombres d'évaluations de l'opérateur
214 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
215 self.__NbCallsAsMatrix,
216 self.__NbCallsAsMethod,
217 self.__NbCallsOfCached,
218 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
219 Operator.NbCallsAsMatrix,
220 Operator.NbCallsAsMethod,
221 Operator.NbCallsOfCached,
223 if which is None: return __nbcalls
224 else: return __nbcalls[which]
226 def __addOneMatrixCall(self):
227 self.__NbCallsAsMatrix += 1 # Decompte local
228 Operator.NbCallsAsMatrix += 1 # Decompte global
230 def __addOneMethodCall(self):
231 self.__NbCallsAsMethod += 1 # Decompte local
232 Operator.NbCallsAsMethod += 1 # Decompte global
234 def __addOneCacheCall(self):
235 self.__NbCallsOfCached += 1 # Decompte local
236 Operator.NbCallsOfCached += 1 # Decompte global
238 # ==============================================================================
241 Classe générale d'interface de type algorithme
243 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
244 d'assimilation, en fournissant un container (dictionnaire) de variables
245 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
247 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
249 def __init__(self, name):
251 L'initialisation présente permet de fabriquer des variables de stockage
252 disponibles de manière générique dans les algorithmes élémentaires. Ces
253 variables de stockage sont ensuite conservées dans un dictionnaire
254 interne à l'objet, mais auquel on accède par la méthode "get".
256 Les variables prévues sont :
257 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
258 - CostFunctionJb : partie ébauche ou background de la fonction-cout
259 - CostFunctionJo : partie observations de la fonction-cout
260 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
261 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
262 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
263 - CurrentState : état courant lors d'itérations
264 - Analysis : l'analyse Xa
265 - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
266 - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
267 - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
268 - Innovation : l'innovation : d = Y - H(X)
269 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
270 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
271 - MahalanobisConsistency : indicateur de consistance des covariances
272 - OMA : Observation moins Analysis : Y - Xa
273 - OMB : Observation moins Background : Y - Xb
274 - AMB : Analysis moins Background : Xa - Xb
275 - APosterioriCovariance : matrice A
276 - APosterioriVariances : variances de la matrice A
277 - APosterioriStandardDeviations : écart-types de la matrice A
278 - APosterioriCorrelations : correlations de la matrice A
279 On peut rajouter des variables à stocker dans l'initialisation de
280 l'algorithme élémentaire qui va hériter de cette classe
282 logging.debug("%s Initialisation"%str(name))
283 self._m = PlatformInfo.SystemUsage()
285 self._name = str( name )
286 self._parameters = {}
287 self.__required_parameters = {}
288 self.StoredVariables = {}
290 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
291 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
292 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
293 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
294 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
295 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
296 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
297 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
298 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
299 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
300 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
301 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
302 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
303 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
304 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
305 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
306 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
307 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
308 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
309 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
310 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
311 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
312 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
315 logging.debug("%s Lancement"%self._name)
316 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
319 def _post_run(self,_oH=None):
320 if self._parameters.has_key("StoreSupplementaryCalculations") and \
321 "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'