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