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, os
  32 import numpy
  33 from daCore import Persistence
  34 from daCore 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 range(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, paire ):
 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         assert len(paire) == 2, "Incorrect number of arguments"
 189         xValue, uValue = paire
 190         if self.__Matrix is not None:
 191             self.__addOneMatrixCall()
 192             return self.__Matrix * xValue
 193         elif uValue is not None:
 194             self.__addOneMethodCall()
 195             return self.__Method( (xValue, uValue) )
 196         else:
 197             self.__addOneMethodCall()
 198             return self.__Method( xValue )
 199
 200     def appliedInXTo(self, paire ):
 201         """
 202         Permet de restituer le résultat de l'application de l'opérateur à un
 203         argument xValue, sachant que l'opérateur est valable en xNominal.
 204         Cette méthode se contente d'appliquer, son argument devant a priori
 205         être du bon type. Si l'opérateur est linéaire car c'est une matrice,
 206         alors il est valable en tout point nominal et il n'est pas nécessaire
 207         d'utiliser xNominal.
 208         Arguments : une liste contenant
 209         - xNominal : argument permettant de donner le point où l'opérateur
 210           est construit pour etre ensuite appliqué
 211         - xValue : argument adapté pour appliquer l'opérateur
 212         """
 213         assert len(paire) == 2, "Incorrect number of arguments"
 214         xNominal, xValue = paire
 215         if self.__Matrix is not None:
 216             self.__addOneMatrixCall()
 217             return self.__Matrix * xValue
 218         else:
 219             self.__addOneMethodCall()
 220             return self.__Method( (xNominal, xValue) )
 221
 222     def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
 223         """
 224         Permet de renvoyer l'opérateur sous la forme d'une matrice
 225         """
 226         if self.__Matrix is not None:
 227             self.__addOneMatrixCall()
 228             return self.__Matrix
 229         elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
 230             self.__addOneMethodCall()
 231             return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
 232         else:
 233             raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
 234
 235     def shape(self):
 236         """
 237         Renvoie la taille sous forme numpy si l'opérateur est disponible sous
 238         la forme d'une matrice
 239         """
 240         if self.__Matrix is not None:
 241             return self.__Matrix.shape
 242         else:
 243             raise ValueError("Matrix form of the operator is not available, nor the shape")
 244
 245     def nbcalls(self, which=None):
 246         """
 247         Renvoie les nombres d'évaluations de l'opérateur
 248         """
 249         __nbcalls = (
 250             self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
 251             self.__NbCallsAsMatrix,
 252             self.__NbCallsAsMethod,
 253             self.__NbCallsOfCached,
 254             Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
 255             Operator.NbCallsAsMatrix,
 256             Operator.NbCallsAsMethod,
 257             Operator.NbCallsOfCached,
 258             )
 259         if which is None: return __nbcalls
 260         else:             return __nbcalls[which]
 261
 262     def __addOneMatrixCall(self):
 263         "Comptabilise un appel"
 264         self.__NbCallsAsMatrix   += 1 # Decompte local
 265         Operator.NbCallsAsMatrix += 1 # Decompte global
 266
 267     def __addOneMethodCall(self):
 268         "Comptabilise un appel"
 269         self.__NbCallsAsMethod   += 1 # Decompte local
 270         Operator.NbCallsAsMethod += 1 # Decompte global
 271
 272     def __addOneCacheCall(self):
 273         "Comptabilise un appel"
 274         self.__NbCallsOfCached   += 1 # Decompte local
 275         Operator.NbCallsOfCached += 1 # Decompte global
 276
 277 # ==============================================================================
 278 class Algorithm(object):
 279     """
 280     Classe générale d'interface de type algorithme
 281
 282     Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
 283     d'assimilation, en fournissant un container (dictionnaire) de variables
 284     persistantes initialisées, et des méthodes d'accès à ces variables stockées.
 285
 286     Une classe élémentaire d'algorithme doit implémenter la méthode "run".
 287     """
 288     def __init__(self, name):
 289         """
 290         L'initialisation présente permet de fabriquer des variables de stockage
 291         disponibles de manière générique dans les algorithmes élémentaires. Ces
 292         variables de stockage sont ensuite conservées dans un dictionnaire
 293         interne à l'objet, mais auquel on accède par la méthode "get".
 294
 295         Les variables prévues sont :
 296             - CostFunctionJ  : fonction-cout globale, somme des deux parties suivantes
 297             - CostFunctionJb : partie ébauche ou background de la fonction-cout
 298             - CostFunctionJo : partie observations de la fonction-cout
 299             - GradientOfCostFunctionJ  : gradient de la fonction-cout globale
 300             - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
 301             - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
 302             - CurrentState : état courant lors d'itérations
 303             - Analysis : l'analyse Xa
 304             - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
 305             - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
 306             - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
 307             - Innovation : l'innovation : d = Y - H(X)
 308             - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
 309             - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
 310             - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
 311             - MahalanobisConsistency : indicateur de consistance des covariances
 312             - OMA : Observation moins Analysis : Y - Xa
 313             - OMB : Observation moins Background : Y - Xb
 314             - AMB : Analysis moins Background : Xa - Xb
 315             - APosterioriCovariance : matrice A
 316             - APosterioriVariances : variances de la matrice A
 317             - APosterioriStandardDeviations : écart-types de la matrice A
 318             - APosterioriCorrelations : correlations de la matrice A
 319             - Residu : dans le cas des algorithmes de vérification
 320         On peut rajouter des variables à stocker dans l'initialisation de
 321         l'algorithme élémentaire qui va hériter de cette classe
 322         """
 323         logging.debug("%s Initialisation", str(name))
 324         self._m = PlatformInfo.SystemUsage()
 325         #
 326         self._name = str( name )
 327         self._parameters = {"StoreSupplementaryCalculations":[]}
 328         self.__required_parameters = {}
 329         self.StoredVariables = {}
 330         #
 331         self.StoredVariables["CostFunctionJ"]                        = Persistence.OneScalar(name = "CostFunctionJ")
 332         self.StoredVariables["CostFunctionJb"]                       = Persistence.OneScalar(name = "CostFunctionJb")
 333         self.StoredVariables["CostFunctionJo"]                       = Persistence.OneScalar(name = "CostFunctionJo")
 334         self.StoredVariables["CostFunctionJAtCurrentOptimum"]        = Persistence.OneScalar(name = "CostFunctionJAtCurrentOptimum")
 335         self.StoredVariables["CostFunctionJbAtCurrentOptimum"]       = Persistence.OneScalar(name = "CostFunctionJbAtCurrentOptimum")
 336         self.StoredVariables["CostFunctionJoAtCurrentOptimum"]       = Persistence.OneScalar(name = "CostFunctionJoAtCurrentOptimum")
 337         self.StoredVariables["GradientOfCostFunctionJ"]              = Persistence.OneVector(name = "GradientOfCostFunctionJ")
 338         self.StoredVariables["GradientOfCostFunctionJb"]             = Persistence.OneVector(name = "GradientOfCostFunctionJb")
 339         self.StoredVariables["GradientOfCostFunctionJo"]             = Persistence.OneVector(name = "GradientOfCostFunctionJo")
 340         self.StoredVariables["CurrentState"]                         = Persistence.OneVector(name = "CurrentState")
 341         self.StoredVariables["Analysis"]                             = Persistence.OneVector(name = "Analysis")
 342         self.StoredVariables["IndexOfOptimum"]                       = Persistence.OneIndex(name = "IndexOfOptimum")
 343         self.StoredVariables["CurrentOptimum"]                       = Persistence.OneVector(name = "CurrentOptimum")
 344         self.StoredVariables["SimulatedObservationAtBackground"]     = Persistence.OneVector(name = "SimulatedObservationAtBackground")
 345         self.StoredVariables["SimulatedObservationAtCurrentState"]   = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
 346         self.StoredVariables["SimulatedObservationAtOptimum"]        = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
 347         self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
 348         self.StoredVariables["Innovation"]                           = Persistence.OneVector(name = "Innovation")
 349         self.StoredVariables["InnovationAtCurrentState"]             = Persistence.OneVector(name = "InnovationAtCurrentState")
 350         self.StoredVariables["SigmaObs2"]                            = Persistence.OneScalar(name = "SigmaObs2")
 351         self.StoredVariables["SigmaBck2"]                            = Persistence.OneScalar(name = "SigmaBck2")
 352         self.StoredVariables["MahalanobisConsistency"]               = Persistence.OneScalar(name = "MahalanobisConsistency")
 353         self.StoredVariables["OMA"]                                  = Persistence.OneVector(name = "OMA")
 354         self.StoredVariables["OMB"]                                  = Persistence.OneVector(name = "OMB")
 355         self.StoredVariables["BMA"]                                  = Persistence.OneVector(name = "BMA")
 356         self.StoredVariables["APosterioriCovariance"]                = Persistence.OneMatrix(name = "APosterioriCovariance")
 357         self.StoredVariables["APosterioriVariances"]                 = Persistence.OneVector(name = "APosterioriVariances")
 358         self.StoredVariables["APosterioriStandardDeviations"]        = Persistence.OneVector(name = "APosterioriStandardDeviations")
 359         self.StoredVariables["APosterioriCorrelations"]              = Persistence.OneMatrix(name = "APosterioriCorrelations")
 360         self.StoredVariables["SimulationQuantiles"]                  = Persistence.OneMatrix(name = "SimulationQuantiles")
 361         self.StoredVariables["Residu"]                               = Persistence.OneScalar(name = "Residu")
 362
 363     def _pre_run(self, Parameters ):
 364         "Pré-calcul"
 365         logging.debug("%s Lancement", self._name)
 366         logging.debug("%s Taille mémoire utilisée de %.0f Mio", self._name, self._m.getUsedMemory("Mio"))
 367         #
 368         # Mise a jour de self._parameters avec Parameters
 369         self.__setParameters(Parameters)
 370         #
 371         # Corrections et complements
 372         if "Bounds" in self._parameters and (type(self._parameters["Bounds"]) is type([]) or type(self._parameters["Bounds"]) is type(())) and (len(self._parameters["Bounds"]) > 0):
 373             logging.debug("%s Prise en compte des bornes effectuee"%(self._name,))
 374         else:
 375             self._parameters["Bounds"] = None
 376         #
 377         if logging.getLogger().level < logging.WARNING:
 378             self._parameters["optiprint"], self._parameters["optdisp"] = 1, 1
 379             if PlatformInfo.has_scipy:
 380                 import scipy.optimize
 381                 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_ALL
 382             else:
 383                 self._parameters["optmessages"] = 15
 384         else:
 385             self._parameters["optiprint"], self._parameters["optdisp"] = -1, 0
 386             if PlatformInfo.has_scipy:
 387                 import scipy.optimize
 388                 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_NONE
 389             else:
 390                 self._parameters["optmessages"] = 15
 391         #
 392         return 0
 393
 394     def _post_run(self,_oH=None):
 395         "Post-calcul"
 396         if ("StoreSupplementaryCalculations" in self._parameters) and \
 397             "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
 398             for _A in self.StoredVariables["APosterioriCovariance"]:
 399                 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
 400                     self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
 401                 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
 402                     self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
 403                 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
 404                     _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
 405                     _C = numpy.dot(_EI, numpy.dot(_A, _EI))
 406                     self.StoredVariables["APosterioriCorrelations"].store( _C )
 407         if _oH is not None:
 408             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))
 409             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))
 410         logging.debug("%s Taille mémoire utilisée de %.0f Mio", self._name, self._m.getUsedMemory("Mio"))
 411         logging.debug("%s Terminé", self._name)
 412         return 0
 413
 414     def get(self, key=None):
 415         """
 416         Renvoie l'une des variables stockées identifiée par la clé, ou le
 417         dictionnaire de l'ensemble des variables disponibles en l'absence de
 418         clé. Ce sont directement les variables sous forme objet qui sont
 419         renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
 420         des classes de persistance.
 421         """
 422         if key is not None:
 423             return self.StoredVariables[key]
 424         else:
 425             return self.StoredVariables
 426
 427     def __contains__(self, key=None):
 428         "D.__contains__(k) -> True if D has a key k, else False"
 429         return key in self.StoredVariables
 430
 431     def keys(self):
 432         "D.keys() -> list of D's keys"
 433         return self.StoredVariables.keys()
 434
 435     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
 436         """
 437         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
 438         sa forme mathématique la plus naturelle possible.
 439         """
 440         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
 441
 442     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
 443         """
 444         Permet de définir dans l'algorithme des paramètres requis et leurs
 445         caractéristiques par défaut.
 446         """
 447         if name is None:
 448             raise ValueError("A name is mandatory to define a required parameter.")
 449         #
 450         self.__required_parameters[name] = {
 451             "default"  : default,
 452             "typecast" : typecast,
 453             "minval"   : minval,
 454             "maxval"   : maxval,
 455             "listval"  : listval,
 456             "message"  : message,
 457             }
 458         logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
 459
 460     def getRequiredParameters(self, noDetails=True):
 461         """
 462         Renvoie la liste des noms de paramètres requis ou directement le
 463         dictionnaire des paramètres requis.
 464         """
 465         if noDetails:
 466             ks = list(self.__required_parameters.keys())
 467             ks.sort()
 468             return ks
 469         else:
 470             return self.__required_parameters
 471
 472     def setParameterValue(self, name=None, value=None):
 473         """
 474         Renvoie la valeur d'un paramètre requis de manière contrôlée
 475         """
 476         default  = self.__required_parameters[name]["default"]
 477         typecast = self.__required_parameters[name]["typecast"]
 478         minval   = self.__required_parameters[name]["minval"]
 479         maxval   = self.__required_parameters[name]["maxval"]
 480         listval  = self.__required_parameters[name]["listval"]
 481         #
 482         if value is None and default is None:
 483             __val = None
 484         elif value is None and default is not None:
 485             if typecast is None: __val = default
 486             else:                __val = typecast( default )
 487         else:
 488             if typecast is None: __val = value
 489             else:                __val = typecast( value )
 490         #
 491         if minval is not None and (numpy.array(__val, float) < minval).any():
 492             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
 493         if maxval is not None and (numpy.array(__val, float) > maxval).any():
 494             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
 495         if listval is not None:
 496             if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
 497                 for v in __val:
 498                     if v not in listval:
 499                         raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
 500             elif __val not in listval:
 501                 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
 502         return __val
 503
 504     def __setParameters(self, fromDico={}):
 505         """
 506         Permet de stocker les paramètres reçus dans le dictionnaire interne.
 507         """
 508         self._parameters.update( fromDico )
 509         for k in self.__required_parameters.keys():
 510             if k in fromDico.keys():
 511                 self._parameters[k] = self.setParameterValue(k,fromDico[k])
 512             else:
 513                 self._parameters[k] = self.setParameterValue(k)
 514             logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
 515
 516 # ==============================================================================
 517 class Diagnostic(object):
 518     """
 519     Classe générale d'interface de type diagnostic
 520
 521     Ce template s'utilise de la manière suivante : il sert de classe "patron" en
 522     même temps que l'une des classes de persistance, comme "OneScalar" par
 523     exemple.
 524
 525     Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
 526     méthode "_formula" pour écrire explicitement et proprement la formule pour
 527     l'écriture mathématique du calcul du diagnostic (méthode interne non
 528     publique), et "calculate" pour activer la précédente tout en ayant vérifié
 529     et préparé les données, et pour stocker les résultats à chaque pas (méthode
 530     externe d'activation).
 531     """
 532     def __init__(self, name = "", parameters = {}):
 533         "Initialisation"
 534         self.name       = str(name)
 535         self.parameters = dict( parameters )
 536
 537     def _formula(self, *args):
 538         """
 539         Doit implémenter l'opération élémentaire de diagnostic sous sa forme
 540         mathématique la plus naturelle possible.
 541         """
 542         raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
 543
 544     def calculate(self, *args):
 545         """
 546         Active la formule de calcul avec les arguments correctement rangés
 547         """
 548         raise NotImplementedError("Diagnostic activation method has not been implemented!")
 549
 550 # ==============================================================================
 551 class ParameterDictionary(object):
 552     """
 553     Classe générale d'interface de type dictionnaire de paramètres
 554     """
 555     def __init__(self,
 556                  name               = "GenericParamDict",
 557                  asAlgorithm        = None,
 558                  asDict             = None,
 559                  asScript           = None,
 560                 ):
 561         """
 562         """
 563         self.__name       = str(name)
 564         self.__A          = None
 565         self.__D          = None
 566         #
 567         if asScript is not None:
 568             __Algo = ImportFromScript(asScript).getvalue( "Algorithm" )
 569             __Dict = ImportFromScript(asScript).getvalue( self.__name, "Parameters" )
 570         else:
 571             __Algo = asAlgorithm
 572             __Dict = asDict
 573         #
 574         if __Algo is not None:
 575             self.__A = str(__Algo)
 576         if __Dict is not None:
 577             self.__D = dict(__Dict)
 578
 579     def __repr__(self):
 580         "x.__repr__() <==> repr(x)"
 581         return repr(self.__A)+"\n"+repr(self.__D)
 582
 583     def __str__(self):
 584         "x.__str__() <==> str(x)"
 585         return str(self.__A)+"\n"+str(self.__D)
 586
 587 # ==============================================================================
 588 class DataObserver(object):
 589     """
 590     Classe générale d'interface de type observer
 591     """
 592     def __init__(self,
 593                  name       = "GenericObserver",
 594                  onVariable = None,
 595                  asTemplate = None,
 596                  asString   = None,
 597                  asScript   = None,
 598                  withInfo   = None,
 599                 ):
 600         self.__name       = str(name)
 601         self.__V          = None
 602         self.__O          = None
 603         self.__I          = None
 604         #
 605         if onVariable is None:
 606             raise ValueError("setting an observer has to be done over a variable name, not over None.")
 607         else:
 608             self.__V = str(onVariable)
 609         #
 610         if withInfo is None:
 611             self.__I = str(onVariable)
 612         else:
 613             self.__I = str(withInfo)
 614         #
 615         if asString is not None:
 616             __FunctionText = asString
 617         elif (asTemplate is not None) and (asTemplate in ObserverTemplates):
 618             __FunctionText = ObserverTemplates[asTemplate]
 619         elif asScript is not None:
 620             __FunctionText = ImportFromScript(asScript).getstring()
 621         else:
 622             __FunctionText = ""
 623         __Function = ObserverF(__FunctionText)
 624         #
 625         self.__O = __Function.getfunc()
 626         #
 627         return {self.__V:(self.__O, self.__I)}
 628
 629     def __repr__(self):
 630         "x.__repr__() <==> repr(x)"
 631         return repr(self.__V)+"\n"+repr(self.__O)
 632
 633     def __str__(self):
 634         "x.__str__() <==> str(x)"
 635         return str(self.__V)+"\n"+str(self.__O)
 636
 637 # ==============================================================================
 638 class State(object):
 639     """
 640     Classe générale d'interface de type état
 641     """
 642     def __init__(self,
 643                  name               = "GenericVector",
 644                  asVector           = None,
 645                  asPersistentVector = None,
 646                  asScript           = None,
 647                  Scheduler          = None,
 648                  toBeChecked        = False,
 649                 ):
 650         """
 651         Permet de définir un vecteur :
 652         - asVector : entrée des données, comme un vecteur compatible avec le
 653           constructeur de numpy.matrix, ou "True" si entrée par script.
 654         - asPersistentVector : entrée des données, comme une série de vecteurs
 655           compatible avec le constructeur de numpy.matrix, ou comme un objet de
 656           type Persistence, ou "True" si entrée par script.
 657         - asScript : si un script valide est donné contenant une variable
 658           nommée "name", la variable est de type "asVector" (par défaut) ou
 659           "asPersistentVector" selon que l'une de ces variables est placée à
 660           "True".
 661         """
 662         self.__name       = str(name)
 663         self.__check      = bool(toBeChecked)
 664         #
 665         self.__V          = None
 666         self.__T          = None
 667         self.__is_vector  = False
 668         self.__is_series  = False
 669         #
 670         if asScript is not None:
 671             __Vector, __Series = None, None
 672             if asPersistentVector:
 673                 __Series = ImportFromScript(asScript).getvalue( self.__name )
 674             else:
 675                 __Vector = ImportFromScript(asScript).getvalue( self.__name )
 676         else:
 677             __Vector, __Series = asVector, asPersistentVector
 678         #
 679         if __Vector is not None:
 680             self.__is_vector = True
 681             self.__V         = numpy.matrix( numpy.ravel(numpy.matrix(asVector)), numpy.float ).T
 682             self.shape       = self.__V.shape
 683             self.size        = self.__V.size
 684         elif __Series is not None:
 685             self.__is_series  = True
 686             if type(__Series) in [tuple, list, numpy.ndarray, numpy.matrix]:
 687                 self.__V = Persistence.OneVector(self.__name, basetype=numpy.matrix)
 688                 for member in __Series:
 689                     self.__V.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
 690                 import sys ; sys.stdout.flush()
 691             else:
 692                 self.__V = __Series
 693             if type(self.__V.shape) in (tuple, list):
 694                 self.shape       = self.__V.shape
 695             else:
 696                 self.shape       = self.__V.shape()
 697             self.size        = self.shape[0] * self.shape[1]
 698         else:
 699             raise ValueError("The %s object is improperly defined, it requires at minima either a vector, a list/tuple of vectors or a persistent object. Please check your vector input."%self.__name)
 700         #
 701         if Scheduler is not None:
 702             self.__T = Scheduler
 703
 704     def isvector(self):
 705         "Vérification du type interne"
 706         return self.__is_vector
 707
 708     def isseries(self):
 709         "Vérification du type interne"
 710         return self.__is_series
 711
 712     def __repr__(self):
 713         "x.__repr__() <==> repr(x)"
 714         return repr(self.__V)
 715
 716     def __str__(self):
 717         "x.__str__() <==> str(x)"
 718         return str(self.__V)
 719
 720 # ==============================================================================
 721 class Covariance(object):
 722     """
 723     Classe générale d'interface de type covariance
 724     """
 725     def __init__(self,
 726                  name          = "GenericCovariance",
 727                  asCovariance  = None,
 728                  asEyeByScalar = None,
 729                  asEyeByVector = None,
 730                  asCovObject   = None,
 731                  asScript      = None,
 732                  toBeChecked   = False,
 733                 ):
 734         """
 735         Permet de définir une covariance :
 736         - asCovariance : entrée des données, comme une matrice compatible avec
 737           le constructeur de numpy.matrix
 738         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 739           multiplicatif d'une matrice de corrélation identité, aucune matrice
 740           n'étant donc explicitement à donner
 741         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 742           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 743           n'étant donc explicitement à donner
 744         - asCovObject : entrée des données comme un objet python, qui a les
 745           methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
 746           "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
 747           "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
 748         - toBeChecked : booléen indiquant si le caractère SDP de la matrice
 749           pleine doit être vérifié
 750         """
 751         self.__name       = str(name)
 752         self.__check      = bool(toBeChecked)
 753         #
 754         self.__C          = None
 755         self.__is_scalar  = False
 756         self.__is_vector  = False
 757         self.__is_matrix  = False
 758         self.__is_object  = False
 759         #
 760         if asScript is not None:
 761             __Matrix, __Scalar, __Vector, __Object = None, None, None, None
 762             if asEyeByScalar:
 763                 __Scalar = _ImportFromScript(asScript).getvalue( "BackgroundError" )
 764             elif asEyeByVector:
 765                 __Vector = _ImportFromScript(asScript).getvalue( "BackgroundError" )
 766             elif asCovObject:
 767                 __Object = _ImportFromScript(asScript).getvalue( "BackgroundError" )
 768             else:
 769                 __Matrix = _ImportFromScript(asScript).getvalue( "BackgroundError" )
 770         else:
 771             __Matrix, __Scalar, __Vector, __Object = asCovariance, asEyeByScalar, asEyeByVector, asCovObject
 772         #
 773         if __Scalar is not None:
 774             if numpy.matrix(__Scalar).size != 1:
 775                 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(__Scalar).size)
 776             self.__is_scalar = True
 777             self.__C         = numpy.abs( float(__Scalar) )
 778             self.shape       = (0,0)
 779             self.size        = 0
 780         elif __Vector is not None:
 781             self.__is_vector = True
 782             self.__C         = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(__Vector, float ) ) ) )
 783             self.shape       = (self.__C.size,self.__C.size)
 784             self.size        = self.__C.size**2
 785         elif __Matrix is not None:
 786             self.__is_matrix = True
 787             self.__C         = numpy.matrix( __Matrix, float )
 788             self.shape       = self.__C.shape
 789             self.size        = self.__C.size
 790         elif __Object is not None:
 791             self.__is_object = True
 792             self.__C         = __Object
 793             for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
 794                 if not hasattr(self.__C,at):
 795                     raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
 796             if hasattr(self.__C,"shape"):
 797                 self.shape       = self.__C.shape
 798             else:
 799                 self.shape       = (0,0)
 800             if hasattr(self.__C,"size"):
 801                 self.size        = self.__C.size
 802             else:
 803                 self.size        = 0
 804         else:
 805             pass
 806             # 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)
 807         #
 808         self.__validate()
 809
 810     def __validate(self):
 811         "Validation"
 812         if self.ismatrix() and min(self.shape) != max(self.shape):
 813             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))
 814         if self.isobject() and min(self.shape) != max(self.shape):
 815             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))
 816         if self.isscalar() and self.__C <= 0:
 817             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
 818         if self.isvector() and (self.__C <= 0).any():
 819             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
 820         if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
 821             try:
 822                 L = numpy.linalg.cholesky( self.__C )
 823             except:
 824                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 825         if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
 826             try:
 827                 L = self.__C.cholesky()
 828             except:
 829                 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
 830
 831     def isscalar(self):
 832         "Vérification du type interne"
 833         return self.__is_scalar
 834
 835     def isvector(self):
 836         "Vérification du type interne"
 837         return self.__is_vector
 838
 839     def ismatrix(self):
 840         "Vérification du type interne"
 841         return self.__is_matrix
 842
 843     def isobject(self):
 844         "Vérification du type interne"
 845         return self.__is_object
 846
 847     def getI(self):
 848         "Inversion"
 849         if   self.ismatrix():
 850             return Covariance(self.__name+"I", asCovariance  = self.__C.I )
 851         elif self.isvector():
 852             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
 853         elif self.isscalar():
 854             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
 855         elif self.isobject():
 856             return Covariance(self.__name+"I", asCovObject   = self.__C.getI() )
 857         else:
 858             return None # Indispensable
 859
 860     def getT(self):
 861         "Transposition"
 862         if   self.ismatrix():
 863             return Covariance(self.__name+"T", asCovariance  = self.__C.T )
 864         elif self.isvector():
 865             return Covariance(self.__name+"T", asEyeByVector = self.__C )
 866         elif self.isscalar():
 867             return Covariance(self.__name+"T", asEyeByScalar = self.__C )
 868         elif self.isobject():
 869             return Covariance(self.__name+"T", asCovObject   = self.__C.getT() )
 870
 871     def cholesky(self):
 872         "Décomposition de Cholesky"
 873         if   self.ismatrix():
 874             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__C) )
 875         elif self.isvector():
 876             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
 877         elif self.isscalar():
 878             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
 879         elif self.isobject() and hasattr(self.__C,"cholesky"):
 880             return Covariance(self.__name+"C", asCovObject   = self.__C.cholesky() )
 881
 882     def choleskyI(self):
 883         "Inversion de la décomposition de Cholesky"
 884         if   self.ismatrix():
 885             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__C).I )
 886         elif self.isvector():
 887             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
 888         elif self.isscalar():
 889             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
 890         elif self.isobject() and hasattr(self.__C,"choleskyI"):
 891             return Covariance(self.__name+"H", asCovObject   = self.__C.choleskyI() )
 892
 893     def diag(self, msize=None):
 894         "Diagonale de la matrice"
 895         if   self.ismatrix():
 896             return numpy.diag(self.__C)
 897         elif self.isvector():
 898             return self.__C
 899         elif self.isscalar():
 900             if msize is None:
 901                 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,))
 902             else:
 903                 return self.__C * numpy.ones(int(msize))
 904         elif self.isobject():
 905             return self.__C.diag()
 906
 907     def asfullmatrix(self, msize=None):
 908         "Matrice pleine"
 909         if   self.ismatrix():
 910             return self.__C
 911         elif self.isvector():
 912             return numpy.matrix( numpy.diag(self.__C), float )
 913         elif self.isscalar():
 914             if msize is None:
 915                 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,))
 916             else:
 917                 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
 918         elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
 919             return self.__C.asfullmatrix()
 920
 921     def trace(self, msize=None):
 922         "Trace de la matrice"
 923         if   self.ismatrix():
 924             return numpy.trace(self.__C)
 925         elif self.isvector():
 926             return float(numpy.sum(self.__C))
 927         elif self.isscalar():
 928             if msize is None:
 929                 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,))
 930             else:
 931                 return self.__C * int(msize)
 932         elif self.isobject():
 933             return self.__C.trace()
 934
 935     def __repr__(self):
 936         "x.__repr__() <==> repr(x)"
 937         return repr(self.__C)
 938
 939     def __str__(self):
 940         "x.__str__() <==> str(x)"
 941         return str(self.__C)
 942
 943     def __add__(self, other):
 944         "x.__add__(y) <==> x+y"
 945         if   self.ismatrix() or self.isobject():
 946             return self.__C + numpy.asmatrix(other)
 947         elif self.isvector() or self.isscalar():
 948             _A = numpy.asarray(other)
 949             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
 950             return numpy.asmatrix(_A)
 951
 952     def __radd__(self, other):
 953         "x.__radd__(y) <==> y+x"
 954         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
 955
 956     def __sub__(self, other):
 957         "x.__sub__(y) <==> x-y"
 958         if   self.ismatrix() or self.isobject():
 959             return self.__C - numpy.asmatrix(other)
 960         elif self.isvector() or self.isscalar():
 961             _A = numpy.asarray(other)
 962             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
 963             return numpy.asmatrix(_A)
 964
 965     def __rsub__(self, other):
 966         "x.__rsub__(y) <==> y-x"
 967         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
 968
 969     def __neg__(self):
 970         "x.__neg__() <==> -x"
 971         return - self.__C
 972
 973     def __mul__(self, other):
 974         "x.__mul__(y) <==> x*y"
 975         if   self.ismatrix() and isinstance(other,numpy.matrix):
 976             return self.__C * other
 977         elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
 978                                or isinstance(other,list) \
 979                                or isinstance(other,tuple)):
 980             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 981                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
 982             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 983                 return self.__C * numpy.asmatrix(other)
 984             else:
 985                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
 986         elif self.isvector() and (isinstance(other,numpy.matrix) \
 987                                or isinstance(other,numpy.ndarray) \
 988                                or isinstance(other,list) \
 989                                or isinstance(other,tuple)):
 990             if numpy.ravel(other).size == self.shape[1]: # Vecteur
 991                 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
 992             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
 993                 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
 994             else:
 995                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
 996         elif self.isscalar() and isinstance(other,numpy.matrix):
 997             return self.__C * other
 998         elif self.isscalar() and (isinstance(other,numpy.ndarray) \
 999                                or isinstance(other,list) \
1000                                or isinstance(other,tuple)):
1001             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
1002                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
1003             else:
1004                 return self.__C * numpy.asmatrix(other)
1005         elif self.isobject():
1006             return self.__C.__mul__(other)
1007         else:
1008             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
1009
1010     def __rmul__(self, other):
1011         "x.__rmul__(y) <==> y*x"
1012         if self.ismatrix() and isinstance(other,numpy.matrix):
1013             return other * self.__C
1014         elif self.isvector() and isinstance(other,numpy.matrix):
1015             if numpy.ravel(other).size == self.shape[0]: # Vecteur
1016                 return numpy.asmatrix(numpy.ravel(other) * self.__C)
1017             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
1018                 return numpy.asmatrix(numpy.array(other) * self.__C)
1019             else:
1020                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
1021         elif self.isscalar() and isinstance(other,numpy.matrix):
1022             return other * self.__C
1023         elif self.isobject():
1024             return self.__C.__rmul__(other)
1025         else:
1026             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
1027
1028     def __len__(self):
1029         "x.__len__() <==> len(x)"
1030         return self.shape[0]
1031
1032 # ==============================================================================
1033 class ObserverF(object):
1034     """
1035     Creation d'une fonction d'observateur a partir de son texte
1036     """
1037     def __init__(self, corps=""):
1038         self.__corps = corps
1039     def func(self,var,info):
1040         "Fonction d'observation"
1041         exec(self.__corps)
1042     def getfunc(self):
1043         "Restitution du pointeur de fonction dans l'objet"
1044         return self.func
1045
1046 # ==============================================================================
1047 class ImportFromScript(object):
1048     """
1049     Obtention d'une variable nommee depuis un fichier script importe
1050     """
1051     def __init__(self, __filename=None):
1052         "Verifie l'existence et importe le script"
1053         __filename = __filename.rstrip(".py")
1054         if __filename is None:
1055             raise ValueError("The name of the file containing the variable to be imported has to be specified.")
1056         if not os.path.isfile(str(__filename)+".py"):
1057             raise ValueError("The file containing the variable to be imported doesn't seem to exist. The given file name is:\n  \"%s\""%__filename)
1058         self.__scriptfile = __import__(__filename, globals(), locals(), [])
1059         self.__scriptstring = open(__filename+".py",'r').read()
1060     def getvalue(self, __varname=None, __synonym=None ):
1061         "Renvoie la variable demandee"
1062         if __varname is None:
1063             raise ValueError("The name of the variable to be imported has to be specified.")
1064         if not hasattr(self.__scriptfile, __varname):
1065             if __synonym is None:
1066                 raise ValueError("The imported script file doesn't contain the specified variable \"%s\"."%__varname)
1067             elif not hasattr(self.__scriptfile, __synonym):
1068                 raise ValueError("The imported script file doesn't contain the specified variable \"%s\"."%__synonym)
1069             else:
1070                 return getattr(self.__scriptfile, __synonym)
1071         else:
1072             return getattr(self.__scriptfile, __varname)
1073     def getstring(self):
1074         "Renvoie le script complet"
1075         return self.__scriptstring
1076
1077 # ==============================================================================
1078 def CostFunction3D(_x,
1079                    _Hm  = None,  # Pour simuler Hm(x) : HO["Direct"].appliedTo
1080                    _HmX = None,  # Simulation déjà faite de Hm(x)
1081                    _arg = None,  # Arguments supplementaires pour Hm, sous la forme d'un tuple
1082                    _BI  = None,
1083                    _RI  = None,
1084                    _Xb  = None,
1085                    _Y   = None,
1086                    _SIV = False, # A résorber pour la 8.0
1087                    _SSC = [],    # self._parameters["StoreSupplementaryCalculations"]
1088                    _nPS = 0,     # nbPreviousSteps
1089                    _QM  = "DA",  # QualityMeasure
1090                    _SSV = {},    # Entrée et/ou sortie : self.StoredVariables
1091                    _fRt = False, # Restitue ou pas la sortie étendue
1092                    _sSc = True,  # Stocke ou pas les SSC
1093                   ):
1094     """
1095     Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
1096     et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
1097     DFO, QuantileRegression
1098     """
1099     if not _sSc:
1100         _SIV = False
1101         _SSC = {}
1102     else:
1103         for k in ["CostFunctionJ",
1104                   "CostFunctionJb",
1105                   "CostFunctionJo",
1106                   "CurrentOptimum",
1107                   "CurrentState",
1108                   "IndexOfOptimum",
1109                   "SimulatedObservationAtCurrentOptimum",
1110                   "SimulatedObservationAtCurrentState",
1111                  ]:
1112             if k not in _SSV:
1113                 _SSV[k] = []
1114             if hasattr(_SSV[k],"store"):
1115                 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
1116     #
1117     _X  = numpy.asmatrix(numpy.ravel( _x )).T
1118     if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
1119         _SSV["CurrentState"].append( _X )
1120     #
1121     if _HmX is not None:
1122         _HX = _HmX
1123     else:
1124         if _Hm is None:
1125             raise ValueError("COSTFUNCTION3D Operator has to be defined.")
1126         if _arg is None:
1127             _HX = _Hm( _X )
1128         else:
1129             _HX = _Hm( _X, *_arg )
1130     _HX = numpy.asmatrix(numpy.ravel( _HX )).T
1131     #
1132     if "SimulatedObservationAtCurrentState" in _SSC or \
1133        "SimulatedObservationAtCurrentOptimum" in _SSC:
1134         _SSV["SimulatedObservationAtCurrentState"].append( _HX )
1135     #
1136     if numpy.any(numpy.isnan(_HX)):
1137         Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
1138     else:
1139         _Y   = numpy.asmatrix(numpy.ravel( _Y )).T
1140         if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
1141             if _BI is None or _RI is None:
1142                 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
1143             _Xb  = numpy.asmatrix(numpy.ravel( _Xb )).T
1144             Jb  = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
1145             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
1146         elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
1147             if _RI is None:
1148                 raise ValueError("Observation error covariance matrix has to be properly defined!")
1149             Jb  = 0.
1150             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
1151         elif _QM in ["LeastSquares", "LS", "L2"]:
1152             Jb  = 0.
1153             Jo  = 0.5 * (_Y - _HX).T * (_Y - _HX)
1154         elif _QM in ["AbsoluteValue", "L1"]:
1155             Jb  = 0.
1156             Jo  = numpy.sum( numpy.abs(_Y - _HX) )
1157         elif _QM in ["MaximumError", "ME"]:
1158             Jb  = 0.
1159             Jo  = numpy.max( numpy.abs(_Y - _HX) )
1160         elif _QM in ["QR", "Null"]:
1161             Jb  = 0.
1162             Jo  = 0.
1163         else:
1164             raise ValueError("Unknown asked quality measure!")
1165         #
1166         J   = float( Jb ) + float( Jo )
1167     #
1168     if _sSc:
1169         _SSV["CostFunctionJb"].append( Jb )
1170         _SSV["CostFunctionJo"].append( Jo )
1171         _SSV["CostFunctionJ" ].append( J )
1172     #
1173     if "IndexOfOptimum" in _SSC or \
1174        "CurrentOptimum" in _SSC or \
1175        "SimulatedObservationAtCurrentOptimum" in _SSC:
1176         IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
1177     if "IndexOfOptimum" in _SSC:
1178         _SSV["IndexOfOptimum"].append( IndexMin )
1179     if "CurrentOptimum" in _SSC:
1180         _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
1181     if "SimulatedObservationAtCurrentOptimum" in _SSC:
1182         _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
1183     #
1184     if _fRt:
1185         return _SSV
1186     else:
1187         if _QM in ["QR"]: # Pour le QuantileRegression
1188             return _HX
1189         else:
1190             return J
1191
1192 # ==============================================================================
1193 if __name__ == "__main__":
1194     print('\n AUTODIAGNOSTIC \n')