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