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 = {}
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 self._parameters.has_key("StoreSupplementaryCalculations") and \
322 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
323 for _A in self.StoredVariables["APosterioriCovariance"]:
324 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
325 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
326 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
327 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
328 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
329 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
330 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
331 self.StoredVariables["APosterioriCorrelations"].store( _C )
333 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)))
334 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)))
335 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
336 logging.debug("%s Terminé"%self._name)
339 def get(self, key=None):
341 Renvoie l'une des variables stockées identifiée par la clé, ou le
342 dictionnaire de l'ensemble des variables disponibles en l'absence de
343 clé. Ce sont directement les variables sous forme objet qui sont
344 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
345 des classes de persistance.
348 return self.StoredVariables[key]
350 return self.StoredVariables
352 def has_key(self, key=None):
354 Vérifie si l'une des variables stockées est identifiée par la clé.
356 return self.StoredVariables.has_key(key)
360 Renvoie la liste des clés de variables stockées.
362 return self.StoredVariables.keys()
364 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
366 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
367 sa forme mathématique la plus naturelle possible.
369 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
371 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
373 Permet de définir dans l'algorithme des paramètres requis et leurs
374 caractéristiques par défaut.
377 raise ValueError("A name is mandatory to define a required parameter.")
379 self.__required_parameters[name] = {
381 "typecast" : typecast,
387 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
389 def getRequiredParameters(self, noDetails=True):
391 Renvoie la liste des noms de paramètres requis ou directement le
392 dictionnaire des paramètres requis.
395 ks = self.__required_parameters.keys()
399 return self.__required_parameters
401 def setParameterValue(self, name=None, value=None):
403 Renvoie la valeur d'un paramètre requis de manière contrôlée
405 default = self.__required_parameters[name]["default"]
406 typecast = self.__required_parameters[name]["typecast"]
407 minval = self.__required_parameters[name]["minval"]
408 maxval = self.__required_parameters[name]["maxval"]
409 listval = self.__required_parameters[name]["listval"]
411 if value is None and default is None:
413 elif value is None and default is not None:
414 if typecast is None: __val = default
415 else: __val = typecast( default )
417 if typecast is None: __val = value
418 else: __val = typecast( value )
420 if minval is not None and (numpy.array(__val) < minval).any():
421 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
422 if maxval is not None and (numpy.array(__val) > maxval).any():
423 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
424 if listval is not None:
425 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
428 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
429 elif __val not in listval:
430 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
433 def setParameters(self, fromDico={}):
435 Permet de stocker les paramètres reçus dans le dictionnaire interne.
437 self._parameters.update( fromDico )
438 for k in self.__required_parameters.keys():
439 if k in fromDico.keys():
440 self._parameters[k] = self.setParameterValue(k,fromDico[k])
442 self._parameters[k] = self.setParameterValue(k)
443 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
445 # ==============================================================================
448 Classe générale d'interface de type diagnostic
450 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
451 même temps que l'une des classes de persistance, comme "OneScalar" par
454 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
455 méthode "_formula" pour écrire explicitement et proprement la formule pour
456 l'écriture mathématique du calcul du diagnostic (méthode interne non
457 publique), et "calculate" pour activer la précédente tout en ayant vérifié
458 et préparé les données, et pour stocker les résultats à chaque pas (méthode
459 externe d'activation).
461 def __init__(self, name = "", parameters = {}):
462 self.name = str(name)
463 self.parameters = dict( parameters )
465 def _formula(self, *args):
467 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
468 mathématique la plus naturelle possible.
470 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
472 def calculate(self, *args):
474 Active la formule de calcul avec les arguments correctement rangés
476 raise NotImplementedError("Diagnostic activation method has not been implemented!")
478 # ==============================================================================
481 Classe générale d'interface de type covariance
484 name = "GenericCovariance",
486 asEyeByScalar = None,
487 asEyeByVector = None,
491 Permet de définir une covariance :
492 - asCovariance : entrée des données, comme une matrice compatible avec
493 le constructeur de numpy.matrix
494 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
495 multiplicatif d'une matrice de corrélation identité, aucune matrice
496 n'étant donc explicitement à donner
497 - asEyeByVector : entrée des données comme un seul vecteur de variance,
498 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
499 n'étant donc explicitement à donner
500 - asCovObject : entrée des données comme un objet python, qui a les
501 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
502 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
503 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
505 self.__name = str(name)
508 self.__is_scalar = False
509 self.__is_vector = False
510 self.__is_matrix = False
511 self.__is_object = False
512 if asEyeByScalar is not None:
513 self.__is_scalar = True
514 self.__C = numpy.abs( float(asEyeByScalar) )
517 elif asEyeByVector is not None:
518 self.__is_vector = True
519 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
520 self.shape = (self.__C.size,self.__C.size)
521 self.size = self.__C.size**2
522 elif asCovariance is not None:
523 self.__is_matrix = True
524 self.__C = numpy.matrix( asCovariance, float )
525 self.shape = self.__C.shape
526 self.size = self.__C.size
527 elif asCovObject is not None:
528 self.__is_object = True
529 self.__C = asCovObject
530 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
531 if not hasattr(self.__C,at):
532 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
533 if hasattr(self.__C,"shape"):
534 self.shape = self.__C.shape
537 if hasattr(self.__C,"size"):
538 self.size = self.__C.size
543 # 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)
547 def __validate(self):
548 if self.ismatrix() and min(self.shape) != max(self.shape):
549 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))
550 if self.isobject() and min(self.shape) != max(self.shape):
551 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))
552 if self.isscalar() and self.__C <= 0:
553 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
554 if self.isvector() and (self.__C <= 0).any():
555 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
556 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
558 L = numpy.linalg.cholesky( self.__C )
560 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
563 return self.__is_scalar
566 return self.__is_vector
569 return self.__is_matrix
572 return self.__is_object
576 return Covariance(self.__name+"I", asCovariance = self.__C.I )
577 elif self.isvector():
578 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
579 elif self.isscalar():
580 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
581 elif self.isobject():
582 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
584 return None # Indispensable
588 return Covariance(self.__name+"T", asCovariance = self.__C.T )
589 elif self.isvector():
590 return Covariance(self.__name+"T", asEyeByVector = self.__C )
591 elif self.isscalar():
592 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
593 elif self.isobject():
594 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
598 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
599 elif self.isvector():
600 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
601 elif self.isscalar():
602 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
603 elif self.isobject() and hasattr(self.__C,"cholesky"):
604 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
608 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
609 elif self.isvector():
610 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
611 elif self.isscalar():
612 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
613 elif self.isobject() and hasattr(self.__C,"choleskyI"):
614 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
616 def diag(self, msize=None):
618 return numpy.diag(self.__C)
619 elif self.isvector():
621 elif self.isscalar():
623 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,))
625 return self.__C * numpy.ones(int(msize))
626 elif self.isobject():
627 return self.__C.diag()
629 def asfullmatrix(self, msize=None):
632 elif self.isvector():
633 return numpy.matrix( numpy.diag(self.__C), float )
634 elif self.isscalar():
636 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,))
638 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
639 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
640 return self.__C.asfullmatrix()
642 def trace(self, msize=None):
644 return numpy.trace(self.__C)
645 elif self.isvector():
646 return float(numpy.sum(self.__C))
647 elif self.isscalar():
649 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,))
651 return self.__C * int(msize)
652 elif self.isobject():
653 return self.__C.trace()
656 return repr(self.__C)
661 def __add__(self, other):
662 if self.ismatrix() or self.isobject():
663 return self.__C + numpy.asmatrix(other)
664 elif self.isvector() or self.isscalar():
665 _A = numpy.asarray(other)
666 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
667 return numpy.asmatrix(_A)
669 def __radd__(self, other):
670 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
672 def __sub__(self, other):
673 if self.ismatrix() or self.isobject():
674 return self.__C - numpy.asmatrix(other)
675 elif self.isvector() or self.isscalar():
676 _A = numpy.asarray(other)
677 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
678 return numpy.asmatrix(_A)
680 def __rsub__(self, other):
681 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
686 def __mul__(self, other):
687 if self.ismatrix() and isinstance(other,numpy.matrix):
688 return self.__C * other
689 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
690 or isinstance(other,list) \
691 or isinstance(other,tuple)):
692 if numpy.ravel(other).size == self.shape[1]: # Vecteur
693 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
694 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
695 return self.__C * numpy.asmatrix(other)
697 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
698 elif self.isvector() and (isinstance(other,numpy.matrix) \
699 or isinstance(other,numpy.ndarray) \
700 or isinstance(other,list) \
701 or isinstance(other,tuple)):
702 if numpy.ravel(other).size == self.shape[1]: # Vecteur
703 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
704 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
705 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
707 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
708 elif self.isscalar() and isinstance(other,numpy.matrix):
709 return self.__C * other
710 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
711 or isinstance(other,list) \
712 or isinstance(other,tuple)):
713 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
714 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
716 return self.__C * numpy.asmatrix(other)
717 elif self.isobject():
718 return self.__C.__mul__(other)
720 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
722 def __rmul__(self, other):
723 if self.ismatrix() and isinstance(other,numpy.matrix):
724 return other * self.__C
725 elif self.isvector() and isinstance(other,numpy.matrix):
726 if numpy.ravel(other).size == self.shape[0]: # Vecteur
727 return numpy.asmatrix(numpy.ravel(other) * self.__C)
728 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
729 return numpy.asmatrix(numpy.array(other) * self.__C)
731 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
732 elif self.isscalar() and isinstance(other,numpy.matrix):
733 return other * self.__C
734 elif self.isobject():
735 return self.__C.__rmul__(other)
737 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
742 # ==============================================================================
743 if __name__ == "__main__":
744 print '\n AUTODIAGNOSTIC \n'