src/daComposant/daCore/BasicObjects.py

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