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