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