src/daComposant/daCore/BasicObjects.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2015 EDF R&D
   4 #
   5 #  This library is free software; you can redistribute it and/or
   6 #  modify it under the terms of the GNU Lesser General Public
   7 #  License as published by the Free Software Foundation; either
   8 #  version 2.1 of the License.
   9 #
  10 #  This library is distributed in the hope that it will be useful,
  11 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 #  Lesser General Public License for more details.
  14 #
  15 #  You should have received a copy of the GNU Lesser General Public
  16 #  License along with this library; if not, write to the Free Software
  17 #  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 #  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 #  Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 """
  24     Définit les outils généraux élémentaires.
  25
  26     Ce module est destiné à etre appelée par AssimilationStudy pour constituer
  27     les objets élémentaires de l'algorithme.
  28 """
  29 __author__ = "Jean-Philippe ARGAUD"
  30 __all__ = []
  31
  32 import logging, copy
  33 import numpy
  34 import Persistence
  35 import PlatformInfo
  36
  37 # ==============================================================================
  38 class CacheManager:
  39     """
  40     Classe générale de gestion d'un cache de calculs
  41     """
  42     def __init__(self,
  43             toleranceInRedundancy = 1.e-18,
  44             lenghtOfRedundancy    = -1,
  45             ):
  46         """
  47         Les caractéristiques de tolérance peuvent être modifées à la création.
  48         """
  49         self.__tolerBP  = float(toleranceInRedundancy)
  50         self.__lenghtOR = int(lenghtOfRedundancy)
  51         self.__initlnOR = self.__lenghtOR
  52         self.clearCache()
  53
  54     def clearCache(self):
  55         self.__listOPCV = [] # Operator Previous Calculated Points, Results, Point Norms
  56         # logging.debug("CM Tolerance de determination des doublons : %.2e"%self.__tolerBP)
  57
  58     def wasCalculatedIn(self, xValue ): #, info="" ):
  59         __alc = False
  60         __HxV = None
  61         for i in xrange(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
  62             if xValue.size != self.__listOPCV[i][0].size:
  63                 # logging.debug("CM Différence de la taille %s de X et de celle %s du point %i déjà calculé"%(xValue.shape,i,self.__listOPCP[i].shape))
  64                 continue
  65             if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
  66                 __alc  = True
  67                 __HxV = self.__listOPCV[i][1]
  68                 # logging.debug("CM Cas%s déja calculé, portant le numéro %i"%(info,i))
  69                 break
  70         return __alc, __HxV
  71
  72     def storeValueInX(self, xValue, HxValue ):
  73         if self.__lenghtOR < 0:
  74             self.__lenghtOR = 2 * xValue.size + 2
  75             self.__initlnOR = self.__lenghtOR
  76         while len(self.__listOPCV) > self.__lenghtOR:
  77             # logging.debug("CM Réduction de la liste des cas à %i éléments par suppression du premier"%self.__lenghtOR)
  78             self.__listOPCV.pop(0)
  79         self.__listOPCV.append( (
  80             copy.copy(numpy.ravel(xValue)),
  81             copy.copy(HxValue),
  82             numpy.linalg.norm(xValue),
  83             ) )
  84
  85     def disable(self):
  86         self.__initlnOR = self.__lenghtOR
  87         self.__lenghtOR = 0
  88
  89     def enable(self):
  90         self.__lenghtOR = self.__initlnOR
  91
  92 # ==============================================================================
  93 class Operator:
  94     """
  95     Classe générale d'interface de type opérateur
  96     """
  97     NbCallsAsMatrix = 0
  98     NbCallsAsMethod = 0
  99     NbCallsOfCached = 0
 100     CM = CacheManager()
 101     #
 102     def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
 103         """
 104         On construit un objet de ce type en fournissant à l'aide de l'un des
 105         deux mots-clé, soit une fonction python, soit une matrice.
 106         Arguments :
 107         - fromMethod : argument de type fonction Python
 108         - fromMatrix : argument adapté au constructeur numpy.matrix
 109         - avoidingRedundancy : évite ou pas les calculs redondants
 110         """
 111         self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
 112         self.__AvoidRC = bool( avoidingRedundancy )
 113         if   fromMethod is not None:
 114             self.__Method = fromMethod
 115             self.__Matrix = None
 116             self.__Type   = "Method"
 117         elif fromMatrix is not None:
 118             self.__Method = None
 119             self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
 120             self.__Type   = "Matrix"
 121         else:
 122             self.__Method = None
 123             self.__Matrix = None
 124             self.__Type   = None
 125
 126     def disableAvoidingRedundancy(self):
 127         Operator.CM.disable()
 128
 129     def enableAvoidingRedundancy(self):
 130         if self.__AvoidRC:
 131             Operator.CM.enable()
 132         else:
 133             Operator.CM.disable()
 134
 135     def isType(self):
 136         return self.__Type
 137
 138     def appliedTo(self, xValue):
 139         """
 140         Permet de restituer le résultat de l'application de l'opérateur à un
 141         argument xValue. Cette méthode se contente d'appliquer, son argument
 142         devant a priori être du bon type.
 143         Arguments :
 144         - xValue : argument adapté pour appliquer l'opérateur
 145         """
 146         if self.__AvoidRC:
 147             __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
 148         else:
 149             __alreadyCalculated = False
 150         #
 151         if __alreadyCalculated:
 152             self.__addOneCacheCall()
 153             HxValue = __HxV
 154         else:
 155             if self.__Matrix is not None:
 156                 self.__addOneMatrixCall()
 157                 HxValue = self.__Matrix * xValue
 158             else:
 159                 self.__addOneMethodCall()
 160                 HxValue = self.__Method( xValue )
 161             if self.__AvoidRC:
 162                 Operator.CM.storeValueInX(xValue,HxValue)
 163         #
 164         return HxValue
 165
 166     def appliedControledFormTo(self, (xValue, uValue) ):
 167         """
 168         Permet de restituer le résultat de l'application de l'opérateur à une
 169         paire (xValue, uValue). Cette méthode se contente d'appliquer, son
 170         argument devant a priori être du bon type. Si la uValue est None,
 171         on suppose que l'opérateur ne s'applique qu'à xValue.
 172         Arguments :
 173         - xValue : argument X adapté pour appliquer l'opérateur
 174         - uValue : argument U adapté pour appliquer l'opérateur
 175         """
 176         if self.__Matrix is not None:
 177             self.__addOneMatrixCall()
 178             return self.__Matrix * xValue
 179         elif uValue is not None:
 180             self.__addOneMethodCall()
 181             return self.__Method( (xValue, uValue) )
 182         else:
 183             self.__addOneMethodCall()
 184             return self.__Method( xValue )
 185
 186     def appliedInXTo(self, (xNominal, xValue) ):
 187         """
 188         Permet de restituer le résultat de l'application de l'opérateur à un
 189         argument xValue, sachant que l'opérateur est valable en xNominal.
 190         Cette méthode se contente d'appliquer, son argument devant a priori
 191         être du bon type. Si l'opérateur est linéaire car c'est une matrice,
 192         alors il est valable en tout point nominal et il n'est pas nécessaire
 193         d'utiliser xNominal.
 194         Arguments : une liste contenant
 195         - xNominal : argument permettant de donner le point où l'opérateur
 196           est construit pour etre ensuite appliqué
 197         - xValue : argument adapté pour appliquer l'opérateur
 198         """
 199         if self.__Matrix is not None:
 200             self.__addOneMatrixCall()
 201             return self.__Matrix * xValue
 202         else:
 203             self.__addOneMethodCall()
 204             return self.__Method( (xNominal, xValue) )
 205
 206     def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
 207         """
 208         Permet de renvoyer l'opérateur sous la forme d'une matrice
 209         """
 210         if self.__Matrix is not None:
 211             self.__addOneMatrixCall()
 212             return self.__Matrix
 213         elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
 214             self.__addOneMethodCall()
 215             return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
 216         else:
 217             raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
 218
 219     def shape(self):
 220         """
 221         Renvoie la taille sous forme numpy si l'opérateur est disponible sous
 222         la forme d'une matrice
 223         """
 224         if self.__Matrix is not None:
 225             return self.__Matrix.shape
 226         else:
 227             raise ValueError("Matrix form of the operator is not available, nor the shape")
 228
 229     def nbcalls(self, which=None):
 230         """
 231         Renvoie les nombres d'évaluations de l'opérateur
 232         """
 233         __nbcalls = (
 234             self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
 235             self.__NbCallsAsMatrix,
 236             self.__NbCallsAsMethod,
 237             self.__NbCallsOfCached,
 238             Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
 239             Operator.NbCallsAsMatrix,
 240             Operator.NbCallsAsMethod,
 241             Operator.NbCallsOfCached,
 242             )
 243         if which is None: return __nbcalls
 244         else:             return __nbcalls[which]
 245
 246     def __addOneMatrixCall(self):
 247         self.__NbCallsAsMatrix   += 1 # Decompte local
 248         Operator.NbCallsAsMatrix += 1 # Decompte global
 249
 250     def __addOneMethodCall(self):
 251         self.__NbCallsAsMethod   += 1 # Decompte local
 252         Operator.NbCallsAsMethod += 1 # Decompte global
 253
 254     def __addOneCacheCall(self):
 255         self.__NbCallsOfCached   += 1 # Decompte local
 256         Operator.NbCallsOfCached += 1 # Decompte global
 257
 258 # ==============================================================================
 259 class Algorithm:
 260     """
 261     Classe générale d'interface de type algorithme
 262
 263     Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
 264     d'assimilation, en fournissant un container (dictionnaire) de variables
 265     persistantes initialisées, et des méthodes d'accès à ces variables stockées.
 266
 267     Une classe élémentaire d'algorithme doit implémenter la méthode "run".
 268     """
 269     def __init__(self, name):
 270         """
 271         L'initialisation présente permet de fabriquer des variables de stockage
 272         disponibles de manière générique dans les algorithmes élémentaires. Ces
 273         variables de stockage sont ensuite conservées dans un dictionnaire
 274         interne à l'objet, mais auquel on accède par la méthode "get".
 275
 276         Les variables prévues sont :
 277             - CostFunctionJ  : fonction-cout globale, somme des deux parties suivantes
 278             - CostFunctionJb : partie ébauche ou background de la fonction-cout
 279             - CostFunctionJo : partie observations de la fonction-cout
 280             - GradientOfCostFunctionJ  : gradient de la fonction-cout globale
 281             - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
 282             - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
 283             - CurrentState : état courant lors d'itérations
 284             - Analysis : l'analyse Xa
 285             - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
 286             - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
 287             - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
 288             - Innovation : l'innovation : d = Y - H(X)
 289             - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
 290             - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
 291             - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
 292             - MahalanobisConsistency : indicateur de consistance des covariances
 293             - OMA : Observation moins Analysis : Y - Xa
 294             - OMB : Observation moins Background : Y - Xb
 295             - AMB : Analysis moins Background : Xa - Xb
 296             - APosterioriCovariance : matrice A
 297             - APosterioriVariances : variances de la matrice A
 298             - APosterioriStandardDeviations : écart-types de la matrice A
 299             - APosterioriCorrelations : correlations de la matrice A
 300         On peut rajouter des variables à stocker dans l'initialisation de
 301         l'algorithme élémentaire qui va hériter de cette classe
 302         """
 303         logging.debug("%s Initialisation"%str(name))
 304         self._m = PlatformInfo.SystemUsage()
 305         #
 306         self._name = str( name )
 307         self._parameters = {"StoreSupplementaryCalculations":[]}
 308         self.__required_parameters = {}
 309         self.StoredVariables = {}
 310         #
 311         self.StoredVariables["CostFunctionJ"]                        = Persistence.OneScalar(name = "CostFunctionJ")
 312         self.StoredVariables["CostFunctionJb"]                       = Persistence.OneScalar(name = "CostFunctionJb")
 313         self.StoredVariables["CostFunctionJo"]                       = Persistence.OneScalar(name = "CostFunctionJo")
 314         self.StoredVariables["GradientOfCostFunctionJ"]              = Persistence.OneVector(name = "GradientOfCostFunctionJ")
 315         self.StoredVariables["GradientOfCostFunctionJb"]             = Persistence.OneVector(name = "GradientOfCostFunctionJb")
 316         self.StoredVariables["GradientOfCostFunctionJo"]             = Persistence.OneVector(name = "GradientOfCostFunctionJo")
 317         self.StoredVariables["CurrentState"]                         = Persistence.OneVector(name = "CurrentState")
 318         self.StoredVariables["Analysis"]                             = Persistence.OneVector(name = "Analysis")
 319         self.StoredVariables["IndexOfOptimum"]                       = Persistence.OneIndex(name = "IndexOfOptimum")
 320         self.StoredVariables["CurrentOptimum"]                       = Persistence.OneVector(name = "CurrentOptimum")
 321         self.StoredVariables["SimulatedObservationAtBackground"]     = Persistence.OneVector(name = "SimulatedObservationAtBackground")
 322         self.StoredVariables["SimulatedObservationAtCurrentState"]   = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
 323         self.StoredVariables["SimulatedObservationAtOptimum"]        = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
 324         self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
 325         self.StoredVariables["Innovation"]                           = Persistence.OneVector(name = "Innovation")
 326         self.StoredVariables["InnovationAtCurrentState"]             = Persistence.OneVector(name = "InnovationAtCurrentState")
 327         self.StoredVariables["SigmaObs2"]                            = Persistence.OneScalar(name = "SigmaObs2")
 328         self.StoredVariables["SigmaBck2"]                            = Persistence.OneScalar(name = "SigmaBck2")
 329         self.StoredVariables["MahalanobisConsistency"]               = Persistence.OneScalar(name = "MahalanobisConsistency")
 330         self.StoredVariables["OMA"]                                  = Persistence.OneVector(name = "OMA")
 331         self.StoredVariables["OMB"]                                  = Persistence.OneVector(name = "OMB")
 332         self.StoredVariables["BMA"]                                  = Persistence.OneVector(name = "BMA")
 333         self.StoredVariables["APosterioriCovariance"]                = Persistence.OneMatrix(name = "APosterioriCovariance")
 334         self.StoredVariables["APosterioriVariances"]                 = Persistence.OneVector(name = "APosterioriVariances")
 335         self.StoredVariables["APosterioriStandardDeviations"]        = Persistence.OneVector(name = "APosterioriStandardDeviations")
 336         self.StoredVariables["APosterioriCorrelations"]              = Persistence.OneMatrix(name = "APosterioriCorrelations")
 337         self.StoredVariables["SimulationQuantiles"]                  = Persistence.OneMatrix(name = "SimulationQuantiles")
 338
 339     def _pre_run(self):
 340         logging.debug("%s Lancement"%self._name)
 341         logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
 342         return 0
 343
 344     def _post_run(self,_oH=None):
 345         if ("StoreSupplementaryCalculations" in self._parameters) and \
 346             "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
 347             for _A in self.StoredVariables["APosterioriCovariance"]:
 348                 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
 349                     self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
 350                 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
 351                     self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
 352                 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
 353                     _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
 354                     _C = numpy.dot(_EI, numpy.dot(_A, _EI))
 355                     self.StoredVariables["APosterioriCorrelations"].store( _C )
 356         if _oH is not None:
 357             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)))
 358             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)))
 359         logging.debug("%s Taille mémoire utilisée de %.1f Mio"%(self._name, self._m.getUsedMemory("Mio")))
 360         logging.debug("%s Terminé"%self._name)
 361         return 0
 362
 363     def get(self, key=None):
 364         """
 365         Renvoie l'une des variables stockées identifiée par la clé, ou le
 366         dictionnaire de l'ensemble des variables disponibles en l'absence de
 367         clé. Ce sont directement les variables sous forme objet qui sont
 368         renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
 369         des classes de persistance.
 370         """
 371         if key is not None:
 372             return self.StoredVariables[key]
 373         else:
 374             return self.StoredVariables
 375
 376     def __contains__(self, key=None):
 377         """
 378         Vérifie si l'une des variables stockées est identifiée par la clé.
 379         """
 380         return (key in self.StoredVariables)
 381
 382     def keys(self):
 383         """
 384         Renvoie la liste des clés de variables stockées.
 385         """
 386         return self.StoredVariables.keys()
 387
 388     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
 389         """
 390         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
 391         sa forme mathématique la plus naturelle possible.
 392         """
 393         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
 394
 395     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
 396         """
 397         Permet de définir dans l'algorithme des paramètres requis et leurs
 398         caractéristiques par défaut.
 399         """
 400         if name is None:
 401             raise ValueError("A name is mandatory to define a required parameter.")
 402         #
 403         self.__required_parameters[name] = {
 404             "default"  : default,
 405             "typecast" : typecast,
 406             "minval"   : minval,
 407             "maxval"   : maxval,
 408             "listval"  : listval,
 409             "message"  : message,
 410             }
 411         logging.debug("%s %s (valeur par défaut = %s)"%(self._name, message, self.setParameterValue(name)))
 412
 413     def getRequiredParameters(self, noDetails=True):
 414         """
 415         Renvoie la liste des noms de paramètres requis ou directement le
 416         dictionnaire des paramètres requis.
 417         """
 418         if noDetails:
 419             ks = self.__required_parameters.keys()
 420             ks.sort()
 421             return ks
 422         else:
 423             return self.__required_parameters
 424
 425     def setParameterValue(self, name=None, value=None):
 426         """
 427         Renvoie la valeur d'un paramètre requis de manière contrôlée
 428         """
 429         default  = self.__required_parameters[name]["default"]
 430         typecast = self.__required_parameters[name]["typecast"]
 431         minval   = self.__required_parameters[name]["minval"]
 432         maxval   = self.__required_parameters[name]["maxval"]
 433         listval  = self.__required_parameters[name]["listval"]
 434         #
 435         if value is None and default is None:
 436             __val = None
 437         elif value is None and default is not None:
 438             if typecast is None: __val = default
 439             else:                __val = typecast( default )
 440         else:
 441             if typecast is None: __val = value
 442             else:                __val = typecast( value )
 443         #
 444         if minval is not None and (numpy.array(__val, float) < minval).any():
 445             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
 446         if maxval is not None and (numpy.array(__val, float) > maxval).any():
 447             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
 448         if listval is not None:
 449             if typecast is list or typecast is tuple or type(__val) is list or type(__val) is tuple:
 450                 for v in __val:
 451                     if v not in listval:
 452                         raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
 453             elif __val not in listval:
 454                 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
 455         return __val
 456
 457     def setParameters(self, fromDico={}):
 458         """
 459         Permet de stocker les paramètres reçus dans le dictionnaire interne.
 460         """
 461         self._parameters.update( fromDico )
 462         for k in self.__required_parameters.keys():
 463             if k in fromDico.keys():
 464                 self._parameters[k] = self.setParameterValue(k,fromDico[k])
 465             else:
 466                 self._parameters[k] = self.setParameterValue(k)
 467             logging.debug("%s %s : %s"%(self._name, self.__required_parameters[k]["message"], self._parameters[k]))
 468
 469 # ==============================================================================
 470 class Diagnostic:
 471     """
 472     Classe générale d'interface de type diagnostic
 473
 474     Ce template s'utilise de la manière suivante : il sert de classe "patron" en
 475     même temps que l'une des classes de persistance, comme "OneScalar" par
 476     exemple.
 477
 478     Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
 479     méthode "_formula" pour écrire explicitement et proprement la formule pour
 480     l'écriture mathématique du calcul du diagnostic (méthode interne non
 481     publique), et "calculate" pour activer la précédente tout en ayant vérifié
 482     et préparé les données, et pour stocker les résultats à chaque pas (méthode
 483     externe d'activation).
 484     """
 485     def __init__(self, name = "", parameters = {}):
 486         self.name       = str(name)
 487         self.parameters = dict( parameters )
 488
 489     def _formula(self, *args):
 490         """
 491         Doit implémenter l'opération élémentaire de diagnostic sous sa forme
 492         mathématique la plus naturelle possible.
 493         """
 494         raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
 495
 496     def calculate(self, *args):
 497         """
 498         Active la formule de calcul avec les arguments correctement rangés
 499         """
 500         raise NotImplementedError("Diagnostic activation method has not been implemented!")
 501
 502 # ==============================================================================
 503 class Covariance:
 504     """
 505     Classe générale d'interface de type covariance
 506     """
 507     def __init__(self,
 508             name          = "GenericCovariance",
 509             asCovariance  = None,
 510             asEyeByScalar = None,
 511             asEyeByVector = None,
 512             asCovObject   = None,
 513             toBeChecked   = False,
 514             ):
 515         """
 516         Permet de définir une covariance :
 517         - asCovariance : entrée des données, comme une matrice compatible avec
 518           le constructeur de numpy.matrix
 519         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 520           multiplicatif d'une matrice de corrélation identité, aucune matrice
 521           n'étant donc explicitement à donner
 522         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 523           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 524           n'étant donc explicitement à donner
 525         - asCovObject : entrée des données comme un objet python, qui a les
 526           methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
 527           "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
 528           "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
 529         - toBeChecked : booléen indiquant si le caractère SDP de la matrice
 530           pleine doit être vérifié
 531         """
 532         self.__name       = str(name)
 533         self.__check      = bool(toBeChecked)
 534         #
 535         self.__C          = None
 536         self.__is_scalar  = False
 537         self.__is_vector  = False
 538         self.__is_matrix  = False
 539         self.__is_object  = False
 540         if asEyeByScalar is not None:
 541             self.__is_scalar = True
 542             self.__C         = numpy.abs( float(asEyeByScalar) )
 543             self.shape       = (0,0)
 544             self.size        = 0
 545         elif asEyeByVector is not None:
 546             self.__is_vector = True
 547             self.__C         = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
 548             self.shape       = (self.__C.size,self.__C.size)
 549             self.size        = self.__C.size**2
 550         elif asCovariance is not None:
 551             self.__is_matrix = True
 552             self.__C         = numpy.matrix( asCovariance, float )
 553             self.shape       = self.__C.shape
 554             self.size        = self.__C.size
 555         elif asCovObject is not None:
 556             self.__is_object = True
 557             self.__C         = asCovObject
 558             for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
 559                 if not hasattr(self.__C,at):
 560                     raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
 561             if hasattr(self.__C,"shape"):
 562                 self.shape       = self.__C.shape
 563             else:
 564                 self.shape       = (0,0)
 565             if hasattr(self.__C,"size"):
 566                 self.size        = self.__C.size
 567             else:
 568                 self.size        = 0
 569         else:
 570             pass
 571             # 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)
 572         #
 573         self.__validate()
 574
 575     def __validate(self):
 576         if self.ismatrix() and min(self.shape) != max(self.shape):
 577             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))
 578         if self.isobject() and min(self.shape) != max(self.shape):
 579             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))
 580         if self.isscalar() and self.__C <= 0:
 581             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
 582         if self.isvector() and (self.__C <= 0).any():
 583             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
 584         if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
 585             try:
 586                 L = numpy.linalg.cholesky( self.__C )
 587             except:
 588                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 589         if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
 590             try:
 591                 L = self.__C.cholesky()
 592             except:
 593                 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 594
 595     def isscalar(self):
 596         return self.__is_scalar
 597
 598     def isvector(self):
 599         return self.__is_vector
 600
 601     def ismatrix(self):
 602         return self.__is_matrix
 603
 604     def isobject(self):
 605         return self.__is_object
 606
 607     def getI(self):
 608         if   self.ismatrix():
 609             return Covariance(self.__name+"I", asCovariance  = self.__C.I )
 610         elif self.isvector():
 611             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
 612         elif self.isscalar():
 613             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
 614         elif self.isobject():
 615             return Covariance(self.__name+"I", asCovObject   = self.__C.getI() )
 616         else:
 617             return None # Indispensable
 618
 619     def getT(self):
 620         if   self.ismatrix():
 621             return Covariance(self.__name+"T", asCovariance  = self.__C.T )
 622         elif self.isvector():
 623             return Covariance(self.__name+"T", asEyeByVector = self.__C )
 624         elif self.isscalar():
 625             return Covariance(self.__name+"T", asEyeByScalar = self.__C )
 626         elif self.isobject():
 627             return Covariance(self.__name+"T", asCovObject   = self.__C.getT() )
 628
 629     def cholesky(self):
 630         if   self.ismatrix():
 631             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__C) )
 632         elif self.isvector():
 633             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
 634         elif self.isscalar():
 635             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
 636         elif self.isobject() and hasattr(self.__C,"cholesky"):
 637             return Covariance(self.__name+"C", asCovObject   = self.__C.cholesky() )
 638
 639     def choleskyI(self):
 640         if   self.ismatrix():
 641             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__C).I )
 642         elif self.isvector():
 643             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
 644         elif self.isscalar():
 645             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
 646         elif self.isobject() and hasattr(self.__C,"choleskyI"):
 647             return Covariance(self.__name+"H", asCovObject   = self.__C.choleskyI() )
 648
 649     def diag(self, msize=None):
 650         if   self.ismatrix():
 651             return numpy.diag(self.__C)
 652         elif self.isvector():
 653             return self.__C
 654         elif self.isscalar():
 655             if msize is None:
 656                 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,))
 657             else:
 658                 return self.__C * numpy.ones(int(msize))
 659         elif self.isobject():
 660             return self.__C.diag()
 661
 662     def asfullmatrix(self, msize=None):
 663         if   self.ismatrix():
 664             return self.__C
 665         elif self.isvector():
 666             return numpy.matrix( numpy.diag(self.__C), float )
 667         elif self.isscalar():
 668             if msize is None:
 669                 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,))
 670             else:
 671                 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
 672         elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
 673             return self.__C.asfullmatrix()
 674
 675     def trace(self, msize=None):
 676         if   self.ismatrix():
 677             return numpy.trace(self.__C)
 678         elif self.isvector():
 679             return float(numpy.sum(self.__C))
 680         elif self.isscalar():
 681             if msize is None:
 682                 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,))
 683             else:
 684                 return self.__C * int(msize)
 685         elif self.isobject():
 686             return self.__C.trace()
 687
 688     def __repr__(self):
 689         return repr(self.__C)
 690
 691     def __str__(self):
 692         return str(self.__C)
 693
 694     def __add__(self, other):
 695         if   self.ismatrix() or self.isobject():
 696             return self.__C + numpy.asmatrix(other)
 697         elif self.isvector() or self.isscalar():
 698             _A = numpy.asarray(other)
 699             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
 700             return numpy.asmatrix(_A)
 701
 702     def __radd__(self, other):
 703         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
 704
 705     def __sub__(self, other):
 706         if   self.ismatrix() or self.isobject():
 707             return self.__C - numpy.asmatrix(other)
 708         elif self.isvector() or self.isscalar():
 709             _A = numpy.asarray(other)
 710             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
 711             return numpy.asmatrix(_A)
 712
 713     def __rsub__(self, other):
 714         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
 715
 716     def __neg__(self):
 717         return - self.__C
 718
 719     def __mul__(self, other):
 720         if   self.ismatrix() and isinstance(other,numpy.matrix):
 721             return self.__C * other
 722         elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
 723                                or isinstance(other,list) \
 724                                or isinstance(other,tuple)):
 725             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 726                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
 727             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 728                 return self.__C * numpy.asmatrix(other)
 729             else:
 730                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
 731         elif self.isvector() and (isinstance(other,numpy.matrix) \
 732                                or isinstance(other,numpy.ndarray) \
 733                                or isinstance(other,list) \
 734                                or isinstance(other,tuple)):
 735             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 736                 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
 737             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 738                 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
 739             else:
 740                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 741         elif self.isscalar() and isinstance(other,numpy.matrix):
 742             return self.__C * other
 743         elif self.isscalar() and (isinstance(other,numpy.ndarray) \
 744                                or isinstance(other,list) \
 745                                or isinstance(other,tuple)):
 746             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
 747                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
 748             else:
 749                 return self.__C * numpy.asmatrix(other)
 750         elif self.isobject():
 751             return self.__C.__mul__(other)
 752         else:
 753             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
 754
 755     def __rmul__(self, other):
 756         if self.ismatrix() and isinstance(other,numpy.matrix):
 757             return other * self.__C
 758         elif self.isvector() and isinstance(other,numpy.matrix):
 759             if numpy.ravel(other).size == self.shape[0]: # Vecteur
 760                 return numpy.asmatrix(numpy.ravel(other) * self.__C)
 761             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
 762                 return numpy.asmatrix(numpy.array(other) * self.__C)
 763             else:
 764                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 765         elif self.isscalar() and isinstance(other,numpy.matrix):
 766             return other * self.__C
 767         elif self.isobject():
 768             return self.__C.__rmul__(other)
 769         else:
 770             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
 771
 772     def __len__(self):
 773         return self.shape[0]
 774
 775 # ==============================================================================
 776 def CostFunction3D(
 777     _x,
 778     _Hm  = None,  # Pour simuler Hm(x) : HO["Direct"].appliedTo
 779     _HmX = None,  # Simulation déjà faite de Hm(x)
 780     _arg = None,  # Arguments supplementaires pour Hm, sous la forme d'un tuple
 781     _BI  = None,
 782     _RI  = None,
 783     _Xb  = None,
 784     _Y   = None,
 785     _SIV = False, # A résorber pour la 8.0
 786     _SSC = [],    # self._parameters["StoreSupplementaryCalculations"]
 787     _nPS = 0,     # nbPreviousSteps
 788     _QM  = "DA",  # QualityMeasure
 789     _SSV = {},    # Entrée et/ou sortie : self.StoredVariables
 790     _fRt = False, # Restitue ou pas la sortie étendue
 791     _sSc = True,  # Stocke ou pas les SSC
 792     ):
 793     """
 794     Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
 795     et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
 796     DFO, QuantileRegression
 797     """
 798     if not _sSc:
 799         _SIV = False
 800         _SSC = {}
 801     else:
 802         for k in ["CostFunctionJ",
 803                   "CostFunctionJb",
 804                   "CostFunctionJo",
 805                   "CurrentOptimum",
 806                   "CurrentState",
 807                   "IndexOfOptimum",
 808                   "SimulatedObservationAtCurrentOptimum",
 809                   "SimulatedObservationAtCurrentState",
 810                   ]:
 811             if k not in _SSV:
 812                 _SSV[k] = []
 813             if hasattr(_SSV[k],"store"):
 814                 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
 815     #
 816     _X  = numpy.asmatrix(numpy.ravel( _x )).T
 817     if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
 818         _SSV["CurrentState"].append( _X )
 819     #
 820     if _HmX is not None:
 821         _HX = _HmX
 822     else:
 823         if _Hm is None:
 824             raise ValueError("%s Operator has to be defined."%(self.__name,))
 825         if _arg is None:
 826             _HX = _Hm( _X )
 827         else:
 828             _HX = _Hm( _X, *_arg )
 829     _HX = numpy.asmatrix(numpy.ravel( _HX )).T
 830     #
 831     if "SimulatedObservationAtCurrentState" in _SSC or \
 832        "SimulatedObservationAtCurrentOptimum" in _SSC:
 833         _SSV["SimulatedObservationAtCurrentState"].append( _HX )
 834     #
 835     if numpy.any(numpy.isnan(_HX)):
 836         Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
 837     else:
 838         _Y   = numpy.asmatrix(numpy.ravel( _Y )).T
 839         if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
 840             if _BI is None or _RI is None:
 841                 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
 842             _Xb  = numpy.asmatrix(numpy.ravel( _Xb )).T
 843             Jb  = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
 844             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
 845         elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
 846             if _RI is None:
 847                 raise ValueError("Observation error covariance matrix has to be properly defined!")
 848             Jb  = 0.
 849             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
 850         elif _QM in ["LeastSquares", "LS", "L2"]:
 851             Jb  = 0.
 852             Jo  = 0.5 * (_Y - _HX).T * (_Y - _HX)
 853         elif _QM in ["AbsoluteValue", "L1"]:
 854             Jb  = 0.
 855             Jo  = numpy.sum( numpy.abs(_Y - _HX) )
 856         elif _QM in ["MaximumError", "ME"]:
 857             Jb  = 0.
 858             Jo  = numpy.max( numpy.abs(_Y - _HX) )
 859         elif _QM in ["QR", "Null"]:
 860             Jb  = 0.
 861             Jo  = 0.
 862         else:
 863             raise ValueError("Unknown asked quality measure!")
 864         #
 865         J   = float( Jb ) + float( Jo )
 866     #
 867     if _sSc:
 868         _SSV["CostFunctionJb"].append( Jb )
 869         _SSV["CostFunctionJo"].append( Jo )
 870         _SSV["CostFunctionJ" ].append( J )
 871     #
 872     if "IndexOfOptimum" in _SSC or \
 873        "CurrentOptimum" in _SSC or \
 874        "SimulatedObservationAtCurrentOptimum" in _SSC:
 875         IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
 876     if "IndexOfOptimum" in _SSC:
 877         _SSV["IndexOfOptimum"].append( IndexMin )
 878     if "CurrentOptimum" in _SSC:
 879         _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
 880     if "SimulatedObservationAtCurrentOptimum" in _SSC:
 881         _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
 882     #
 883     if _fRt:
 884         return _SSV
 885     else:
 886         if _QM in ["QR"]: # Pour le QuantileRegression
 887             return _HX
 888         else:
 889             return J
 890
 891 # ==============================================================================
 892 if __name__ == "__main__":
 893     print '\n AUTODIAGNOSTIC \n'