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             ):
 514         """
 515         Permet de définir une covariance :
 516         - asCovariance : entrée des données, comme une matrice compatible avec
 517           le constructeur de numpy.matrix
 518         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 519           multiplicatif d'une matrice de corrélation identité, aucune matrice
 520           n'étant donc explicitement à donner
 521         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 522           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 523           n'étant donc explicitement à donner
 524         - asCovObject : entrée des données comme un objet python, qui a les
 525           methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
 526           "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
 527           "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
 528         """
 529         self.__name       = str(name)
 530         #
 531         self.__C          = None
 532         self.__is_scalar  = False
 533         self.__is_vector  = False
 534         self.__is_matrix  = False
 535         self.__is_object  = False
 536         if asEyeByScalar is not None:
 537             self.__is_scalar = True
 538             self.__C         = numpy.abs( float(asEyeByScalar) )
 539             self.shape       = (0,0)
 540             self.size        = 0
 541         elif asEyeByVector is not None:
 542             self.__is_vector = True
 543             self.__C         = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
 544             self.shape       = (self.__C.size,self.__C.size)
 545             self.size        = self.__C.size**2
 546         elif asCovariance is not None:
 547             self.__is_matrix = True
 548             self.__C         = numpy.matrix( asCovariance, float )
 549             self.shape       = self.__C.shape
 550             self.size        = self.__C.size
 551         elif asCovObject is not None:
 552             self.__is_object = True
 553             self.__C         = asCovObject
 554             for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
 555                 if not hasattr(self.__C,at):
 556                     raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
 557             if hasattr(self.__C,"shape"):
 558                 self.shape       = self.__C.shape
 559             else:
 560                 self.shape       = (0,0)
 561             if hasattr(self.__C,"size"):
 562                 self.size        = self.__C.size
 563             else:
 564                 self.size        = 0
 565         else:
 566             pass
 567             # 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)
 568         #
 569         self.__validate()
 570
 571     def __validate(self):
 572         if self.ismatrix() and min(self.shape) != max(self.shape):
 573             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))
 574         if self.isobject() and min(self.shape) != max(self.shape):
 575             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))
 576         if self.isscalar() and self.__C <= 0:
 577             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
 578         if self.isvector() and (self.__C <= 0).any():
 579             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
 580         if self.ismatrix() and logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
 581             try:
 582                 L = numpy.linalg.cholesky( self.__C )
 583             except:
 584                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 585
 586     def isscalar(self):
 587         return self.__is_scalar
 588
 589     def isvector(self):
 590         return self.__is_vector
 591
 592     def ismatrix(self):
 593         return self.__is_matrix
 594
 595     def isobject(self):
 596         return self.__is_object
 597
 598     def getI(self):
 599         if   self.ismatrix():
 600             return Covariance(self.__name+"I", asCovariance  = self.__C.I )
 601         elif self.isvector():
 602             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
 603         elif self.isscalar():
 604             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
 605         elif self.isobject():
 606             return Covariance(self.__name+"I", asCovObject   = self.__C.getI() )
 607         else:
 608             return None # Indispensable
 609
 610     def getT(self):
 611         if   self.ismatrix():
 612             return Covariance(self.__name+"T", asCovariance  = self.__C.T )
 613         elif self.isvector():
 614             return Covariance(self.__name+"T", asEyeByVector = self.__C )
 615         elif self.isscalar():
 616             return Covariance(self.__name+"T", asEyeByScalar = self.__C )
 617         elif self.isobject():
 618             return Covariance(self.__name+"T", asCovObject   = self.__C.getT() )
 619
 620     def cholesky(self):
 621         if   self.ismatrix():
 622             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__C) )
 623         elif self.isvector():
 624             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
 625         elif self.isscalar():
 626             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
 627         elif self.isobject() and hasattr(self.__C,"cholesky"):
 628             return Covariance(self.__name+"C", asCovObject   = self.__C.cholesky() )
 629
 630     def choleskyI(self):
 631         if   self.ismatrix():
 632             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__C).I )
 633         elif self.isvector():
 634             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
 635         elif self.isscalar():
 636             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
 637         elif self.isobject() and hasattr(self.__C,"choleskyI"):
 638             return Covariance(self.__name+"H", asCovObject   = self.__C.choleskyI() )
 639
 640     def diag(self, msize=None):
 641         if   self.ismatrix():
 642             return numpy.diag(self.__C)
 643         elif self.isvector():
 644             return self.__C
 645         elif self.isscalar():
 646             if msize is None:
 647                 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,))
 648             else:
 649                 return self.__C * numpy.ones(int(msize))
 650         elif self.isobject():
 651             return self.__C.diag()
 652
 653     def asfullmatrix(self, msize=None):
 654         if   self.ismatrix():
 655             return self.__C
 656         elif self.isvector():
 657             return numpy.matrix( numpy.diag(self.__C), float )
 658         elif self.isscalar():
 659             if msize is None:
 660                 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,))
 661             else:
 662                 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
 663         elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
 664             return self.__C.asfullmatrix()
 665
 666     def trace(self, msize=None):
 667         if   self.ismatrix():
 668             return numpy.trace(self.__C)
 669         elif self.isvector():
 670             return float(numpy.sum(self.__C))
 671         elif self.isscalar():
 672             if msize is None:
 673                 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,))
 674             else:
 675                 return self.__C * int(msize)
 676         elif self.isobject():
 677             return self.__C.trace()
 678
 679     def __repr__(self):
 680         return repr(self.__C)
 681
 682     def __str__(self):
 683         return str(self.__C)
 684
 685     def __add__(self, other):
 686         if   self.ismatrix() or self.isobject():
 687             return self.__C + numpy.asmatrix(other)
 688         elif self.isvector() or self.isscalar():
 689             _A = numpy.asarray(other)
 690             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
 691             return numpy.asmatrix(_A)
 692
 693     def __radd__(self, other):
 694         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
 695
 696     def __sub__(self, other):
 697         if   self.ismatrix() or self.isobject():
 698             return self.__C - numpy.asmatrix(other)
 699         elif self.isvector() or self.isscalar():
 700             _A = numpy.asarray(other)
 701             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
 702             return numpy.asmatrix(_A)
 703
 704     def __rsub__(self, other):
 705         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
 706
 707     def __neg__(self):
 708         return - self.__C
 709
 710     def __mul__(self, other):
 711         if   self.ismatrix() and isinstance(other,numpy.matrix):
 712             return self.__C * other
 713         elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
 714                                or isinstance(other,list) \
 715                                or isinstance(other,tuple)):
 716             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 717                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
 718             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 719                 return self.__C * numpy.asmatrix(other)
 720             else:
 721                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
 722         elif self.isvector() and (isinstance(other,numpy.matrix) \
 723                                or isinstance(other,numpy.ndarray) \
 724                                or isinstance(other,list) \
 725                                or isinstance(other,tuple)):
 726             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 727                 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
 728             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 729                 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
 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 self.__C * other
 734         elif self.isscalar() and (isinstance(other,numpy.ndarray) \
 735                                or isinstance(other,list) \
 736                                or isinstance(other,tuple)):
 737             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
 738                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
 739             else:
 740                 return self.__C * numpy.asmatrix(other)
 741         elif self.isobject():
 742             return self.__C.__mul__(other)
 743         else:
 744             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
 745
 746     def __rmul__(self, other):
 747         if self.ismatrix() and isinstance(other,numpy.matrix):
 748             return other * self.__C
 749         elif self.isvector() and isinstance(other,numpy.matrix):
 750             if numpy.ravel(other).size == self.shape[0]: # Vecteur
 751                 return numpy.asmatrix(numpy.ravel(other) * self.__C)
 752             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
 753                 return numpy.asmatrix(numpy.array(other) * self.__C)
 754             else:
 755                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 756         elif self.isscalar() and isinstance(other,numpy.matrix):
 757             return other * self.__C
 758         elif self.isobject():
 759             return self.__C.__rmul__(other)
 760         else:
 761             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
 762
 763     def __len__(self):
 764         return self.shape[0]
 765
 766 # ==============================================================================
 767 def CostFunction3D(
 768     _x,
 769     _Hm  = None,  # Pour simuler Hm(x) : HO["Direct"].appliedTo
 770     _HmX = None,  # Simulation déjà faite de Hm(x)
 771     _arg = None,  # Arguments supplementaires pour Hm, sous la forme d'un tuple
 772     _BI  = None,
 773     _RI  = None,
 774     _Xb  = None,
 775     _Y   = None,
 776     _SIV = False, # A résorber pour la 8.0
 777     _SSC = [],    # self._parameters["StoreSupplementaryCalculations"]
 778     _nPS = 0,     # nbPreviousSteps
 779     _QM  = "DA",  # QualityMeasure
 780     _SSV = {},    # Entrée et/ou sortie : self.StoredVariables
 781     _fRt = False, # Restitue ou pas la sortie étendue
 782     _sSc = True,  # Stocke ou pas les SSC
 783     ):
 784     """
 785     Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
 786     et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
 787     DFO, QuantileRegression
 788     """
 789     if not _sSc:
 790         _SIV = False
 791         _SSC = {}
 792     else:
 793         for k in ["CostFunctionJ",
 794                   "CostFunctionJb",
 795                   "CostFunctionJo",
 796                   "CurrentOptimum",
 797                   "CurrentState",
 798                   "IndexOfOptimum",
 799                   "SimulatedObservationAtCurrentOptimum",
 800                   "SimulatedObservationAtCurrentState",
 801                   ]:
 802             if k not in _SSV:
 803                 _SSV[k] = []
 804             if hasattr(_SSV[k],"store"):
 805                 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
 806     #
 807     _X  = numpy.asmatrix(numpy.ravel( _x )).T
 808     if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
 809         _SSV["CurrentState"].append( _X )
 810     #
 811     if _HmX is not None:
 812         _HX = _HmX
 813     else:
 814         if _Hm is None:
 815             raise ValueError("%s Operator has to be defined."%(self.__name,))
 816         if _arg is None:
 817             _HX = _Hm( _X )
 818         else:
 819             _HX = _Hm( _X, *_arg )
 820     _HX = numpy.asmatrix(numpy.ravel( _HX )).T
 821     #
 822     if "SimulatedObservationAtCurrentState" in _SSC or \
 823        "SimulatedObservationAtCurrentOptimum" in _SSC:
 824         _SSV["SimulatedObservationAtCurrentState"].append( _HX )
 825     #
 826     if numpy.any(numpy.isnan(_HX)):
 827         Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
 828     else:
 829         _Y   = numpy.asmatrix(numpy.ravel( _Y )).T
 830         if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
 831             if _BI is None or _RI is None:
 832                 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
 833             _Xb  = numpy.asmatrix(numpy.ravel( _Xb )).T
 834             Jb  = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
 835             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
 836         elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
 837             if _RI is None:
 838                 raise ValueError("Observation error covariance matrix has to be properly defined!")
 839             Jb  = 0.
 840             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
 841         elif _QM in ["LeastSquares", "LS", "L2"]:
 842             Jb  = 0.
 843             Jo  = 0.5 * (_Y - _HX).T * (_Y - _HX)
 844         elif _QM in ["AbsoluteValue", "L1"]:
 845             Jb  = 0.
 846             Jo  = numpy.sum( numpy.abs(_Y - _HX) )
 847         elif _QM in ["MaximumError", "ME"]:
 848             Jb  = 0.
 849             Jo  = numpy.max( numpy.abs(_Y - _HX) )
 850         elif _QM in ["QR", "Null"]:
 851             Jb  = 0.
 852             Jo  = 0.
 853         else:
 854             raise ValueError("Unknown asked quality measure!")
 855         #
 856         J   = float( Jb ) + float( Jo )
 857     #
 858     if _sSc:
 859         _SSV["CostFunctionJb"].append( Jb )
 860         _SSV["CostFunctionJo"].append( Jo )
 861         _SSV["CostFunctionJ" ].append( J )
 862     #
 863     if "IndexOfOptimum" in _SSC or \
 864        "CurrentOptimum" in _SSC or \
 865        "SimulatedObservationAtCurrentOptimum" in _SSC:
 866         IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
 867     if "IndexOfOptimum" in _SSC:
 868         _SSV["IndexOfOptimum"].append( IndexMin )
 869     if "CurrentOptimum" in _SSC:
 870         _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
 871     if "SimulatedObservationAtCurrentOptimum" in _SSC:
 872         _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
 873     #
 874     if _fRt:
 875         return _SSV
 876     else:
 877         if _QM in ["QR"]: # Pour le QuantileRegression
 878             return _HX
 879         else:
 880             return J
 881
 882 # ==============================================================================
 883 if __name__ == "__main__":
 884     print '\n AUTODIAGNOSTIC \n'