src/daComposant/daCore/BasicObjects.py

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