src/daComposant/daCore/BasicObjects.py

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