src/daComposant/daCore/BasicObjects.py

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