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