src/daComposant/daCore/BasicObjects.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2015 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 """
  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         Les caractéristiques de tolérance peuvent être modifées à la création.
  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 Xa
 264             - ObservedState : l'état observé H(X)
 265             - Innovation : l'innovation : d = Y - H(X)
 266             - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
 267             - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
 268             - MahalanobisConsistency : indicateur de consistance des covariances
 269             - OMA : Observation moins Analysis : Y - Xa
 270             - OMB : Observation moins Background : Y - Xb
 271             - AMB : Analysis moins Background : Xa - Xb
 272             - APosterioriCovariance : matrice A
 273         On peut rajouter des variables à stocker dans l'initialisation de
 274         l'algorithme élémentaire qui va hériter de cette classe
 275         """
 276         logging.debug("%s Initialisation"%str(name))
 277         self._m = PlatformInfo.SystemUsage()
 278         #
 279         self._name = str( name )
 280         self._parameters = {}
 281         self.__required_parameters = {}
 282         self.StoredVariables = {}
 283         #
 284         self.StoredVariables["CostFunctionJ"]            = Persistence.OneScalar(name = "CostFunctionJ")
 285         self.StoredVariables["CostFunctionJb"]           = Persistence.OneScalar(name = "CostFunctionJb")
 286         self.StoredVariables["CostFunctionJo"]           = Persistence.OneScalar(name = "CostFunctionJo")
 287         self.StoredVariables["GradientOfCostFunctionJ"]  = Persistence.OneVector(name = "GradientOfCostFunctionJ")
 288         self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
 289         self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
 290         self.StoredVariables["CurrentState"]             = Persistence.OneVector(name = "CurrentState")
 291         self.StoredVariables["Analysis"]                 = Persistence.OneVector(name = "Analysis")
 292         self.StoredVariables["ObservedState"]            = Persistence.OneVector(name = "ObservedState")
 293         self.StoredVariables["Innovation"]               = Persistence.OneVector(name = "Innovation")
 294         self.StoredVariables["SigmaObs2"]                = Persistence.OneScalar(name = "SigmaObs2")
 295         self.StoredVariables["SigmaBck2"]                = Persistence.OneScalar(name = "SigmaBck2")
 296         self.StoredVariables["MahalanobisConsistency"]   = Persistence.OneScalar(name = "MahalanobisConsistency")
 297         self.StoredVariables["OMA"]                      = Persistence.OneVector(name = "OMA")
 298         self.StoredVariables["OMB"]                      = Persistence.OneVector(name = "OMB")
 299         self.StoredVariables["BMA"]                      = Persistence.OneVector(name = "BMA")
 300         self.StoredVariables["APosterioriCovariance"]    = Persistence.OneMatrix(name = "APosterioriCovariance")
 301         self.StoredVariables["SimulationQuantiles"]      = Persistence.OneMatrix(name = "SimulationQuantiles")
 302
 303     def _pre_run(self):
 304         logging.debug("%s Lancement"%self._name)
 305         logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
 306         return 0
 307
 308     def _post_run(self,_oH=None):
 309         if _oH is not None:
 310             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)))
 311             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)))
 312         logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
 313         logging.debug("%s Terminé"%self._name)
 314         return 0
 315
 316     def get(self, key=None):
 317         """
 318         Renvoie l'une des variables stockées identifiée par la clé, ou le
 319         dictionnaire de l'ensemble des variables disponibles en l'absence de
 320         clé. Ce sont directement les variables sous forme objet qui sont
 321         renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
 322         des classes de persistance.
 323         """
 324         if key is not None:
 325             return self.StoredVariables[key]
 326         else:
 327             return self.StoredVariables
 328
 329     def has_key(self, key=None):
 330         """
 331         Vérifie si l'une des variables stockées est identifiée par la clé.
 332         """
 333         return self.StoredVariables.has_key(key)
 334
 335     def keys(self):
 336         """
 337         Renvoie la liste des clés de variables stockées.
 338         """
 339         return self.StoredVariables.keys()
 340
 341     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
 342         """
 343         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
 344         sa forme mathématique la plus naturelle possible.
 345         """
 346         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
 347
 348     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
 349         """
 350         Permet de définir dans l'algorithme des paramètres requis et leurs
 351         caractéristiques par défaut.
 352         """
 353         if name is None:
 354             raise ValueError("A name is mandatory to define a required parameter.")
 355         #
 356         self.__required_parameters[name] = {
 357             "default"  : default,
 358             "typecast" : typecast,
 359             "minval"   : minval,
 360             "maxval"   : maxval,
 361             "listval"  : listval,
 362             "message"  : message,
 363             }
 364         logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
 365
 366     def getRequiredParameters(self, noDetails=True):
 367         """
 368         Renvoie la liste des noms de paramètres requis ou directement le
 369         dictionnaire des paramètres requis.
 370         """
 371         if noDetails:
 372             ks = self.__required_parameters.keys()
 373             ks.sort()
 374             return ks
 375         else:
 376             return self.__required_parameters
 377
 378     def setParameterValue(self, name=None, value=None):
 379         """
 380         Renvoie la valeur d'un paramètre requis de manière contrôlée
 381         """
 382         default  = self.__required_parameters[name]["default"]
 383         typecast = self.__required_parameters[name]["typecast"]
 384         minval   = self.__required_parameters[name]["minval"]
 385         maxval   = self.__required_parameters[name]["maxval"]
 386         listval  = self.__required_parameters[name]["listval"]
 387         #
 388         if value is None and default is None:
 389             __val = None
 390         elif value is None and default is not None:
 391             if typecast is None: __val = default
 392             else:                __val = typecast( default )
 393         else:
 394             if typecast is None: __val = value
 395             else:                __val = typecast( value )
 396         #
 397         if minval is not None and (numpy.array(__val) < minval).any():
 398             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
 399         if maxval is not None and (numpy.array(__val) > maxval).any():
 400             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
 401         if listval is not None:
 402             if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
 403                 for v in __val:
 404                     if v 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."%(v, name, listval))
 406             elif __val not in listval:
 407                 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
 408         return __val
 409
 410     def setParameters(self, fromDico={}):
 411         """
 412         Permet de stocker les paramètres reçus dans le dictionnaire interne.
 413         """
 414         self._parameters.update( fromDico )
 415         for k in self.__required_parameters.keys():
 416             if k in fromDico.keys():
 417                 self._parameters[k] = self.setParameterValue(k,fromDico[k])
 418             else:
 419                 self._parameters[k] = self.setParameterValue(k)
 420             logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
 421
 422 # ==============================================================================
 423 class Diagnostic:
 424     """
 425     Classe générale d'interface de type diagnostic
 426
 427     Ce template s'utilise de la manière suivante : il sert de classe "patron" en
 428     même temps que l'une des classes de persistance, comme "OneScalar" par
 429     exemple.
 430
 431     Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
 432     méthode "_formula" pour écrire explicitement et proprement la formule pour
 433     l'écriture mathématique du calcul du diagnostic (méthode interne non
 434     publique), et "calculate" pour activer la précédente tout en ayant vérifié
 435     et préparé les données, et pour stocker les résultats à chaque pas (méthode
 436     externe d'activation).
 437     """
 438     def __init__(self, name = "", parameters = {}):
 439         self.name       = str(name)
 440         self.parameters = dict( parameters )
 441
 442     def _formula(self, *args):
 443         """
 444         Doit implémenter l'opération élémentaire de diagnostic sous sa forme
 445         mathématique la plus naturelle possible.
 446         """
 447         raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
 448
 449     def calculate(self, *args):
 450         """
 451         Active la formule de calcul avec les arguments correctement rangés
 452         """
 453         raise NotImplementedError("Diagnostic activation method has not been implemented!")
 454
 455 # ==============================================================================
 456 class Covariance:
 457     """
 458     Classe générale d'interface de type covariance
 459     """
 460     def __init__(self,
 461             name          = "GenericCovariance",
 462             asCovariance  = None,
 463             asEyeByScalar = None,
 464             asEyeByVector = None,
 465             ):
 466         """
 467         Permet de définir une covariance :
 468         - asCovariance : entrée des données, comme une matrice compatible avec
 469           le constructeur de numpy.matrix
 470         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 471           multiplicatif d'une matrice de corrélation identité, aucune matrice
 472           n'étant donc explicitement à donner
 473         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 474           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 475           n'étant donc explicitement à donner
 476         """
 477         self.__name       = str(name)
 478         #
 479         self.__B          = None
 480         self.__is_scalar  = False
 481         self.__is_vector  = False
 482         self.__is_matrix  = False
 483         if asEyeByScalar is not None:
 484             self.__is_scalar = True
 485             self.__B         = numpy.abs( float(asEyeByScalar) )
 486             self.shape       = (0,0)
 487             self.size        = 0
 488         elif asEyeByVector is not None:
 489             self.__is_vector = True
 490             self.__B         = numpy.abs( numpy.array( numpy.ravel( asEyeByVector ), float ) )
 491             self.shape       = (self.__B.size,self.__B.size)
 492             self.size        = self.__B.size**2
 493         elif asCovariance is not None:
 494             self.__is_matrix = True
 495             self.__B         = numpy.matrix( asCovariance, float )
 496             self.shape       = self.__B.shape
 497             self.size        = self.__B.size
 498         else:
 499             pass
 500             # 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)
 501         #
 502         self.__validate()
 503
 504     def __validate(self):
 505         if self.ismatrix() and min(self.shape) != max(self.shape):
 506             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))
 507         if self.isscalar() and self.__B <= 0:
 508             raise ValueError("The %s covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__B_scalar))
 509         if self.isvector() and (self.__B <= 0).any():
 510             raise ValueError("The %s covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
 511         if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
 512             try:
 513                 L = numpy.linalg.cholesky( self.__B )
 514             except:
 515                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 516
 517     def isscalar(self):
 518         return self.__is_scalar
 519
 520     def isvector(self):
 521         return self.__is_vector
 522
 523     def ismatrix(self):
 524         return self.__is_matrix
 525
 526     def getI(self):
 527         if   self.ismatrix():
 528             return Covariance(self.__name+"I", asCovariance  = self.__B.I )
 529         elif self.isvector():
 530             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__B )
 531         elif self.isscalar():
 532             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__B )
 533         else:
 534             return None
 535
 536     def getT(self):
 537         if   self.ismatrix():
 538             return Covariance(self.__name+"T", asCovariance  = self.__B.T )
 539         elif self.isvector():
 540             return Covariance(self.__name+"T", asEyeByVector = self.__B )
 541         elif self.isscalar():
 542             return Covariance(self.__name+"T", asEyeByScalar = self.__B )
 543
 544     def cholesky(self):
 545         if   self.ismatrix():
 546             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__B) )
 547         elif self.isvector():
 548             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__B ) )
 549         elif self.isscalar():
 550             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__B ) )
 551
 552     def choleskyI(self):
 553         if   self.ismatrix():
 554             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__B).I )
 555         elif self.isvector():
 556             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__B ) )
 557         elif self.isscalar():
 558             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__B ) )
 559
 560     def diag(self, msize=None):
 561         if   self.ismatrix():
 562             return numpy.diag(self.__B)
 563         elif self.isvector():
 564             return self.__B
 565         elif self.isscalar():
 566             if msize is None:
 567                 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,))
 568             else:
 569                 return self.__B * numpy.ones(int(msize))
 570
 571     def asfullmatrix(self, msize=None):
 572         if   self.ismatrix():
 573             return self.__B
 574         elif self.isvector():
 575             return numpy.matrix( numpy.diag(self.__B), float )
 576         elif self.isscalar():
 577             if msize is None:
 578                 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,))
 579             else:
 580                 return numpy.matrix( self.__B * numpy.eye(int(msize)), float )
 581
 582     def trace(self, msize=None):
 583         if   self.ismatrix():
 584             return numpy.trace(self.__B)
 585         elif self.isvector():
 586             return float(numpy.sum(self.__B))
 587         elif self.isscalar():
 588             if msize is None:
 589                 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,))
 590             else:
 591                 return self.__B * int(msize)
 592
 593     def __repr__(self):
 594         return repr(self.__B)
 595
 596     def __str__(self):
 597         return str(self.__B)
 598
 599     def __add__(self, other):
 600         if   self.ismatrix():
 601             return self.__B + numpy.asmatrix(other)
 602         elif self.isvector() or self.isscalar():
 603             _A = numpy.asarray(other)
 604             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__B
 605             return numpy.asmatrix(_A)
 606
 607     def __radd__(self, other):
 608         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
 609
 610     def __sub__(self, other):
 611         if   self.ismatrix():
 612             return self.__B - numpy.asmatrix(other)
 613         elif self.isvector() or self.isscalar():
 614             _A = numpy.asarray(other)
 615             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__B - _A.reshape(_A.size)[::_A.shape[1]+1]
 616             return numpy.asmatrix(_A)
 617
 618     def __rsub__(self, other):
 619         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
 620
 621     def __neg__(self):
 622         return - self.__B
 623
 624     def __mul__(self, other):
 625         if   self.ismatrix() and isinstance(other,numpy.matrix):
 626             return self.__B * other
 627         elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
 628                                or isinstance(other,list) \
 629                                or isinstance(other,tuple)):
 630             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 631                 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
 632             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 633                 return self.__B * numpy.asmatrix(other)
 634             else:
 635                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
 636         elif self.isvector() and (isinstance(other,numpy.matrix) \
 637                                or isinstance(other,numpy.ndarray) \
 638                                or isinstance(other,list) \
 639                                or isinstance(other,tuple)):
 640             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 641                 return numpy.asmatrix(self.__B * numpy.ravel(other)).T
 642             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 643                 return numpy.asmatrix((self.__B * (numpy.asarray(other).transpose())).transpose())
 644             else:
 645                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 646         elif self.isscalar() and isinstance(other,numpy.matrix):
 647             return self.__B * other
 648         elif self.isscalar() and (isinstance(other,numpy.ndarray) \
 649                                or isinstance(other,list) \
 650                                or isinstance(other,tuple)):
 651             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
 652                 return self.__B * numpy.asmatrix(numpy.ravel(other)).T
 653             else:
 654                 return self.__B * numpy.asmatrix(other)
 655         else:
 656             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
 657
 658     def __rmul__(self, other):
 659         if self.ismatrix() and isinstance(other,numpy.matrix):
 660             return other * self.__B
 661         elif self.isvector() and isinstance(other,numpy.matrix):
 662             if numpy.ravel(other).size == self.shape[0]: # Vecteur
 663                 return numpy.asmatrix(numpy.ravel(other) * self.__B)
 664             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
 665                 return numpy.asmatrix(numpy.array(other) * self.__B)
 666             else:
 667                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 668         elif self.isscalar() and isinstance(other,numpy.matrix):
 669             return other * self.__B
 670         else:
 671             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
 672
 673     def __len__(self):
 674         return self.shape[0]
 675
 676 # ==============================================================================
 677 if __name__ == "__main__":
 678     print '\n AUTODIAGNOSTIC \n'