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é à ê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["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         self.StoredVariables["Residu"]                               = Persistence.OneScalar(name = "Residu")
 350
 351     def _pre_run(self):
 352         "Pré-calcul"
 353         logging.debug("%s Lancement", self._name)
 354         logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
 355         return 0
 356
 357     def _post_run(self,_oH=None):
 358         "Post-calcul"
 359         if ("StoreSupplementaryCalculations" in self._parameters) and \
 360             "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
 361             for _A in self.StoredVariables["APosterioriCovariance"]:
 362                 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
 363                     self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
 364                 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
 365                     self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
 366                 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
 367                     _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
 368                     _C = numpy.dot(_EI, numpy.dot(_A, _EI))
 369                     self.StoredVariables["APosterioriCorrelations"].store( _C )
 370         if _oH is not None:
 371             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))
 372             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))
 373         logging.debug("%s Taille mémoire utilisée de %.1f Mio", self._name, self._m.getUsedMemory("Mio"))
 374         logging.debug("%s Terminé", self._name)
 375         return 0
 376
 377     def get(self, key=None):
 378         """
 379         Renvoie l'une des variables stockées identifiée par la clé, ou le
 380         dictionnaire de l'ensemble des variables disponibles en l'absence de
 381         clé. Ce sont directement les variables sous forme objet qui sont
 382         renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
 383         des classes de persistance.
 384         """
 385         if key is not None:
 386             return self.StoredVariables[key]
 387         else:
 388             return self.StoredVariables
 389
 390     def __contains__(self, key=None):
 391         "D.__contains__(k) -> True if D has a key k, else False"
 392         return key in self.StoredVariables
 393
 394     def keys(self):
 395         "D.keys() -> list of D's keys"
 396         return self.StoredVariables.keys()
 397
 398     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
 399         """
 400         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
 401         sa forme mathématique la plus naturelle possible.
 402         """
 403         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
 404
 405     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
 406         """
 407         Permet de définir dans l'algorithme des paramètres requis et leurs
 408         caractéristiques par défaut.
 409         """
 410         if name is None:
 411             raise ValueError("A name is mandatory to define a required parameter.")
 412         #
 413         self.__required_parameters[name] = {
 414             "default"  : default,
 415             "typecast" : typecast,
 416             "minval"   : minval,
 417             "maxval"   : maxval,
 418             "listval"  : listval,
 419             "message"  : message,
 420             }
 421         logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
 422
 423     def getRequiredParameters(self, noDetails=True):
 424         """
 425         Renvoie la liste des noms de paramètres requis ou directement le
 426         dictionnaire des paramètres requis.
 427         """
 428         if noDetails:
 429             ks = self.__required_parameters.keys()
 430             ks.sort()
 431             return ks
 432         else:
 433             return self.__required_parameters
 434
 435     def setParameterValue(self, name=None, value=None):
 436         """
 437         Renvoie la valeur d'un paramètre requis de manière contrôlée
 438         """
 439         default  = self.__required_parameters[name]["default"]
 440         typecast = self.__required_parameters[name]["typecast"]
 441         minval   = self.__required_parameters[name]["minval"]
 442         maxval   = self.__required_parameters[name]["maxval"]
 443         listval  = self.__required_parameters[name]["listval"]
 444         #
 445         if value is None and default is None:
 446             __val = None
 447         elif value is None and default is not None:
 448             if typecast is None: __val = default
 449             else:                __val = typecast( default )
 450         else:
 451             if typecast is None: __val = value
 452             else:                __val = typecast( value )
 453         #
 454         if minval is not None and (numpy.array(__val, float) < minval).any():
 455             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
 456         if maxval is not None and (numpy.array(__val, float) > maxval).any():
 457             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
 458         if listval is not None:
 459             if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
 460                 for v in __val:
 461                     if v not in listval:
 462                         raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
 463             elif __val not in listval:
 464                 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
 465         return __val
 466
 467     def setParameters(self, fromDico={}):
 468         """
 469         Permet de stocker les paramètres reçus dans le dictionnaire interne.
 470         """
 471         self._parameters.update( fromDico )
 472         for k in self.__required_parameters.keys():
 473             if k in fromDico.keys():
 474                 self._parameters[k] = self.setParameterValue(k,fromDico[k])
 475             else:
 476                 self._parameters[k] = self.setParameterValue(k)
 477             logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
 478
 479 # ==============================================================================
 480 class Diagnostic(object):
 481     """
 482     Classe générale d'interface de type diagnostic
 483
 484     Ce template s'utilise de la manière suivante : il sert de classe "patron" en
 485     même temps que l'une des classes de persistance, comme "OneScalar" par
 486     exemple.
 487
 488     Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
 489     méthode "_formula" pour écrire explicitement et proprement la formule pour
 490     l'écriture mathématique du calcul du diagnostic (méthode interne non
 491     publique), et "calculate" pour activer la précédente tout en ayant vérifié
 492     et préparé les données, et pour stocker les résultats à chaque pas (méthode
 493     externe d'activation).
 494     """
 495     def __init__(self, name = "", parameters = {}):
 496         "Initialisation"
 497         self.name       = str(name)
 498         self.parameters = dict( parameters )
 499
 500     def _formula(self, *args):
 501         """
 502         Doit implémenter l'opération élémentaire de diagnostic sous sa forme
 503         mathématique la plus naturelle possible.
 504         """
 505         raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
 506
 507     def calculate(self, *args):
 508         """
 509         Active la formule de calcul avec les arguments correctement rangés
 510         """
 511         raise NotImplementedError("Diagnostic activation method has not been implemented!")
 512
 513 # ==============================================================================
 514 class Covariance(object):
 515     """
 516     Classe générale d'interface de type covariance
 517     """
 518     def __init__(self,
 519                  name          = "GenericCovariance",
 520                  asCovariance  = None,
 521                  asEyeByScalar = None,
 522                  asEyeByVector = None,
 523                  asCovObject   = None,
 524                  toBeChecked   = False,
 525                 ):
 526         """
 527         Permet de définir une covariance :
 528         - asCovariance : entrée des données, comme une matrice compatible avec
 529           le constructeur de numpy.matrix
 530         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 531           multiplicatif d'une matrice de corrélation identité, aucune matrice
 532           n'étant donc explicitement à donner
 533         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 534           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 535           n'étant donc explicitement à donner
 536         - asCovObject : entrée des données comme un objet python, qui a les
 537           methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
 538           "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
 539           "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
 540         - toBeChecked : booléen indiquant si le caractère SDP de la matrice
 541           pleine doit être vérifié
 542         """
 543         self.__name       = str(name)
 544         self.__check      = bool(toBeChecked)
 545         #
 546         self.__C          = None
 547         self.__is_scalar  = False
 548         self.__is_vector  = False
 549         self.__is_matrix  = False
 550         self.__is_object  = False
 551         if asEyeByScalar is not None:
 552             self.__is_scalar = True
 553             self.__C         = numpy.abs( float(asEyeByScalar) )
 554             self.shape       = (0,0)
 555             self.size        = 0
 556         elif asEyeByVector is not None:
 557             self.__is_vector = True
 558             self.__C         = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
 559             self.shape       = (self.__C.size,self.__C.size)
 560             self.size        = self.__C.size**2
 561         elif asCovariance is not None:
 562             self.__is_matrix = True
 563             self.__C         = numpy.matrix( asCovariance, float )
 564             self.shape       = self.__C.shape
 565             self.size        = self.__C.size
 566         elif asCovObject is not None:
 567             self.__is_object = True
 568             self.__C         = asCovObject
 569             for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
 570                 if not hasattr(self.__C,at):
 571                     raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
 572             if hasattr(self.__C,"shape"):
 573                 self.shape       = self.__C.shape
 574             else:
 575                 self.shape       = (0,0)
 576             if hasattr(self.__C,"size"):
 577                 self.size        = self.__C.size
 578             else:
 579                 self.size        = 0
 580         else:
 581             pass
 582             # 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)
 583         #
 584         self.__validate()
 585
 586     def __validate(self):
 587         "Validation"
 588         if self.ismatrix() and min(self.shape) != max(self.shape):
 589             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))
 590         if self.isobject() and min(self.shape) != max(self.shape):
 591             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))
 592         if self.isscalar() and self.__C <= 0:
 593             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
 594         if self.isvector() and (self.__C <= 0).any():
 595             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
 596         if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
 597             try:
 598                 L = numpy.linalg.cholesky( self.__C )
 599             except:
 600                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 601         if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
 602             try:
 603                 L = self.__C.cholesky()
 604             except:
 605                 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 606
 607     def isscalar(self):
 608         "Vérification du type interne"
 609         return self.__is_scalar
 610
 611     def isvector(self):
 612         "Vérification du type interne"
 613         return self.__is_vector
 614
 615     def ismatrix(self):
 616         "Vérification du type interne"
 617         return self.__is_matrix
 618
 619     def isobject(self):
 620         "Vérification du type interne"
 621         return self.__is_object
 622
 623     def getI(self):
 624         "Inversion"
 625         if   self.ismatrix():
 626             return Covariance(self.__name+"I", asCovariance  = self.__C.I )
 627         elif self.isvector():
 628             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
 629         elif self.isscalar():
 630             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
 631         elif self.isobject():
 632             return Covariance(self.__name+"I", asCovObject   = self.__C.getI() )
 633         else:
 634             return None # Indispensable
 635
 636     def getT(self):
 637         "Transposition"
 638         if   self.ismatrix():
 639             return Covariance(self.__name+"T", asCovariance  = self.__C.T )
 640         elif self.isvector():
 641             return Covariance(self.__name+"T", asEyeByVector = self.__C )
 642         elif self.isscalar():
 643             return Covariance(self.__name+"T", asEyeByScalar = self.__C )
 644         elif self.isobject():
 645             return Covariance(self.__name+"T", asCovObject   = self.__C.getT() )
 646
 647     def cholesky(self):
 648         "Décomposition de Cholesky"
 649         if   self.ismatrix():
 650             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__C) )
 651         elif self.isvector():
 652             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
 653         elif self.isscalar():
 654             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
 655         elif self.isobject() and hasattr(self.__C,"cholesky"):
 656             return Covariance(self.__name+"C", asCovObject   = self.__C.cholesky() )
 657
 658     def choleskyI(self):
 659         "Inversion de la décomposition de Cholesky"
 660         if   self.ismatrix():
 661             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__C).I )
 662         elif self.isvector():
 663             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
 664         elif self.isscalar():
 665             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
 666         elif self.isobject() and hasattr(self.__C,"choleskyI"):
 667             return Covariance(self.__name+"H", asCovObject   = self.__C.choleskyI() )
 668
 669     def diag(self, msize=None):
 670         "Diagonale de la matrice"
 671         if   self.ismatrix():
 672             return numpy.diag(self.__C)
 673         elif self.isvector():
 674             return self.__C
 675         elif self.isscalar():
 676             if msize is None:
 677                 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,))
 678             else:
 679                 return self.__C * numpy.ones(int(msize))
 680         elif self.isobject():
 681             return self.__C.diag()
 682
 683     def asfullmatrix(self, msize=None):
 684         "Matrice pleine"
 685         if   self.ismatrix():
 686             return self.__C
 687         elif self.isvector():
 688             return numpy.matrix( numpy.diag(self.__C), float )
 689         elif self.isscalar():
 690             if msize is None:
 691                 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,))
 692             else:
 693                 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
 694         elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
 695             return self.__C.asfullmatrix()
 696
 697     def trace(self, msize=None):
 698         "Trace de la matrice"
 699         if   self.ismatrix():
 700             return numpy.trace(self.__C)
 701         elif self.isvector():
 702             return float(numpy.sum(self.__C))
 703         elif self.isscalar():
 704             if msize is None:
 705                 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,))
 706             else:
 707                 return self.__C * int(msize)
 708         elif self.isobject():
 709             return self.__C.trace()
 710
 711     def __repr__(self):
 712         "x.__repr__() <==> repr(x)"
 713         return repr(self.__C)
 714
 715     def __str__(self):
 716         "x.__str__() <==> str(x)"
 717         return str(self.__C)
 718
 719     def __add__(self, other):
 720         "x.__add__(y) <==> x+y"
 721         if   self.ismatrix() or self.isobject():
 722             return self.__C + numpy.asmatrix(other)
 723         elif self.isvector() or self.isscalar():
 724             _A = numpy.asarray(other)
 725             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
 726             return numpy.asmatrix(_A)
 727
 728     def __radd__(self, other):
 729         "x.__radd__(y) <==> y+x"
 730         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
 731
 732     def __sub__(self, other):
 733         "x.__sub__(y) <==> x-y"
 734         if   self.ismatrix() or self.isobject():
 735             return self.__C - numpy.asmatrix(other)
 736         elif self.isvector() or self.isscalar():
 737             _A = numpy.asarray(other)
 738             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
 739             return numpy.asmatrix(_A)
 740
 741     def __rsub__(self, other):
 742         "x.__rsub__(y) <==> y-x"
 743         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
 744
 745     def __neg__(self):
 746         "x.__neg__() <==> -x"
 747         return - self.__C
 748
 749     def __mul__(self, other):
 750         "x.__mul__(y) <==> x*y"
 751         if   self.ismatrix() and isinstance(other,numpy.matrix):
 752             return self.__C * other
 753         elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
 754                                or isinstance(other,list) \
 755                                or isinstance(other,tuple)):
 756             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 757                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
 758             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 759                 return self.__C * numpy.asmatrix(other)
 760             else:
 761                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
 762         elif self.isvector() and (isinstance(other,numpy.matrix) \
 763                                or isinstance(other,numpy.ndarray) \
 764                                or isinstance(other,list) \
 765                                or isinstance(other,tuple)):
 766             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 767                 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
 768             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 769                 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
 770             else:
 771                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 772         elif self.isscalar() and isinstance(other,numpy.matrix):
 773             return self.__C * other
 774         elif self.isscalar() and (isinstance(other,numpy.ndarray) \
 775                                or isinstance(other,list) \
 776                                or isinstance(other,tuple)):
 777             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
 778                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
 779             else:
 780                 return self.__C * numpy.asmatrix(other)
 781         elif self.isobject():
 782             return self.__C.__mul__(other)
 783         else:
 784             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
 785
 786     def __rmul__(self, other):
 787         "x.__rmul__(y) <==> y*x"
 788         if self.ismatrix() and isinstance(other,numpy.matrix):
 789             return other * self.__C
 790         elif self.isvector() and isinstance(other,numpy.matrix):
 791             if numpy.ravel(other).size == self.shape[0]: # Vecteur
 792                 return numpy.asmatrix(numpy.ravel(other) * self.__C)
 793             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
 794                 return numpy.asmatrix(numpy.array(other) * self.__C)
 795             else:
 796                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 797         elif self.isscalar() and isinstance(other,numpy.matrix):
 798             return other * self.__C
 799         elif self.isobject():
 800             return self.__C.__rmul__(other)
 801         else:
 802             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
 803
 804     def __len__(self):
 805         "x.__len__() <==> len(x)"
 806         return self.shape[0]
 807
 808 # ==============================================================================
 809 def CostFunction3D(_x,
 810                    _Hm  = None,  # Pour simuler Hm(x) : HO["Direct"].appliedTo
 811                    _HmX = None,  # Simulation déjà faite de Hm(x)
 812                    _arg = None,  # Arguments supplementaires pour Hm, sous la forme d'un tuple
 813                    _BI  = None,
 814                    _RI  = None,
 815                    _Xb  = None,
 816                    _Y   = None,
 817                    _SIV = False, # A résorber pour la 8.0
 818                    _SSC = [],    # self._parameters["StoreSupplementaryCalculations"]
 819                    _nPS = 0,     # nbPreviousSteps
 820                    _QM  = "DA",  # QualityMeasure
 821                    _SSV = {},    # Entrée et/ou sortie : self.StoredVariables
 822                    _fRt = False, # Restitue ou pas la sortie étendue
 823                    _sSc = True,  # Stocke ou pas les SSC
 824                   ):
 825     """
 826     Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
 827     et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
 828     DFO, QuantileRegression
 829     """
 830     if not _sSc:
 831         _SIV = False
 832         _SSC = {}
 833     else:
 834         for k in ["CostFunctionJ",
 835                   "CostFunctionJb",
 836                   "CostFunctionJo",
 837                   "CurrentOptimum",
 838                   "CurrentState",
 839                   "IndexOfOptimum",
 840                   "SimulatedObservationAtCurrentOptimum",
 841                   "SimulatedObservationAtCurrentState",
 842                  ]:
 843             if k not in _SSV:
 844                 _SSV[k] = []
 845             if hasattr(_SSV[k],"store"):
 846                 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
 847     #
 848     _X  = numpy.asmatrix(numpy.ravel( _x )).T
 849     if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
 850         _SSV["CurrentState"].append( _X )
 851     #
 852     if _HmX is not None:
 853         _HX = _HmX
 854     else:
 855         if _Hm is None:
 856             raise ValueError("COSTFUNCTION3D Operator has to be defined.")
 857         if _arg is None:
 858             _HX = _Hm( _X )
 859         else:
 860             _HX = _Hm( _X, *_arg )
 861     _HX = numpy.asmatrix(numpy.ravel( _HX )).T
 862     #
 863     if "SimulatedObservationAtCurrentState" in _SSC or \
 864        "SimulatedObservationAtCurrentOptimum" in _SSC:
 865         _SSV["SimulatedObservationAtCurrentState"].append( _HX )
 866     #
 867     if numpy.any(numpy.isnan(_HX)):
 868         Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
 869     else:
 870         _Y   = numpy.asmatrix(numpy.ravel( _Y )).T
 871         if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
 872             if _BI is None or _RI is None:
 873                 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
 874             _Xb  = numpy.asmatrix(numpy.ravel( _Xb )).T
 875             Jb  = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
 876             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
 877         elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
 878             if _RI is None:
 879                 raise ValueError("Observation error covariance matrix has to be properly defined!")
 880             Jb  = 0.
 881             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
 882         elif _QM in ["LeastSquares", "LS", "L2"]:
 883             Jb  = 0.
 884             Jo  = 0.5 * (_Y - _HX).T * (_Y - _HX)
 885         elif _QM in ["AbsoluteValue", "L1"]:
 886             Jb  = 0.
 887             Jo  = numpy.sum( numpy.abs(_Y - _HX) )
 888         elif _QM in ["MaximumError", "ME"]:
 889             Jb  = 0.
 890             Jo  = numpy.max( numpy.abs(_Y - _HX) )
 891         elif _QM in ["QR", "Null"]:
 892             Jb  = 0.
 893             Jo  = 0.
 894         else:
 895             raise ValueError("Unknown asked quality measure!")
 896         #
 897         J   = float( Jb ) + float( Jo )
 898     #
 899     if _sSc:
 900         _SSV["CostFunctionJb"].append( Jb )
 901         _SSV["CostFunctionJo"].append( Jo )
 902         _SSV["CostFunctionJ" ].append( J )
 903     #
 904     if "IndexOfOptimum" in _SSC or \
 905        "CurrentOptimum" in _SSC or \
 906        "SimulatedObservationAtCurrentOptimum" in _SSC:
 907         IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
 908     if "IndexOfOptimum" in _SSC:
 909         _SSV["IndexOfOptimum"].append( IndexMin )
 910     if "CurrentOptimum" in _SSC:
 911         _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
 912     if "SimulatedObservationAtCurrentOptimum" in _SSC:
 913         _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
 914     #
 915     if _fRt:
 916         return _SSV
 917     else:
 918         if _QM in ["QR"]: # Pour le QuantileRegression
 919             return _HX
 920         else:
 921             return J
 922
 923 # ==============================================================================
 924 if __name__ == "__main__":
 925     print '\n AUTODIAGNOSTIC \n'