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