src/daComposant/daCore/BasicObjects.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2014 EDF R&D
   4 #
   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.
   9 #
  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.
  14 #
  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
  18 #
  19 #  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 #  Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 __doc__ = """
  24     Définit les outils généraux élémentaires.
  25
  26     Ce module est destiné à etre appelée par AssimilationStudy pour constituer
  27     les objets élémentaires de l'algorithme.
  28 """
  29 __author__ = "Jean-Philippe ARGAUD"
  30
  31 import logging, copy
  32 import numpy
  33 import Persistence
  34
  35 # ==============================================================================
  36 class CacheManager:
  37     """
  38     Classe générale de gestion d'un cache de calculs
  39     """
  40     def __init__(self,
  41             toleranceInRedundancy = 1.e-18,
  42             lenghtOfRedundancy    = -1,
  43             ):
  44         """
  45
  46         """
  47         self.__tolerBP  = float(toleranceInRedundancy)
  48         self.__lenghtOR = int(lenghtOfRedundancy)
  49         self.clearCache()
  50
  51     def clearCache(self):
  52         self.__listOPCP = [] # Operator Previous Calculated Points
  53         self.__listOPCR = [] # Operator Previous Calculated Results
  54         self.__listOPCN = [] # Operator Previous Calculated Point Norms
  55         self.__ac       = False
  56         self.__iac      = -1
  57
  58     def wasCalculatedIn(self, xValue, info="" ):
  59         self.__ac, self.__iac = False, -1
  60         for i in xrange(len(self.__listOPCP)-1,-1,-1):
  61             if xValue.size != self.__listOPCP[i].size:
  62                 continue
  63             if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCP[i]) < self.__tolerBP * self.__listOPCN[i]:
  64                 self.__ac, self.__iac = True, i
  65                 break
  66         return self.__ac
  67
  68     def getValueInX(self, info="" ):
  69         if self.__ac and (-1 < self.__iac < len(self.__listOPCR)):
  70             __HX = self.__listOPCR[self.__iac]
  71             self.__ac, self.__iac = False, -1
  72             return __HX
  73         else:
  74             raise ValueError("CM Cas%s non encore disponible, vérifier au préalable son existence"%info)
  75
  76     def storeValueInX(self, xValue, HxValue ):
  77         if self.__lenghtOR < 0: self.__lenghtOR = 2 * xValue.size + 2
  78         if len(self.__listOPCP) > self.__lenghtOR:
  79             self.__listOPCP.pop(0)
  80             self.__listOPCR.pop(0)
  81             self.__listOPCN.pop(0)
  82         self.__listOPCP.append( copy.copy(numpy.ravel(xValue)) )
  83         self.__listOPCR.append( copy.copy(HxValue) )
  84         self.__listOPCN.append( numpy.linalg.norm(xValue) )
  85
  86 # ==============================================================================
  87 class Operator:
  88     """
  89     Classe générale d'interface de type opérateur
  90     """
  91     NbCallsAsMatrix = 0
  92     NbCallsAsMethod = 0
  93     NbCallsOfCached = 0
  94     CM = CacheManager()
  95     #
  96     def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
  97         """
  98         On construit un objet de ce type en fournissant à l'aide de l'un des
  99         deux mots-clé, soit une fonction python, soit matrice.
 100         Arguments :
 101         - fromMethod : argument de type fonction Python
 102         - fromMatrix : argument adapté au constructeur numpy.matrix
 103         """
 104         self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
 105         self.__AvoidRC = bool( avoidingRedundancy )
 106         if   fromMethod is not None:
 107             self.__Method = fromMethod
 108             self.__Matrix = None
 109             self.__Type   = "Method"
 110         elif fromMatrix is not None:
 111             self.__Method = None
 112             self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
 113             self.__Type   = "Matrix"
 114         else:
 115             self.__Method = None
 116             self.__Matrix = None
 117             self.__Type   = None
 118
 119     def isType(self):
 120         return self.__Type
 121
 122     def appliedTo(self, xValue):
 123         """
 124         Permet de restituer le résultat de l'application de l'opérateur à un
 125         argument xValue. Cette méthode se contente d'appliquer, son argument
 126         devant a priori être du bon type.
 127         Arguments :
 128         - xValue : argument adapté pour appliquer l'opérateur
 129         """
 130         if self.__AvoidRC and Operator.CM.wasCalculatedIn(xValue):
 131             __alreadyCalculated = True
 132         else:
 133             __alreadyCalculated = False
 134         #
 135         if __alreadyCalculated:
 136             self.__addOneCacheCall()
 137             HxValue = Operator.CM.getValueInX()
 138         else:
 139             if self.__Matrix is not None:
 140                 self.__addOneMatrixCall()
 141                 HxValue = self.__Matrix * xValue
 142             else:
 143                 self.__addOneMethodCall()
 144                 HxValue = self.__Method( xValue )
 145             if self.__AvoidRC:
 146                 Operator.CM.storeValueInX(xValue,HxValue)
 147         #
 148         return HxValue
 149
 150     def appliedControledFormTo(self, (xValue, uValue) ):
 151         """
 152         Permet de restituer le résultat de l'application de l'opérateur à une
 153         paire (xValue, uValue). Cette méthode se contente d'appliquer, son
 154         argument devant a priori être du bon type. Si la uValue est None,
 155         on suppose que l'opérateur ne s'applique qu'à xValue.
 156         Arguments :
 157         - xValue : argument X adapté pour appliquer l'opérateur
 158         - uValue : argument U adapté pour appliquer l'opérateur
 159         """
 160         if self.__Matrix is not None:
 161             self.__addOneMatrixCall()
 162             return self.__Matrix * xValue
 163         elif uValue is not None:
 164             self.__addOneMethodCall()
 165             return self.__Method( (xValue, uValue) )
 166         else:
 167             self.__addOneMethodCall()
 168             return self.__Method( xValue )
 169
 170     def appliedInXTo(self, (xNominal, xValue) ):
 171         """
 172         Permet de restituer le résultat de l'application de l'opérateur à un
 173         argument xValue, sachant que l'opérateur est valable en xNominal.
 174         Cette méthode se contente d'appliquer, son argument devant a priori
 175         être du bon type. Si l'opérateur est linéaire car c'est une matrice,
 176         alors il est valable en tout point nominal et il n'est pas nécessaire
 177         d'utiliser xNominal.
 178         Arguments : une liste contenant
 179         - xNominal : argument permettant de donner le point où l'opérateur
 180           est construit pour etre ensuite appliqué
 181         - xValue : argument adapté pour appliquer l'opérateur
 182         """
 183         if self.__Matrix is not None:
 184             self.__addOneMatrixCall()
 185             return self.__Matrix * xValue
 186         else:
 187             self.__addOneMethodCall()
 188             return self.__Method( (xNominal, xValue) )
 189
 190     def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
 191         """
 192         Permet de renvoyer l'opérateur sous la forme d'une matrice
 193         """
 194         if self.__Matrix is not None:
 195             self.__addOneMatrixCall()
 196             return self.__Matrix
 197         elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
 198             self.__addOneMethodCall()
 199             return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
 200         else:
 201             raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
 202
 203     def shape(self):
 204         """
 205         Renvoie la taille sous forme numpy si l'opérateur est disponible sous
 206         la forme d'une matrice
 207         """
 208         if self.__Matrix is not None:
 209             return self.__Matrix.shape
 210         else:
 211             raise ValueError("Matrix form of the operator is not available, nor the shape")
 212
 213     def nbcalls(self, which=None):
 214         """
 215         Renvoie les nombres d'évaluations de l'opérateur
 216         """
 217         __nbcalls = (
 218             self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
 219             self.__NbCallsAsMatrix,
 220             self.__NbCallsAsMethod,
 221             self.__NbCallsOfCached,
 222             Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
 223             Operator.NbCallsAsMatrix,
 224             Operator.NbCallsAsMethod,
 225             Operator.NbCallsOfCached,
 226             )
 227         if which is None: return __nbcalls
 228         else:             return __nbcalls[which]
 229
 230     def __addOneMatrixCall(self):
 231         self.__NbCallsAsMatrix   += 1 # Decompte local
 232         Operator.NbCallsAsMatrix += 1 # Decompte global
 233
 234     def __addOneMethodCall(self):
 235         self.__NbCallsAsMethod   += 1 # Decompte local
 236         Operator.NbCallsAsMethod += 1 # Decompte global
 237
 238     def __addOneCacheCall(self):
 239         self.__NbCallsOfCached   += 1 # Decompte local
 240         Operator.NbCallsOfCached += 1 # Decompte global
 241
 242 # ==============================================================================
 243 class Algorithm:
 244     """
 245     Classe générale d'interface de type algorithme
 246
 247     Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
 248     d'assimilation, en fournissant un container (dictionnaire) de variables
 249     persistantes initialisées, et des méthodes d'accès à ces variables stockées.
 250
 251     Une classe élémentaire d'algorithme doit implémenter la méthode "run".
 252     """
 253     def __init__(self, name):
 254         """
 255         L'initialisation présente permet de fabriquer des variables de stockage
 256         disponibles de manière générique dans les algorithmes élémentaires. Ces
 257         variables de stockage sont ensuite conservées dans un dictionnaire
 258         interne à l'objet, mais auquel on accède par la méthode "get".
 259
 260         Les variables prévues sont :
 261             - CostFunctionJ  : fonction-cout globale, somme des deux parties suivantes
 262             - CostFunctionJb : partie ébauche ou background de la fonction-cout
 263             - CostFunctionJo : partie observations de la fonction-cout
 264             - GradientOfCostFunctionJ  : gradient de la fonction-cout globale
 265             - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
 266             - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
 267             - CurrentState : état courant lors d'itérations
 268             - Analysis : l'analyse
 269             - Innovation : l'innovation : d = Y - H Xb
 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         On peut rajouter des variables à stocker dans l'initialisation de
 278         l'algorithme élémentaire qui va hériter de cette classe
 279         """
 280         logging.debug("%s Initialisation"%str(name))
 281         self._name = str( name )
 282         self._parameters = {}
 283         self.__required_parameters = {}
 284         self.StoredVariables = {}
 285         #
 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["Innovation"]               = Persistence.OneVector(name = "Innovation")
 295         self.StoredVariables["SigmaObs2"]                = Persistence.OneScalar(name = "SigmaObs2")
 296         self.StoredVariables["SigmaBck2"]                = Persistence.OneScalar(name = "SigmaBck2")
 297         self.StoredVariables["MahalanobisConsistency"]   = Persistence.OneScalar(name = "MahalanobisConsistency")
 298         self.StoredVariables["OMA"]                      = Persistence.OneVector(name = "OMA")
 299         self.StoredVariables["OMB"]                      = Persistence.OneVector(name = "OMB")
 300         self.StoredVariables["BMA"]                      = Persistence.OneVector(name = "BMA")
 301         self.StoredVariables["APosterioriCovariance"]    = Persistence.OneMatrix(name = "APosterioriCovariance")
 302         self.StoredVariables["SimulationQuantiles"]      = Persistence.OneMatrix(name = "SimulationQuantiles")
 303
 304     def get(self, key=None):
 305         """
 306         Renvoie l'une des variables stockées identifiée par la clé, ou le
 307         dictionnaire de l'ensemble des variables disponibles en l'absence de
 308         clé. Ce sont directement les variables sous forme objet qui sont
 309         renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
 310         des classes de persistance.
 311         """
 312         if key is not None:
 313             return self.StoredVariables[key]
 314         else:
 315             return self.StoredVariables
 316
 317     def has_key(self, key=None):
 318         """
 319         Vérifie si l'une des variables stockées est identifiée par la clé.
 320         """
 321         return self.StoredVariables.has_key(key)
 322
 323     def keys(self):
 324         """
 325         Renvoie la liste des clés de variables stockées.
 326         """
 327         return self.StoredVariables.keys()
 328
 329     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
 330         """
 331         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
 332         sa forme mathématique la plus naturelle possible.
 333         """
 334         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
 335
 336     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
 337         """
 338         Permet de définir dans l'algorithme des paramètres requis et leurs
 339         caractéristiques par défaut.
 340         """
 341         if name is None:
 342             raise ValueError("A name is mandatory to define a required parameter.")
 343         #
 344         self.__required_parameters[name] = {
 345             "default"  : default,
 346             "typecast" : typecast,
 347             "minval"   : minval,
 348             "maxval"   : maxval,
 349             "listval"  : listval,
 350             "message"  : message,
 351             }
 352         logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
 353
 354     def getRequiredParameters(self, noDetails=True):
 355         """
 356         Renvoie la liste des noms de paramètres requis ou directement le
 357         dictionnaire des paramètres requis.
 358         """
 359         if noDetails:
 360             ks = self.__required_parameters.keys()
 361             ks.sort()
 362             return ks
 363         else:
 364             return self.__required_parameters
 365
 366     def setParameterValue(self, name=None, value=None):
 367         """
 368         Renvoie la valeur d'un paramètre requis de manière contrôlée
 369         """
 370         default  = self.__required_parameters[name]["default"]
 371         typecast = self.__required_parameters[name]["typecast"]
 372         minval   = self.__required_parameters[name]["minval"]
 373         maxval   = self.__required_parameters[name]["maxval"]
 374         listval  = self.__required_parameters[name]["listval"]
 375         #
 376         if value is None and default is None:
 377             __val = None
 378         elif value is None and default is not None:
 379             if typecast is None: __val = default
 380             else:                __val = typecast( default )
 381         else:
 382             if typecast is None: __val = value
 383             else:                __val = typecast( value )
 384         #
 385         if minval is not None and __val < minval:
 386             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
 387         if maxval is not None and __val > maxval:
 388             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
 389         if listval is not None:
 390             if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
 391                 for v in __val:
 392                     if v not in listval:
 393                         raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
 394             elif __val not in listval:
 395                 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
 396         return __val
 397
 398     def setParameters(self, fromDico={}):
 399         """
 400         Permet de stocker les paramètres reçus dans le dictionnaire interne.
 401         """
 402         self._parameters.update( fromDico )
 403         for k in self.__required_parameters.keys():
 404             if k in fromDico.keys():
 405                 self._parameters[k] = self.setParameterValue(k,fromDico[k])
 406             else:
 407                 self._parameters[k] = self.setParameterValue(k)
 408             logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
 409
 410 # ==============================================================================
 411 class Diagnostic:
 412     """
 413     Classe générale d'interface de type diagnostic
 414
 415     Ce template s'utilise de la manière suivante : il sert de classe "patron" en
 416     même temps que l'une des classes de persistance, comme "OneScalar" par
 417     exemple.
 418
 419     Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
 420     méthode "_formula" pour écrire explicitement et proprement la formule pour
 421     l'écriture mathématique du calcul du diagnostic (méthode interne non
 422     publique), et "calculate" pour activer la précédente tout en ayant vérifié
 423     et préparé les données, et pour stocker les résultats à chaque pas (méthode
 424     externe d'activation).
 425     """
 426     def __init__(self, name = "", parameters = {}):
 427         self.name       = str(name)
 428         self.parameters = dict( parameters )
 429
 430     def _formula(self, *args):
 431         """
 432         Doit implémenter l'opération élémentaire de diagnostic sous sa forme
 433         mathématique la plus naturelle possible.
 434         """
 435         raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
 436
 437     def calculate(self, *args):
 438         """
 439         Active la formule de calcul avec les arguments correctement rangés
 440         """
 441         raise NotImplementedError("Diagnostic activation method has not been implemented!")
 442
 443 # ==============================================================================
 444 class Covariance:
 445     """
 446     Classe générale d'interface de type covariance
 447     """
 448     def __init__(self,
 449             name          = "GenericCovariance",
 450             asCovariance  = None,
 451             asEyeByScalar = None,
 452             asEyeByVector = None,
 453             ):
 454         """
 455         Permet de définir une covariance :
 456         - asCovariance : entrée des données, comme une matrice compatible avec
 457           le constructeur de numpy.matrix
 458         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 459           multiplicatif d'une matrice de corrélation identité, aucune matrice
 460           n'étant donc explicitement à donner
 461         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 462           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 463           n'étant donc explicitement à donner
 464         """
 465         self.__name       = str(name)
 466         #
 467         self.__B          = None
 468         self.__is_scalar  = False
 469         self.__is_vector  = False
 470         self.__is_matrix  = False
 471         if asEyeByScalar is not None:
 472             self.__is_scalar = True
 473             self.__B         = numpy.abs( float(asEyeByScalar) )
 474             self.shape       = (0,0)
 475             self.size        = 0
 476         elif asEyeByVector is not None:
 477             self.__is_vector = True
 478             self.__B         = numpy.abs( numpy.array( numpy.ravel( asEyeByVector ), float ) )
 479             self.shape       = (self.__B.size,self.__B.size)
 480             self.size        = self.__B.size**2
 481         elif asCovariance is not None:
 482             self.__is_matrix = True
 483             self.__B         = numpy.matrix( asCovariance, float )
 484             self.shape       = self.__B.shape
 485             self.size        = self.__B.size
 486         else:
 487             pass
 488             # 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)
 489         #
 490         self.__validate()
 491
 492     def __validate(self):
 493         if self.ismatrix() and min(self.shape) != max(self.shape):
 494             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))
 495         if self.isscalar() and self.__B <= 0:
 496             raise ValueError("The %s covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__B_scalar))
 497         if self.isvector() and (self.__B <= 0).any():
 498             raise ValueError("The %s covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
 499         if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
 500             try:
 501                 L = numpy.linalg.cholesky( self.__B )
 502             except:
 503                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 504
 505     def isscalar(self):
 506         return self.__is_scalar
 507
 508     def isvector(self):
 509         return self.__is_vector
 510
 511     def ismatrix(self):
 512         return self.__is_matrix
 513
 514     def getI(self):
 515         if   self.ismatrix():
 516             return Covariance(self.__name+"I", asCovariance  = self.__B.I )
 517         elif self.isvector():
 518             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__B )
 519         elif self.isscalar():
 520             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__B )
 521         else:
 522             return None
 523
 524     def getT(self):
 525         if   self.ismatrix():
 526             return Covariance(self.__name+"T", asCovariance  = self.__B.T )
 527         elif self.isvector():
 528             return Covariance(self.__name+"T", asEyeByVector = self.__B )
 529         elif self.isscalar():
 530             return Covariance(self.__name+"T", asEyeByScalar = self.__B )
 531
 532     def cholesky(self):
 533         if   self.ismatrix():
 534             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__B) )
 535         elif self.isvector():
 536             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__B ) )
 537         elif self.isscalar():
 538             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__B ) )
 539
 540     def choleskyI(self):
 541         if   self.ismatrix():
 542             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__B).I )
 543         elif self.isvector():
 544             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__B ) )
 545         elif self.isscalar():
 546             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__B ) )
 547
 548     def diag(self, msize=None):
 549         if   self.ismatrix():
 550             return numpy.diag(self.__B)
 551         elif self.isvector():
 552             return self.__B
 553         elif self.isscalar():
 554             if msize is None:
 555                 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,))
 556             else:
 557                 return self.__B * numpy.ones(int(msize))
 558
 559     def asfullmatrix(self, msize=None):
 560         if   self.ismatrix():
 561             return self.__B
 562         elif self.isvector():
 563             return numpy.matrix( numpy.diag(self.__B), float )
 564         elif self.isscalar():
 565             if msize is None:
 566                 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,))
 567             else:
 568                 return numpy.matrix( self.__B * numpy.eye(int(msize)), float )
 569
 570     def trace(self, msize=None):
 571         if   self.ismatrix():
 572             return numpy.trace(self.__B)
 573         elif self.isvector():
 574             return float(numpy.sum(self.__B))
 575         elif self.isscalar():
 576             if msize is None:
 577                 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,))
 578             else:
 579                 return self.__B * int(msize)
 580
 581     def __repr__(self):
 582         return repr(self.__B)
 583
 584     def __str__(self):
 585         return str(self.__B)
 586
 587     def __add__(self, other):
 588         if   self.ismatrix():
 589             return self.__B + numpy.asmatrix(other)
 590         elif self.isvector() or self.isscalar():
 591             _A = numpy.asarray(other)
 592             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__B
 593             return numpy.asmatrix(_A)
 594
 595     def __radd__(self, other):
 596         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
 597
 598     def __sub__(self, other):
 599         if   self.ismatrix():
 600             return self.__B - numpy.asmatrix(other)
 601         elif self.isvector() or self.isscalar():
 602             _A = numpy.asarray(other)
 603             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__B - _A.reshape(_A.size)[::_A.shape[1]+1]
 604             return numpy.asmatrix(_A)
 605
 606     def __rsub__(self, other):
 607         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
 608
 609     def __neg__(self):
 610         return - self.__B
 611
 612     def __mul__(self, other):
 613         if   self.ismatrix() and isinstance(other,numpy.matrix):
 614             return self.__B * other
 615         elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
 616                                or isinstance(other,list) \
 617                                or isinstance(other,tuple)):
 618             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 619                 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
 620             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 621                 return self.__B * numpy.asmatrix(other)
 622             else:
 623                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
 624         elif self.isvector() and (isinstance(other,numpy.matrix) \
 625                                or isinstance(other,numpy.ndarray) \
 626                                or isinstance(other,list) \
 627                                or isinstance(other,tuple)):
 628             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 629                 return numpy.asmatrix(self.__B * numpy.ravel(other)).T
 630             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 631                 return numpy.asmatrix((self.__B * (numpy.asarray(other).transpose())).transpose())
 632             else:
 633                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 634         elif self.isscalar() and isinstance(other,numpy.matrix):
 635             return self.__B * other
 636         elif self.isscalar() and (isinstance(other,numpy.ndarray) \
 637                                or isinstance(other,list) \
 638                                or isinstance(other,tuple)):
 639             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
 640                 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
 641             else:
 642                 return self.__B * numpy.asmatrix(other)
 643         else:
 644             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
 645
 646     def __rmul__(self, other):
 647         if self.ismatrix() and isinstance(other,numpy.matrix):
 648             return other * self.__B
 649         elif self.isvector() and isinstance(other,numpy.matrix):
 650             if numpy.ravel(other).size == self.shape[0]: # Vecteur
 651                 return numpy.asmatrix(numpy.ravel(other) * self.__B)
 652             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
 653                 return numpy.asmatrix(numpy.array(other) * self.__B)
 654             else:
 655                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 656         elif self.isscalar() and isinstance(other,numpy.matrix):
 657             return other * self.__B
 658         else:
 659             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
 660
 661     def __len__(self):
 662         return self.shape[0]
 663
 664 # ==============================================================================
 665 if __name__ == "__main__":
 666     print '\n AUTODIAGNOSTIC \n'