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 On peut rajouter des variables à stocker dans l'initialisation de
277 l'algorithme élémentaire qui va hériter de cette classe
279 logging.debug("%s Initialisation"%str(name))
280 self._m = PlatformInfo.SystemUsage()
282 self._name = str( name )
283 self._parameters = {}
284 self.__required_parameters = {}
285 self.StoredVariables = {}
287 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
288 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
289 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
290 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
291 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
292 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
293 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
294 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
295 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
296 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
297 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
298 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
299 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
300 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
301 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
302 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
303 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
304 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
305 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
306 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
309 logging.debug("%s Lancement"%self._name)
310 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
313 def _post_run(self,_oH=None):
315 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)))
316 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)))
317 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
318 logging.debug("%s Terminé"%self._name)
321 def get(self, key=None):
323 Renvoie l'une des variables stockées identifiée par la clé, ou le
324 dictionnaire de l'ensemble des variables disponibles en l'absence de
325 clé. Ce sont directement les variables sous forme objet qui sont
326 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
327 des classes de persistance.
330 return self.StoredVariables[key]
332 return self.StoredVariables
334 def has_key(self, key=None):
336 Vérifie si l'une des variables stockées est identifiée par la clé.
338 return self.StoredVariables.has_key(key)
342 Renvoie la liste des clés de variables stockées.
344 return self.StoredVariables.keys()
346 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
348 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
349 sa forme mathématique la plus naturelle possible.
351 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
353 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
355 Permet de définir dans l'algorithme des paramètres requis et leurs
356 caractéristiques par défaut.
359 raise ValueError("A name is mandatory to define a required parameter.")
361 self.__required_parameters[name] = {
363 "typecast" : typecast,
369 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
371 def getRequiredParameters(self, noDetails=True):
373 Renvoie la liste des noms de paramètres requis ou directement le
374 dictionnaire des paramètres requis.
377 ks = self.__required_parameters.keys()
381 return self.__required_parameters
383 def setParameterValue(self, name=None, value=None):
385 Renvoie la valeur d'un paramètre requis de manière contrôlée
387 default = self.__required_parameters[name]["default"]
388 typecast = self.__required_parameters[name]["typecast"]
389 minval = self.__required_parameters[name]["minval"]
390 maxval = self.__required_parameters[name]["maxval"]
391 listval = self.__required_parameters[name]["listval"]
393 if value is None and default is None:
395 elif value is None and default is not None:
396 if typecast is None: __val = default
397 else: __val = typecast( default )
399 if typecast is None: __val = value
400 else: __val = typecast( value )
402 if minval is not None and (numpy.array(__val) < minval).any():
403 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
404 if maxval is not None and (numpy.array(__val) > maxval).any():
405 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
406 if listval is not None:
407 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
410 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
411 elif __val not in listval:
412 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
415 def setParameters(self, fromDico={}):
417 Permet de stocker les paramètres reçus dans le dictionnaire interne.
419 self._parameters.update( fromDico )
420 for k in self.__required_parameters.keys():
421 if k in fromDico.keys():
422 self._parameters[k] = self.setParameterValue(k,fromDico[k])
424 self._parameters[k] = self.setParameterValue(k)
425 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
427 # ==============================================================================
430 Classe générale d'interface de type diagnostic
432 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
433 même temps que l'une des classes de persistance, comme "OneScalar" par
436 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
437 méthode "_formula" pour écrire explicitement et proprement la formule pour
438 l'écriture mathématique du calcul du diagnostic (méthode interne non
439 publique), et "calculate" pour activer la précédente tout en ayant vérifié
440 et préparé les données, et pour stocker les résultats à chaque pas (méthode
441 externe d'activation).
443 def __init__(self, name = "", parameters = {}):
444 self.name = str(name)
445 self.parameters = dict( parameters )
447 def _formula(self, *args):
449 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
450 mathématique la plus naturelle possible.
452 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
454 def calculate(self, *args):
456 Active la formule de calcul avec les arguments correctement rangés
458 raise NotImplementedError("Diagnostic activation method has not been implemented!")
460 # ==============================================================================
463 Classe générale d'interface de type covariance
466 name = "GenericCovariance",
468 asEyeByScalar = None,
469 asEyeByVector = None,
473 Permet de définir une covariance :
474 - asCovariance : entrée des données, comme une matrice compatible avec
475 le constructeur de numpy.matrix
476 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
477 multiplicatif d'une matrice de corrélation identité, aucune matrice
478 n'étant donc explicitement à donner
479 - asEyeByVector : entrée des données comme un seul vecteur de variance,
480 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
481 n'étant donc explicitement à donner
482 - asCovObject : entrée des données comme un objet python, qui a les
483 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
484 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
485 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
487 self.__name = str(name)
490 self.__is_scalar = False
491 self.__is_vector = False
492 self.__is_matrix = False
493 self.__is_object = False
494 if asEyeByScalar is not None:
495 self.__is_scalar = True
496 self.__C = numpy.abs( float(asEyeByScalar) )
499 elif asEyeByVector is not None:
500 self.__is_vector = True
501 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
502 self.shape = (self.__C.size,self.__C.size)
503 self.size = self.__C.size**2
504 elif asCovariance is not None:
505 self.__is_matrix = True
506 self.__C = numpy.matrix( asCovariance, float )
507 self.shape = self.__C.shape
508 self.size = self.__C.size
509 elif asCovObject is not None:
510 self.__is_object = True
511 self.__C = asCovObject
512 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
513 if not hasattr(self.__C,at):
514 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
515 if hasattr(self.__C,"shape"):
516 self.shape = self.__C.shape
519 if hasattr(self.__C,"size"):
520 self.size = self.__C.size
525 # 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)
529 def __validate(self):
530 if self.ismatrix() and min(self.shape) != max(self.shape):
531 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))
532 if self.isobject() and min(self.shape) != max(self.shape):
533 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))
534 if self.isscalar() and self.__C <= 0:
535 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
536 if self.isvector() and (self.__C <= 0).any():
537 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
538 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
540 L = numpy.linalg.cholesky( self.__C )
542 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
545 return self.__is_scalar
548 return self.__is_vector
551 return self.__is_matrix
554 return self.__is_object
558 return Covariance(self.__name+"I", asCovariance = self.__C.I )
559 elif self.isvector():
560 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
561 elif self.isscalar():
562 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
563 elif self.isobject():
564 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
566 return None # Indispensable
570 return Covariance(self.__name+"T", asCovariance = self.__C.T )
571 elif self.isvector():
572 return Covariance(self.__name+"T", asEyeByVector = self.__C )
573 elif self.isscalar():
574 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
575 elif self.isobject():
576 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
580 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
581 elif self.isvector():
582 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
583 elif self.isscalar():
584 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
585 elif self.isobject() and hasattr(self.__C,"cholesky"):
586 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
590 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
591 elif self.isvector():
592 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
593 elif self.isscalar():
594 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
595 elif self.isobject() and hasattr(self.__C,"choleskyI"):
596 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
598 def diag(self, msize=None):
600 return numpy.diag(self.__C)
601 elif self.isvector():
603 elif self.isscalar():
605 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,))
607 return self.__C * numpy.ones(int(msize))
608 elif self.isobject():
609 return self.__C.diag()
611 def asfullmatrix(self, msize=None):
614 elif self.isvector():
615 return numpy.matrix( numpy.diag(self.__C), float )
616 elif self.isscalar():
618 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,))
620 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
621 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
622 return self.__C.asfullmatrix()
624 def trace(self, msize=None):
626 return numpy.trace(self.__C)
627 elif self.isvector():
628 return float(numpy.sum(self.__C))
629 elif self.isscalar():
631 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,))
633 return self.__C * int(msize)
634 elif self.isobject():
635 return self.__C.trace()
638 return repr(self.__C)
643 def __add__(self, other):
644 if self.ismatrix() or self.isobject():
645 return self.__C + numpy.asmatrix(other)
646 elif self.isvector() or self.isscalar():
647 _A = numpy.asarray(other)
648 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
649 return numpy.asmatrix(_A)
651 def __radd__(self, other):
652 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
654 def __sub__(self, other):
655 if self.ismatrix() or self.isobject():
656 return self.__C - numpy.asmatrix(other)
657 elif self.isvector() or self.isscalar():
658 _A = numpy.asarray(other)
659 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
660 return numpy.asmatrix(_A)
662 def __rsub__(self, other):
663 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
668 def __mul__(self, other):
669 if self.ismatrix() and isinstance(other,numpy.matrix):
670 return self.__C * other
671 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
672 or isinstance(other,list) \
673 or isinstance(other,tuple)):
674 if numpy.ravel(other).size == self.shape[1]: # Vecteur
675 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
676 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
677 return self.__C * numpy.asmatrix(other)
679 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
680 elif self.isvector() and (isinstance(other,numpy.matrix) \
681 or isinstance(other,numpy.ndarray) \
682 or isinstance(other,list) \
683 or isinstance(other,tuple)):
684 if numpy.ravel(other).size == self.shape[1]: # Vecteur
685 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
686 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
687 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
689 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
690 elif self.isscalar() and isinstance(other,numpy.matrix):
691 return self.__C * other
692 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
693 or isinstance(other,list) \
694 or isinstance(other,tuple)):
695 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
696 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
698 return self.__C * numpy.asmatrix(other)
699 elif self.isobject():
700 return self.__C.__mul__(other)
702 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
704 def __rmul__(self, other):
705 if self.ismatrix() and isinstance(other,numpy.matrix):
706 return other * self.__C
707 elif self.isvector() and isinstance(other,numpy.matrix):
708 if numpy.ravel(other).size == self.shape[0]: # Vecteur
709 return numpy.asmatrix(numpy.ravel(other) * self.__C)
710 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
711 return numpy.asmatrix(numpy.array(other) * self.__C)
713 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
714 elif self.isscalar() and isinstance(other,numpy.matrix):
715 return other * self.__C
716 elif self.isobject():
717 return self.__C.__rmul__(other)
719 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
724 # ==============================================================================
725 if __name__ == "__main__":
726 print '\n AUTODIAGNOSTIC \n'