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