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