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