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"
36 # ==============================================================================
39 Classe générale de gestion d'un cache de calculs
42 toleranceInRedundancy = 1.e-18,
43 lenghtOfRedundancy = -1,
46 Les caractéristiques de tolérance peuvent être modifées à la création.
48 self.__tolerBP = float(toleranceInRedundancy)
49 self.__lenghtOR = int(lenghtOfRedundancy)
53 self.__listOPCV = [] # Operator Previous Calculated Points, Results, Point Norms
54 # logging.debug("CM Tolerance de determination des doublons : %.2e"%self.__tolerBP)
56 def wasCalculatedIn(self, xValue ): #, info="" ):
59 for i in xrange(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
60 if xValue.size != self.__listOPCV[i][0].size:
61 # 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))
63 if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
65 __HxV = self.__listOPCV[i][1]
66 # logging.debug("CM Cas%s déja calculé, portant le numéro %i"%(info,i))
70 def storeValueInX(self, xValue, HxValue ):
71 if self.__lenghtOR < 0: self.__lenghtOR = 2 * xValue.size + 2
72 while len(self.__listOPCV) > self.__lenghtOR:
73 # logging.debug("CM Réduction de la liste des cas à %i éléments par suppression du premier"%self.__lenghtOR)
74 self.__listOPCV.pop(0)
75 self.__listOPCV.append( (
76 copy.copy(numpy.ravel(xValue)),
78 numpy.linalg.norm(xValue),
81 # ==============================================================================
84 Classe générale d'interface de type opérateur
91 def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
93 On construit un objet de ce type en fournissant à l'aide de l'un des
94 deux mots-clé, soit une fonction python, soit matrice.
96 - fromMethod : argument de type fonction Python
97 - fromMatrix : argument adapté au constructeur numpy.matrix
99 self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
100 self.__AvoidRC = bool( avoidingRedundancy )
101 if fromMethod is not None:
102 self.__Method = fromMethod
104 self.__Type = "Method"
105 elif fromMatrix is not None:
107 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
108 self.__Type = "Matrix"
117 def appliedTo(self, xValue):
119 Permet de restituer le résultat de l'application de l'opérateur à un
120 argument xValue. Cette méthode se contente d'appliquer, son argument
121 devant a priori être du bon type.
123 - xValue : argument adapté pour appliquer l'opérateur
126 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
128 __alreadyCalculated = False
130 if __alreadyCalculated:
131 self.__addOneCacheCall()
134 if self.__Matrix is not None:
135 self.__addOneMatrixCall()
136 HxValue = self.__Matrix * xValue
138 self.__addOneMethodCall()
139 HxValue = self.__Method( xValue )
141 Operator.CM.storeValueInX(xValue,HxValue)
145 def appliedControledFormTo(self, (xValue, uValue) ):
147 Permet de restituer le résultat de l'application de l'opérateur à une
148 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
149 argument devant a priori être du bon type. Si la uValue est None,
150 on suppose que l'opérateur ne s'applique qu'à xValue.
152 - xValue : argument X adapté pour appliquer l'opérateur
153 - uValue : argument U adapté pour appliquer l'opérateur
155 if self.__Matrix is not None:
156 self.__addOneMatrixCall()
157 return self.__Matrix * xValue
158 elif uValue is not None:
159 self.__addOneMethodCall()
160 return self.__Method( (xValue, uValue) )
162 self.__addOneMethodCall()
163 return self.__Method( xValue )
165 def appliedInXTo(self, (xNominal, xValue) ):
167 Permet de restituer le résultat de l'application de l'opérateur à un
168 argument xValue, sachant que l'opérateur est valable en xNominal.
169 Cette méthode se contente d'appliquer, son argument devant a priori
170 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
171 alors il est valable en tout point nominal et il n'est pas nécessaire
173 Arguments : une liste contenant
174 - xNominal : argument permettant de donner le point où l'opérateur
175 est construit pour etre ensuite appliqué
176 - xValue : argument adapté pour appliquer l'opérateur
178 if self.__Matrix is not None:
179 self.__addOneMatrixCall()
180 return self.__Matrix * xValue
182 self.__addOneMethodCall()
183 return self.__Method( (xNominal, xValue) )
185 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
187 Permet de renvoyer l'opérateur sous la forme d'une matrice
189 if self.__Matrix is not None:
190 self.__addOneMatrixCall()
192 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
193 self.__addOneMethodCall()
194 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
196 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
200 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
201 la forme d'une matrice
203 if self.__Matrix is not None:
204 return self.__Matrix.shape
206 raise ValueError("Matrix form of the operator is not available, nor the shape")
208 def nbcalls(self, which=None):
210 Renvoie les nombres d'évaluations de l'opérateur
213 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
214 self.__NbCallsAsMatrix,
215 self.__NbCallsAsMethod,
216 self.__NbCallsOfCached,
217 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
218 Operator.NbCallsAsMatrix,
219 Operator.NbCallsAsMethod,
220 Operator.NbCallsOfCached,
222 if which is None: return __nbcalls
223 else: return __nbcalls[which]
225 def __addOneMatrixCall(self):
226 self.__NbCallsAsMatrix += 1 # Decompte local
227 Operator.NbCallsAsMatrix += 1 # Decompte global
229 def __addOneMethodCall(self):
230 self.__NbCallsAsMethod += 1 # Decompte local
231 Operator.NbCallsAsMethod += 1 # Decompte global
233 def __addOneCacheCall(self):
234 self.__NbCallsOfCached += 1 # Decompte local
235 Operator.NbCallsOfCached += 1 # Decompte global
237 # ==============================================================================
240 Classe générale d'interface de type algorithme
242 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
243 d'assimilation, en fournissant un container (dictionnaire) de variables
244 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
246 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
248 def __init__(self, name):
250 L'initialisation présente permet de fabriquer des variables de stockage
251 disponibles de manière générique dans les algorithmes élémentaires. Ces
252 variables de stockage sont ensuite conservées dans un dictionnaire
253 interne à l'objet, mais auquel on accède par la méthode "get".
255 Les variables prévues sont :
256 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
257 - CostFunctionJb : partie ébauche ou background de la fonction-cout
258 - CostFunctionJo : partie observations de la fonction-cout
259 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
260 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
261 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
262 - CurrentState : état courant lors d'itérations
263 - Analysis : l'analyse Xa
264 - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
265 - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
266 - ObservedState : l'état observé H(X)
267 - Innovation : l'innovation : d = Y - H(X)
268 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
269 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
270 - MahalanobisConsistency : indicateur de consistance des covariances
271 - OMA : Observation moins Analysis : Y - Xa
272 - OMB : Observation moins Background : Y - Xb
273 - AMB : Analysis moins Background : Xa - Xb
274 - APosterioriCovariance : matrice A
275 On peut rajouter des variables à stocker dans l'initialisation de
276 l'algorithme élémentaire qui va hériter de cette classe
278 logging.debug("%s Initialisation"%str(name))
279 self._m = PlatformInfo.SystemUsage()
281 self._name = str( name )
282 self._parameters = {}
283 self.__required_parameters = {}
284 self.StoredVariables = {}
286 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
287 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
288 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
289 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
290 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
291 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
292 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
293 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
294 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
295 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
296 self.StoredVariables["ObservedState"] = Persistence.OneVector(name = "ObservedState")
297 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
298 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
299 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
300 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
301 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
302 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
303 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
304 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
305 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
308 logging.debug("%s Lancement"%self._name)
309 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
312 def _post_run(self,_oH=None):
314 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)))
315 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)))
316 logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
317 logging.debug("%s Terminé"%self._name)
320 def get(self, key=None):
322 Renvoie l'une des variables stockées identifiée par la clé, ou le
323 dictionnaire de l'ensemble des variables disponibles en l'absence de
324 clé. Ce sont directement les variables sous forme objet qui sont
325 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
326 des classes de persistance.
329 return self.StoredVariables[key]
331 return self.StoredVariables
333 def has_key(self, key=None):
335 Vérifie si l'une des variables stockées est identifiée par la clé.
337 return self.StoredVariables.has_key(key)
341 Renvoie la liste des clés de variables stockées.
343 return self.StoredVariables.keys()
345 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
347 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
348 sa forme mathématique la plus naturelle possible.
350 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
352 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
354 Permet de définir dans l'algorithme des paramètres requis et leurs
355 caractéristiques par défaut.
358 raise ValueError("A name is mandatory to define a required parameter.")
360 self.__required_parameters[name] = {
362 "typecast" : typecast,
368 logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
370 def getRequiredParameters(self, noDetails=True):
372 Renvoie la liste des noms de paramètres requis ou directement le
373 dictionnaire des paramètres requis.
376 ks = self.__required_parameters.keys()
380 return self.__required_parameters
382 def setParameterValue(self, name=None, value=None):
384 Renvoie la valeur d'un paramètre requis de manière contrôlée
386 default = self.__required_parameters[name]["default"]
387 typecast = self.__required_parameters[name]["typecast"]
388 minval = self.__required_parameters[name]["minval"]
389 maxval = self.__required_parameters[name]["maxval"]
390 listval = self.__required_parameters[name]["listval"]
392 if value is None and default is None:
394 elif value is None and default is not None:
395 if typecast is None: __val = default
396 else: __val = typecast( default )
398 if typecast is None: __val = value
399 else: __val = typecast( value )
401 if minval is not None and (numpy.array(__val) < minval).any():
402 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
403 if maxval is not None and (numpy.array(__val) > maxval).any():
404 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
405 if listval is not None:
406 if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
409 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
410 elif __val not in listval:
411 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
414 def setParameters(self, fromDico={}):
416 Permet de stocker les paramètres reçus dans le dictionnaire interne.
418 self._parameters.update( fromDico )
419 for k in self.__required_parameters.keys():
420 if k in fromDico.keys():
421 self._parameters[k] = self.setParameterValue(k,fromDico[k])
423 self._parameters[k] = self.setParameterValue(k)
424 logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
426 # ==============================================================================
429 Classe générale d'interface de type diagnostic
431 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
432 même temps que l'une des classes de persistance, comme "OneScalar" par
435 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
436 méthode "_formula" pour écrire explicitement et proprement la formule pour
437 l'écriture mathématique du calcul du diagnostic (méthode interne non
438 publique), et "calculate" pour activer la précédente tout en ayant vérifié
439 et préparé les données, et pour stocker les résultats à chaque pas (méthode
440 externe d'activation).
442 def __init__(self, name = "", parameters = {}):
443 self.name = str(name)
444 self.parameters = dict( parameters )
446 def _formula(self, *args):
448 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
449 mathématique la plus naturelle possible.
451 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
453 def calculate(self, *args):
455 Active la formule de calcul avec les arguments correctement rangés
457 raise NotImplementedError("Diagnostic activation method has not been implemented!")
459 # ==============================================================================
462 Classe générale d'interface de type covariance
465 name = "GenericCovariance",
467 asEyeByScalar = None,
468 asEyeByVector = None,
472 Permet de définir une covariance :
473 - asCovariance : entrée des données, comme une matrice compatible avec
474 le constructeur de numpy.matrix
475 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
476 multiplicatif d'une matrice de corrélation identité, aucune matrice
477 n'étant donc explicitement à donner
478 - asEyeByVector : entrée des données comme un seul vecteur de variance,
479 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
480 n'étant donc explicitement à donner
481 - asCovObject : entrée des données comme un objet python, qui a les
482 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
483 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
484 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
486 self.__name = str(name)
489 self.__is_scalar = False
490 self.__is_vector = False
491 self.__is_matrix = False
492 self.__is_object = False
493 if asEyeByScalar is not None:
494 self.__is_scalar = True
495 self.__C = numpy.abs( float(asEyeByScalar) )
498 elif asEyeByVector is not None:
499 self.__is_vector = True
500 self.__C = numpy.abs( numpy.array( numpy.ravel( asEyeByVector ), float ) )
501 self.shape = (self.__C.size,self.__C.size)
502 self.size = self.__C.size**2
503 elif asCovariance is not None:
504 self.__is_matrix = True
505 self.__C = numpy.matrix( asCovariance, float )
506 self.shape = self.__C.shape
507 self.size = self.__C.size
508 elif asCovObject is not None:
509 self.__is_object = True
510 self.__C = asCovObject
511 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
512 if not hasattr(self.__C,at):
513 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
514 if hasattr(self.__C,"shape"):
515 self.shape = self.__C.shape
518 if hasattr(self.__C,"size"):
519 self.size = self.__C.size
524 # 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)
528 def __validate(self):
529 if self.ismatrix() and min(self.shape) != max(self.shape):
530 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))
531 if self.isobject() and min(self.shape) != max(self.shape):
532 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))
533 if self.isscalar() and self.__C <= 0:
534 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
535 if self.isvector() and (self.__C <= 0).any():
536 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
537 if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
539 L = numpy.linalg.cholesky( self.__C )
541 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
544 return self.__is_scalar
547 return self.__is_vector
550 return self.__is_matrix
553 return self.__is_object
557 return Covariance(self.__name+"I", asCovariance = self.__C.I )
558 elif self.isvector():
559 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
560 elif self.isscalar():
561 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
562 elif self.isobject():
563 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
565 return None # Indispensable
569 return Covariance(self.__name+"T", asCovariance = self.__C.T )
570 elif self.isvector():
571 return Covariance(self.__name+"T", asEyeByVector = self.__C )
572 elif self.isscalar():
573 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
574 elif self.isobject():
575 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
579 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
580 elif self.isvector():
581 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
582 elif self.isscalar():
583 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
584 elif self.isobject() and hasattr(self.__C,"cholesky"):
585 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
589 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
590 elif self.isvector():
591 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
592 elif self.isscalar():
593 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
594 elif self.isobject() and hasattr(self.__C,"choleskyI"):
595 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
597 def diag(self, msize=None):
599 return numpy.diag(self.__C)
600 elif self.isvector():
602 elif self.isscalar():
604 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,))
606 return self.__C * numpy.ones(int(msize))
607 elif self.isobject():
608 return self.__C.diag()
610 def asfullmatrix(self, msize=None):
613 elif self.isvector():
614 return numpy.matrix( numpy.diag(self.__C), float )
615 elif self.isscalar():
617 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,))
619 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
620 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
621 return self.__C.asfullmatrix()
623 def trace(self, msize=None):
625 return numpy.trace(self.__C)
626 elif self.isvector():
627 return float(numpy.sum(self.__C))
628 elif self.isscalar():
630 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,))
632 return self.__C * int(msize)
633 elif self.isobject():
634 return self.__C.trace()
637 return repr(self.__C)
642 def __add__(self, other):
643 if self.ismatrix() or self.isobject():
644 return self.__C + numpy.asmatrix(other)
645 elif self.isvector() or self.isscalar():
646 _A = numpy.asarray(other)
647 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
648 return numpy.asmatrix(_A)
650 def __radd__(self, other):
651 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
653 def __sub__(self, other):
654 if self.ismatrix() or self.isobject():
655 return self.__C - numpy.asmatrix(other)
656 elif self.isvector() or self.isscalar():
657 _A = numpy.asarray(other)
658 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
659 return numpy.asmatrix(_A)
661 def __rsub__(self, other):
662 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
667 def __mul__(self, other):
668 if self.ismatrix() and isinstance(other,numpy.matrix):
669 return self.__C * other
670 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
671 or isinstance(other,list) \
672 or isinstance(other,tuple)):
673 if numpy.ravel(other).size == self.shape[1]: # Vecteur
674 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
675 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
676 return self.__C * numpy.asmatrix(other)
678 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
679 elif self.isvector() and (isinstance(other,numpy.matrix) \
680 or isinstance(other,numpy.ndarray) \
681 or isinstance(other,list) \
682 or isinstance(other,tuple)):
683 if numpy.ravel(other).size == self.shape[1]: # Vecteur
684 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
685 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
686 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
688 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
689 elif self.isscalar() and isinstance(other,numpy.matrix):
690 return self.__C * other
691 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
692 or isinstance(other,list) \
693 or isinstance(other,tuple)):
694 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
695 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
697 return self.__C * numpy.asmatrix(other)
698 elif self.isobject():
699 return self.__C.__mul__(other)
701 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
703 def __rmul__(self, other):
704 if self.ismatrix() and isinstance(other,numpy.matrix):
705 return other * self.__C
706 elif self.isvector() and isinstance(other,numpy.matrix):
707 if numpy.ravel(other).size == self.shape[0]: # Vecteur
708 return numpy.asmatrix(numpy.ravel(other) * self.__C)
709 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
710 return numpy.asmatrix(numpy.array(other) * self.__C)
712 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
713 elif self.isscalar() and isinstance(other,numpy.matrix):
714 return other * self.__C
715 elif self.isobject():
716 return self.__C.__rmul__(other)
718 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
723 # ==============================================================================
724 if __name__ == "__main__":
725 print '\n AUTODIAGNOSTIC \n'