src/daComposant/daCore/AssimilationStudy.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2012 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 __doc__ = """
  24     Définit les outils généraux élémentaires.
  25
  26     Ce module est destiné à etre appelée par AssimilationStudy pour constituer
  27     les objets élémentaires de l'étude.
  28 """
  29 __author__ = "Jean-Philippe ARGAUD"
  30
  31 import os, sys
  32 import numpy
  33 import Logging ; Logging.Logging() # A importer en premier
  34 import Persistence
  35 from BasicObjects import Operator
  36 from PlatformInfo import uniq
  37
  38 # ==============================================================================
  39 class AssimilationStudy:
  40     """
  41     Cette classe sert d'interface pour l'utilisation de l'assimilation de
  42     données. Elle contient les méthodes ou accesseurs nécessaires à la
  43     construction d'un calcul d'assimilation.
  44     """
  45     def __init__(self, name=""):
  46         """
  47         Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
  48         élémentaire. Ces variables sont ensuite disponibles pour implémenter un
  49         algorithme élémentaire particulier.
  50
  51         Background............: vecteur Xb
  52         Observation...........: vecteur Y (potentiellement temporel)
  53             d'observations
  54         State.................: vecteur d'état dont une partie est le vecteur de
  55             contrôle. Cette information n'est utile que si l'on veut faire des
  56             calculs sur l'état complet, mais elle n'est pas indispensable pour
  57             l'assimilation.
  58         Control...............: vecteur X contenant toutes les variables de
  59             contrôle, i.e. les paramètres ou l'état dont on veut estimer la
  60             valeur pour obtenir les observations
  61         ObservationOperator...: opérateur d'observation H
  62
  63         Les observations présentent une erreur dont la matrice de covariance est
  64         R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
  65         de covariance est B.
  66         """
  67         self.__name = str(name)
  68         self.__Xb = None
  69         self.__Y  = None
  70         self.__B  = None
  71         self.__R  = None
  72         self.__Q  = None
  73         self.__H  = {}
  74         self.__M  = {}
  75         #
  76         self.__X  = Persistence.OneVector()
  77         self.__Parameters        = {}
  78         self.__StoredDiagnostics = {}
  79         self.__StoredInputs      = {}
  80         #
  81         self.__B_scalar = None
  82         self.__R_scalar = None
  83         self.__Q_scalar = None
  84         self.__Parameters["B_scalar"] = None
  85         self.__Parameters["R_scalar"] = None
  86         self.__Parameters["Q_scalar"] = None
  87         #
  88         # Variables temporaires
  89         self.__algorithm         = {}
  90         self.__algorithmFile     = None
  91         self.__algorithmName     = None
  92         self.__diagnosticFile    = None
  93         #
  94         # Récupère le chemin du répertoire parent et l'ajoute au path
  95         # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
  96         # qui est activée dans Persistence)
  97         self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
  98         sys.path.insert(0, self.__parent)
  99         sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
 100
 101     # ---------------------------------------------------------
 102     def setBackground(self,
 103             asVector           = None,
 104             asPersistentVector = None,
 105             Scheduler          = None,
 106             toBeStored         = False,
 107             ):
 108         """
 109         Permet de définir l'estimation a priori :
 110         - asVector : entrée des données, comme un vecteur compatible avec le
 111           constructeur de numpy.matrix
 112         - asPersistentVector : entrée des données, comme un vecteur de type
 113           persistent contruit avec la classe ad-hoc "Persistence"
 114         - Scheduler est le contrôle temporel des données
 115         - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
 116           être rendue disponible au même titre que les variables de calcul
 117         """
 118         if asVector is not None:
 119             if type( asVector ) is type( numpy.matrix([]) ):
 120                 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
 121             else:
 122                 self.__Xb = numpy.matrix( asVector,    numpy.float ).T
 123         elif asPersistentVector is not None:
 124             if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
 125                 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
 126                 for member in asPersistentVector:
 127                     self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
 128             else:
 129                 self.__Xb = asPersistentVector
 130         else:
 131             raise ValueError("Error: improperly defined background")
 132         if toBeStored:
 133            self.__StoredInputs["Background"] = self.__Xb
 134         return 0
 135
 136     def setBackgroundError(self,
 137             asCovariance  = None,
 138             asEyeByScalar = None,
 139             asEyeByVector = None,
 140             toBeStored    = False,
 141             ):
 142         """
 143         Permet de définir la covariance des erreurs d'ébauche :
 144         - asCovariance : entrée des données, comme une matrice compatible avec
 145           le constructeur de numpy.matrix
 146         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 147           multiplicatif d'une matrice de corrélation identité, aucune matrice
 148           n'étant donc explicitement à donner
 149         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 150           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 151           n'étant donc explicitement à donner
 152         - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
 153           être rendue disponible au même titre que les variables de calcul
 154         """
 155         if asEyeByScalar is not None:
 156             self.__B_scalar = float(asEyeByScalar)
 157             self.__B        = None
 158         elif asEyeByVector is not None:
 159             self.__B_scalar = None
 160             self.__B        = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
 161         elif asCovariance is not None:
 162             self.__B_scalar = None
 163             self.__B        = numpy.matrix( asCovariance, numpy.float )
 164         else:
 165             raise ValueError("Background error covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix")
 166         #
 167         self.__Parameters["B_scalar"] = self.__B_scalar
 168         if toBeStored:
 169             self.__StoredInputs["BackgroundError"] = self.__B
 170         return 0
 171
 172     # -----------------------------------------------------------
 173     def setObservation(self,
 174             asVector           = None,
 175             asPersistentVector = None,
 176             Scheduler          = None,
 177             toBeStored         = False,
 178             ):
 179         """
 180         Permet de définir les observations :
 181         - asVector : entrée des données, comme un vecteur compatible avec le
 182           constructeur de numpy.matrix
 183         - asPersistentVector : entrée des données, comme un vecteur de type
 184           persistent contruit avec la classe ad-hoc "Persistence"
 185         - Scheduler est le contrôle temporel des données disponibles
 186         - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
 187           être rendue disponible au même titre que les variables de calcul
 188         """
 189         if asVector is not None:
 190             if type( asVector ) is type( numpy.matrix([]) ):
 191                 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
 192             else:
 193                 self.__Y = numpy.matrix( asVector,    numpy.float ).T
 194         elif asPersistentVector is not None:
 195             if type( asPersistentVector ) is list or type( asPersistentVector ) is tuple:
 196                 from Persistence import OneVector
 197                 self.__Y = OneVector("Observation", basetype=numpy.array)
 198                 for y in asPersistentVector:
 199                     self.__Y.store( y )
 200             else:
 201                 self.__Y = asPersistentVector
 202         else:
 203             raise ValueError("Error: improperly defined observations")
 204         if toBeStored:
 205             self.__StoredInputs["Observation"] = self.__Y
 206         return 0
 207
 208     def setObservationError(self,
 209             asCovariance  = None,
 210             asEyeByScalar = None,
 211             asEyeByVector = None,
 212             toBeStored    = False,
 213             ):
 214         """
 215         Permet de définir la covariance des erreurs d'observations :
 216         - asCovariance : entrée des données, comme une matrice compatible avec
 217           le constructeur de numpy.matrix
 218         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 219           multiplicatif d'une matrice de corrélation identité, aucune matrice
 220           n'étant donc explicitement à donner
 221         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 222           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 223           n'étant donc explicitement à donner
 224         - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
 225           être rendue disponible au même titre que les variables de calcul
 226         """
 227         if asEyeByScalar is not None:
 228             self.__R_scalar = float(asEyeByScalar)
 229             self.__R        = None
 230         elif asEyeByVector is not None:
 231             self.__R_scalar = None
 232             self.__R        = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
 233         elif asCovariance is not None:
 234             self.__R_scalar = None
 235             self.__R        = numpy.matrix( asCovariance, numpy.float )
 236         else:
 237             raise ValueError("Observation error covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix")
 238         #
 239         self.__Parameters["R_scalar"] = self.__R_scalar
 240         if toBeStored:
 241             self.__StoredInputs["ObservationError"] = self.__R
 242         return 0
 243
 244     def setObservationOperator(self,
 245             asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
 246                           "useApproximatedDerivatives":False,
 247                           "withCenteredDF"            :False,
 248                           "withIncrement"             :0.01,
 249                           "withdX"                    :None,
 250                          },
 251             asMatrix   = None,
 252             appliedToX = None,
 253             toBeStored = False,
 254             ):
 255         """
 256         Permet de définir un opérateur d'observation H. L'ordre de priorité des
 257         définitions et leur sens sont les suivants :
 258         - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
 259           alors on définit l'opérateur à l'aide de fonctions. Si la fonction
 260           "Direct" n'est pas définie, on prend la fonction "Tangent".
 261           Si "useApproximatedDerivatives" est vrai, on utilise une approximation
 262           des opérateurs tangents et adjoints. On utilise par défaut des
 263           différences finies non centrées ou centrées (si "withCenteredDF" est
 264           vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
 265           du point courant ou sur le point fixe "withdX".
 266         - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
 267           None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
 268           la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
 269           matrice fournie doit être sous une forme compatible avec le
 270           constructeur de numpy.matrix.
 271         - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
 272           X divers, l'opérateur par sa valeur appliquée à cet X particulier,
 273           sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
 274           L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
 275         - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
 276           être rendue disponible au même titre que les variables de calcul
 277         """
 278         if (type(asFunction) is type({})) and \
 279                 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
 280                 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
 281             if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
 282             if not asFunction.has_key("withIncrement"):  asFunction["withIncrement"]  = 0.01
 283             if not asFunction.has_key("withdX"):         asFunction["withdX"]         = None
 284             from daNumerics.ApproximatedDerivatives import FDApproximation
 285             FDA = FDApproximation(
 286                 FunctionH  = asFunction["Direct"],
 287                 centeredDF = asFunction["withCenteredDF"],
 288                 increment  = asFunction["withIncrement"],
 289                 dX         = asFunction["withdX"] )
 290             self.__H["Direct"]  = Operator( fromMethod = FDA.FunctionH )
 291             self.__H["Tangent"] = Operator( fromMethod = FDA.TangentH  )
 292             self.__H["Adjoint"] = Operator( fromMethod = FDA.AdjointH  )
 293         elif (type(asFunction) is type({})) and \
 294                 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
 295                 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
 296             if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
 297                 self.__H["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
 298             else:
 299                 self.__H["Direct"] = Operator( fromMethod = asFunction["Direct"]  )
 300             self.__H["Tangent"]    = Operator( fromMethod = asFunction["Tangent"] )
 301             self.__H["Adjoint"]    = Operator( fromMethod = asFunction["Adjoint"] )
 302         elif asMatrix is not None:
 303             matrice = numpy.matrix( asMatrix, numpy.float )
 304             self.__H["Direct"]  = Operator( fromMatrix = matrice )
 305             self.__H["Tangent"] = Operator( fromMatrix = matrice )
 306             self.__H["Adjoint"] = Operator( fromMatrix = matrice.T )
 307         else:
 308             raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
 309         #
 310         if appliedToX is not None:
 311             self.__H["AppliedToX"] = {}
 312             if type(appliedToX) is not dict:
 313                 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
 314             for key in appliedToX.keys():
 315                 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
 316                     # Pour le cas où l'on a une vraie matrice
 317                     self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
 318                 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
 319                     # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
 320                     self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
 321                 else:
 322                     self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key],    numpy.float ).T
 323         else:
 324             self.__H["AppliedToX"] = None
 325         #
 326         if toBeStored:
 327             self.__StoredInputs["ObservationOperator"] = self.__H
 328         return 0
 329
 330     # -----------------------------------------------------------
 331     def setEvolutionModel(self,
 332             asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
 333                           "useApproximatedDerivatives":False,
 334                           "withCenteredDF"            :False,
 335                           "withIncrement"             :0.01,
 336                           "withdX"                    :None,
 337                          },
 338             asMatrix   = None,
 339             Scheduler  = None,
 340             toBeStored = False,
 341             ):
 342         """
 343         Permet de définir un opérateur d'évolution M. L'ordre de priorité des
 344         définitions et leur sens sont les suivants :
 345         - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
 346           alors on définit l'opérateur à l'aide de fonctions. Si la fonction
 347           "Direct" n'est pas définie, on prend la fonction "Tangent".
 348           Si "useApproximatedDerivatives" est vrai, on utilise une approximation
 349           des opérateurs tangents et adjoints. On utilise par défaut des
 350           différences finies non centrées ou centrées (si "withCenteredDF" est
 351           vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
 352           du point courant ou sur le point fixe "withdX".
 353         - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
 354           None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
 355           la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
 356           matrice fournie doit être sous une forme compatible avec le
 357           constructeur de numpy.matrix.
 358         - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
 359           être rendue disponible au même titre que les variables de calcul
 360         """
 361         if (type(asFunction) is type({})) and \
 362                 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
 363                 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
 364             if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
 365             if not asFunction.has_key("withIncrement"):  asFunction["withIncrement"]  = 0.01
 366             if not asFunction.has_key("withdX"):         asFunction["withdX"]         = None
 367             from daNumerics.ApproximatedDerivatives import FDApproximation
 368             FDA = FDApproximation(
 369                 FunctionH  = asFunction["Direct"],
 370                 centeredDF = asFunction["withCenteredDF"],
 371                 increment  = asFunction["withIncrement"],
 372                 dX         = asFunction["withdX"] )
 373             self.__M["Direct"]  = Operator( fromMethod = FDA.FunctionH )
 374             self.__M["Tangent"] = Operator( fromMethod = FDA.TangentH  )
 375             self.__M["Adjoint"] = Operator( fromMethod = FDA.AdjointH  )
 376         elif (type(asFunction) is type({})) and \
 377                 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
 378                 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
 379             if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
 380                 self.__M["Direct"] = Operator( fromMethod = asFunction["Tangent"]  )
 381             else:
 382                 self.__M["Direct"] = Operator( fromMethod = asFunction["Direct"]  )
 383             self.__M["Tangent"]    = Operator( fromMethod = asFunction["Tangent"] )
 384             self.__M["Adjoint"]    = Operator( fromMethod = asFunction["Adjoint"] )
 385         elif asMatrix is not None:
 386             matrice = numpy.matrix( asMatrix, numpy.float )
 387             self.__M["Direct"]  = Operator( fromMatrix = matrice )
 388             self.__M["Tangent"] = Operator( fromMatrix = matrice )
 389             self.__M["Adjoint"] = Operator( fromMatrix = matrice.T )
 390         else:
 391             raise ValueError("Improperly defined evolution operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
 392         #
 393         if toBeStored:
 394             self.__StoredInputs["EvolutionModel"] = self.__M
 395         return 0
 396
 397     def setEvolutionError(self,
 398             asCovariance  = None,
 399             asEyeByScalar = None,
 400             asEyeByVector = None,
 401             toBeStored    = False,
 402             ):
 403         """
 404         Permet de définir la covariance des erreurs de modèle :
 405         - asCovariance : entrée des données, comme une matrice compatible avec
 406           le constructeur de numpy.matrix
 407         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
 408           multiplicatif d'une matrice de corrélation identité, aucune matrice
 409           n'étant donc explicitement à donner
 410         - asEyeByVector : entrée des données comme un seul vecteur de variance,
 411           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
 412           n'étant donc explicitement à donner
 413         - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
 414           être rendue disponible au même titre que les variables de calcul
 415         """
 416         if asEyeByScalar is not None:
 417             self.__Q_scalar = float(asEyeByScalar)
 418             self.__Q        = None
 419         elif asEyeByVector is not None:
 420             self.__Q_scalar = None
 421             self.__Q        = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
 422         elif asCovariance is not None:
 423             self.__Q_scalar = None
 424             self.__Q        = numpy.matrix( asCovariance, numpy.float )
 425         else:
 426             raise ValueError("Evolution error covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix")
 427         #
 428         self.__Parameters["Q_scalar"] = self.__Q_scalar
 429         if toBeStored:
 430             self.__StoredInputs["EvolutionError"] = self.__Q
 431         return 0
 432
 433     # -----------------------------------------------------------
 434     def setControls (self,
 435             asVector = None,
 436             toBeStored   = False,
 437             ):
 438         """
 439         Permet de définir la valeur initiale du vecteur X contenant toutes les
 440         variables de contrôle, i.e. les paramètres ou l'état dont on veut
 441         estimer la valeur pour obtenir les observations. C'est utile pour un
 442         algorithme itératif/incrémental.
 443         - asVector : entrée des données, comme un vecteur compatible avec le
 444           constructeur de numpy.matrix.
 445         - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
 446           être rendue disponible au même titre que les variables de calcul
 447         """
 448         if asVector is not None:
 449             self.__X.store( asVector )
 450         if toBeStored:
 451             self.__StoredInputs["Controls"] = self.__X
 452         return 0
 453
 454     # -----------------------------------------------------------
 455     def setAlgorithm(self, choice = None ):
 456         """
 457         Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
 458         d'assimilation. L'argument est un champ caractère se rapportant au nom
 459         d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
 460         d'assimilation sur les arguments fixes.
 461         """
 462         if choice is None:
 463             raise ValueError("Error: algorithm choice has to be given")
 464         if self.__algorithmName is not None:
 465             raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
 466         daDirectory = "daAlgorithms"
 467         #
 468         # Recherche explicitement le fichier complet
 469         # ------------------------------------------
 470         module_path = None
 471         for directory in sys.path:
 472             if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
 473                 module_path = os.path.abspath(os.path.join(directory, daDirectory))
 474         if module_path is None:
 475             raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n             The search path is %s"%(choice, daDirectory, sys.path))
 476         #
 477         # Importe le fichier complet comme un module
 478         # ------------------------------------------
 479         try:
 480             sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
 481             self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
 482             self.__algorithmName = str(choice)
 483             sys.path = sys_path_tmp ; del sys_path_tmp
 484         except ImportError, e:
 485             raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n             The import error message is: %s"%(choice,e))
 486         #
 487         # Instancie un objet du type élémentaire du fichier
 488         # -------------------------------------------------
 489         self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
 490         self.__StoredInputs["AlgorithmName"] = self.__algorithmName
 491         return 0
 492
 493     def setAlgorithmParameters(self, asDico=None):
 494         """
 495         Permet de définir les paramètres de l'algorithme, sous la forme d'un
 496         dictionnaire.
 497         """
 498         if asDico is not None:
 499             self.__Parameters.update( dict( asDico ) )
 500         #
 501         self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
 502         return 0
 503
 504     def getAlgorithmParameters(self, noDetails=True):
 505         """
 506         Renvoie la liste des paramètres requis selon l'algorithme
 507         """
 508         return self.__algorithm.getRequiredParameters(noDetails)
 509
 510     # -----------------------------------------------------------
 511     def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
 512         """
 513         Permet de sélectionner un diagnostic a effectuer.
 514         """
 515         if choice is None:
 516             raise ValueError("Error: diagnostic choice has to be given")
 517         daDirectory = "daDiagnostics"
 518         #
 519         # Recherche explicitement le fichier complet
 520         # ------------------------------------------
 521         module_path = None
 522         for directory in sys.path:
 523             if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
 524                 module_path = os.path.abspath(os.path.join(directory, daDirectory))
 525         if module_path is None:
 526             raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n             The search path is %s"%(choice, daDirectory, sys.path))
 527         #
 528         # Importe le fichier complet comme un module
 529         # ------------------------------------------
 530         try:
 531             sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
 532             self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
 533             sys.path = sys_path_tmp ; del sys_path_tmp
 534         except ImportError, e:
 535             raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n             The import error message is: %s"%(choice,e))
 536         #
 537         # Instancie un objet du type élémentaire du fichier
 538         # -------------------------------------------------
 539         if self.__StoredInputs.has_key(name):
 540             raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
 541         elif self.__StoredDiagnostics.has_key(name):
 542             raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
 543         else:
 544             self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
 545                 name       = name,
 546                 unit       = unit,
 547                 basetype   = basetype,
 548                 parameters = parameters )
 549         return 0
 550
 551     # -----------------------------------------------------------
 552     def shape_validate(self):
 553         """
 554         Validation de la correspondance correcte des tailles des variables et
 555         des matrices s'il y en a.
 556         """
 557         if self.__Xb is None:                  __Xb_shape = (0,)
 558         elif hasattr(self.__Xb,"shape"):
 559             if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
 560             else:                              __Xb_shape = self.__Xb.shape()
 561         else: raise TypeError("Xb has no attribute of shape: problem !")
 562         #
 563         if self.__Y is None:                  __Y_shape = (0,)
 564         elif hasattr(self.__Y,"shape"):
 565             if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
 566             else:                             __Y_shape = self.__Y.shape()
 567         else: raise TypeError("Y has no attribute of shape: problem !")
 568         #
 569         if self.__B is None and self.__B_scalar is None:
 570             __B_shape = (0,0)
 571         elif self.__B is None and self.__B_scalar is not None:
 572             __B_shape = (max(__Xb_shape),max(__Xb_shape))
 573         elif hasattr(self.__B,"shape"):
 574             if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
 575             else:                             __B_shape = self.__B.shape()
 576         else: raise TypeError("B has no attribute of shape: problem !")
 577         #
 578         if self.__R is None and self.__R_scalar is None:
 579             __R_shape = (0,0)
 580         elif self.__R is None and self.__R_scalar is not None:
 581             __R_shape = (max(__Y_shape),max(__Y_shape))
 582         elif hasattr(self.__R,"shape"):
 583             if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
 584             else:                             __R_shape = self.__R.shape()
 585         else: raise TypeError("R has no attribute of shape: problem !")
 586         #
 587         if self.__Q is None:                  __Q_shape = (0,0)
 588         elif hasattr(self.__Q,"shape"):
 589             if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
 590             else:                             __Q_shape = self.__Q.shape()
 591         else: raise TypeError("Q has no attribute of shape: problem !")
 592         #
 593         if len(self.__H) == 0:                          __H_shape = (0,0)
 594         elif type(self.__H) is type({}):                __H_shape = (0,0)
 595         elif hasattr(self.__H["Direct"],"shape"):
 596             if type(self.__H["Direct"].shape) is tuple: __H_shape = self.__H["Direct"].shape
 597             else:                                       __H_shape = self.__H["Direct"].shape()
 598         else: raise TypeError("H has no attribute of shape: problem !")
 599         #
 600         if len(self.__M) == 0:                          __M_shape = (0,0)
 601         elif type(self.__M) is type({}):                __M_shape = (0,0)
 602         elif hasattr(self.__M["Direct"],"shape"):
 603             if type(self.__M["Direct"].shape) is tuple: __M_shape = self.__M["Direct"].shape
 604             else:                                       __M_shape = self.__M["Direct"].shape()
 605         else: raise TypeError("M has no attribute of shape: problem !")
 606         #
 607         # Vérification des conditions
 608         # ---------------------------
 609         if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
 610             raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
 611         if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
 612             raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
 613         #
 614         if not( min(__B_shape) == max(__B_shape) ):
 615             raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
 616         if not( min(__R_shape) == max(__R_shape) ):
 617             raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
 618         if not( min(__Q_shape) == max(__Q_shape) ):
 619             raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
 620         if not( min(__M_shape) == max(__M_shape) ):
 621             raise ValueError("Shape characteristic of M is incorrect: \"%s\""%(__M_shape,))
 622         #
 623         if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[1] == max(__Xb_shape) ):
 624             raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__H_shape,__Xb_shape))
 625         if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[0] == max(__Y_shape) ):
 626             raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__H_shape,__Y_shape))
 627         if len(self.__H) > 0 and not(type(self.__H) is type({})) and len(self.__B) > 0 and not( __H_shape[1] == __B_shape[0] ):
 628             raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__H_shape,__B_shape))
 629         if len(self.__H) > 0 and not(type(self.__H) is type({})) and len(self.__R) > 0 and not( __H_shape[0] == __R_shape[1] ):
 630             raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__H_shape,__R_shape))
 631         #
 632         if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
 633             if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
 634                 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
 635                 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
 636                 for member in asPersistentVector:
 637                     self.__Xb.store( numpy.matrix( numpy.asarray(member).flatten(), numpy.float ).T )
 638                 __Xb_shape = min(__B_shape)
 639             else:
 640                 raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
 641         #
 642         if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
 643             raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
 644         #
 645         if self.__M is not None and len(self.__M) > 0 and not(type(self.__M) is type({})) and not( __M_shape[1] == max(__Xb_shape) ):
 646             raise ValueError("Shape characteristic of M \"%s\" and X \"%s\" are incompatible"%(__M_shape,__Xb_shape))
 647         #
 648         return 1
 649
 650     # -----------------------------------------------------------
 651     def analyze(self):
 652         """
 653         Permet de lancer le calcul d'assimilation.
 654
 655         Le nom de la méthode à activer est toujours "run". Les paramètres en
 656         arguments de la méthode sont fixés. En sortie, on obtient les résultats
 657         dans la variable de type dictionnaire "StoredVariables", qui contient en
 658         particulier des objets de Persistence pour les analyses, OMA...
 659         """
 660         self.shape_validate()
 661         #
 662         self.__algorithm.run(
 663             Xb         = self.__Xb,
 664             Y          = self.__Y,
 665             H          = self.__H,
 666             M          = self.__M,
 667             R          = self.__R,
 668             B          = self.__B,
 669             Q          = self.__Q,
 670             Parameters = self.__Parameters,
 671             )
 672         return 0
 673
 674     # -----------------------------------------------------------
 675     def get(self, key=None):
 676         """
 677         Renvoie les résultats disponibles après l'exécution de la méthode
 678         d'assimilation, ou les diagnostics disponibles. Attention, quand un
 679         diagnostic porte le même nom qu'une variable stockée, c'est la variable
 680         stockée qui est renvoyée, et le diagnostic est inatteignable.
 681         """
 682         if key is not None:
 683             if self.__algorithm.has_key(key):
 684                 return self.__algorithm.get( key )
 685             elif self.__StoredInputs.has_key(key):
 686                 return self.__StoredInputs[key]
 687             elif self.__StoredDiagnostics.has_key(key):
 688                 return self.__StoredDiagnostics[key]
 689             else:
 690                 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
 691         else:
 692             allvariables = self.__algorithm.get()
 693             allvariables.update( self.__StoredDiagnostics )
 694             allvariables.update( self.__StoredInputs )
 695             return allvariables
 696
 697     def get_available_variables(self):
 698         """
 699         Renvoie les variables potentiellement utilisables pour l'étude,
 700         initialement stockées comme données d'entrées ou dans les algorithmes,
 701         identifiés par les chaînes de caractères. L'algorithme doit avoir été
 702         préalablement choisi sinon la méthode renvoie "None".
 703         """
 704         if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
 705             return None
 706         else:
 707             variables = []
 708             if len( self.__algorithm.keys()) > 0:
 709                 variables.extend( self.__algorithm.get().keys() )
 710             if len( self.__StoredInputs.keys() ) > 0:
 711                 variables.extend( self.__StoredInputs.keys() )
 712             variables.sort()
 713             return variables
 714
 715     def get_available_algorithms(self):
 716         """
 717         Renvoie la liste des algorithmes potentiellement utilisables, identifiés
 718         par les chaînes de caractères.
 719         """
 720         files = []
 721         for directory in sys.path:
 722             if os.path.isdir(os.path.join(directory,"daAlgorithms")):
 723                 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
 724                     root, ext = os.path.splitext(fname)
 725                     if ext == '.py' and root != '__init__':
 726                         files.append(root)
 727         files.sort()
 728         return files
 729
 730     def get_available_diagnostics(self):
 731         """
 732         Renvoie la liste des diagnostics potentiellement utilisables, identifiés
 733         par les chaînes de caractères.
 734         """
 735         files = []
 736         for directory in sys.path:
 737             if os.path.isdir(os.path.join(directory,"daDiagnostics")):
 738                 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
 739                     root, ext = os.path.splitext(fname)
 740                     if ext == '.py' and root != '__init__':
 741                         files.append(root)
 742         files.sort()
 743         return files
 744
 745     # -----------------------------------------------------------
 746     def get_algorithms_main_path(self):
 747         """
 748         Renvoie le chemin pour le répertoire principal contenant les algorithmes
 749         dans un sous-répertoire "daAlgorithms"
 750         """
 751         return self.__parent
 752
 753     def add_algorithms_path(self, asPath=None):
 754         """
 755         Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
 756         se trouve un sous-répertoire "daAlgorithms"
 757
 758         Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
 759         pas indispensable de le rajouter ici.
 760         """
 761         if not os.path.isdir(asPath):
 762             raise ValueError("The given "+asPath+" argument must exist as a directory")
 763         if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
 764             raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
 765         if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
 766             raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
 767         sys.path.insert(0, os.path.abspath(asPath))
 768         sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
 769         return 1
 770
 771     def get_diagnostics_main_path(self):
 772         """
 773         Renvoie le chemin pour le répertoire principal contenant les diagnostics
 774         dans un sous-répertoire "daDiagnostics"
 775         """
 776         return self.__parent
 777
 778     def add_diagnostics_path(self, asPath=None):
 779         """
 780         Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
 781         se trouve un sous-répertoire "daDiagnostics"
 782
 783         Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
 784         pas indispensable de le rajouter ici.
 785         """
 786         if not os.path.isdir(asPath):
 787             raise ValueError("The given "+asPath+" argument must exist as a directory")
 788         if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
 789             raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
 790         if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
 791             raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
 792         sys.path.insert(0, os.path.abspath(asPath))
 793         sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
 794         return 1
 795
 796     # -----------------------------------------------------------
 797     def setDataObserver(self,
 798             VariableName   = None,
 799             HookFunction   = None,
 800             HookParameters = None,
 801             Scheduler      = None,
 802             ):
 803         """
 804         Permet d'associer un observer à une ou des variables nommées gérées en
 805         interne, activable selon des règles définies dans le Scheduler. A chaque
 806         pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
 807         les arguments (variable persistante VariableName, paramètres HookParameters).
 808         """
 809         #
 810         if type( self.__algorithm ) is dict:
 811             raise ValueError("No observer can be build before choosing an algorithm.")
 812         #
 813         # Vérification du nom de variable et typage
 814         # -----------------------------------------
 815         if type( VariableName ) is str:
 816             VariableNames = [VariableName,]
 817         elif type( VariableName ) is list:
 818             VariableNames = map( str, VariableName )
 819         else:
 820             raise ValueError("The observer requires a name or a list of names of variables.")
 821         #
 822         # Association interne de l'observer à la variable
 823         # -----------------------------------------------
 824         for n in VariableNames:
 825             if not self.__algorithm.has_key( n ):
 826                 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
 827         else:
 828             self.__algorithm.StoredVariables[ n ].setDataObserver(
 829                 Scheduler      = Scheduler,
 830                 HookFunction   = HookFunction,
 831                 HookParameters = HookParameters,
 832                 )
 833
 834     def removeDataObserver(self,
 835             VariableName   = None,
 836             HookFunction   = None,
 837             ):
 838         """
 839         Permet de retirer un observer à une ou des variable nommée.
 840         """
 841         #
 842         if type( self.__algorithm ) is dict:
 843             raise ValueError("No observer can be removed before choosing an algorithm.")
 844         #
 845         # Vérification du nom de variable et typage
 846         # -----------------------------------------
 847         if type( VariableName ) is str:
 848             VariableNames = [VariableName,]
 849         elif type( VariableName ) is list:
 850             VariableNames = map( str, VariableName )
 851         else:
 852             raise ValueError("The observer requires a name or a list of names of variables.")
 853         #
 854         # Association interne de l'observer à la variable
 855         # -----------------------------------------------
 856         for n in VariableNames:
 857             if not self.__algorithm.has_key( n ):
 858                 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
 859         else:
 860             self.__algorithm.StoredVariables[ n ].removeDataObserver(
 861                 HookFunction   = HookFunction,
 862                 )
 863
 864     # -----------------------------------------------------------
 865     def setDebug(self, level=10):
 866         """
 867         Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
 868         appel pour changer le niveau de verbosité, avec :
 869         NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
 870         """
 871         import logging
 872         log = logging.getLogger()
 873         log.setLevel( level )
 874
 875     def unsetDebug(self):
 876         """
 877         Remet le logger au niveau par défaut
 878         """
 879         import logging
 880         log = logging.getLogger()
 881         log.setLevel( logging.WARNING )
 882
 883     def prepare_to_pickle(self):
 884         self.__algorithmFile = None
 885         self.__diagnosticFile = None
 886         self.__H  = {}
 887         self.__M  = {}
 888
 889 # ==============================================================================
 890 if __name__ == "__main__":
 891     print '\n AUTODIAGNOSTIC \n'
 892
 893     ADD = AssimilationStudy("Ma premiere etude BLUE")
 894
 895     ADD.setBackground         (asVector     = [0, 1, 2])
 896     ADD.setBackgroundError    (asCovariance = "1 0 0;0 1 0;0 0 1")
 897     ADD.setObservation        (asVector     = [0.5, 1.5, 2.5])
 898     ADD.setObservationError   (asCovariance = "1 0 0;0 1 0;0 0 1")
 899     ADD.setObservationOperator(asMatrix     = "1 0 0;0 1 0;0 0 1")
 900
 901     ADD.setAlgorithm(choice="Blue")
 902
 903     ADD.analyze()
 904
 905     print "Nombre d'analyses  :", ADD.get("Analysis").stepnumber()
 906     print "Ebauche            :", [0, 1, 2]
 907     print "Observation        :", [0.5, 1.5, 2.5]
 908     print "Demi-somme         :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
 909     print "  qui doit être identique à :"
 910     print "Analyse résultante :", ADD.get("Analysis").valueserie(0)
 911     print
 912
 913     print "Algorithmes disponibles.......................:", ADD.get_available_algorithms()
 914     # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path()
 915     print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics()
 916     # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path()
 917     print "Variables disponibles.........................:", ADD.get_available_variables()
 918     print
 919
 920     print "Paramètres requis par algorithme :"
 921     for algo in ADD.get_available_algorithms():
 922         tmpADD = AssimilationStudy("Un algorithme")
 923         tmpADD.setAlgorithm(choice=algo)
 924         print "  %25s : %s"%(algo,tmpADD.getAlgorithmParameters())
 925         del tmpADD
 926     print
 927
 928     ADD.setDiagnostic("RMS", "Ma RMS")
 929
 930     liste = ADD.get().keys()
 931     liste.sort()
 932     print "Variables et diagnostics nommés disponibles...:", liste
 933
 934     print
 935     print "Exemple de mise en place d'un observeur :"
 936     def obs(var=None,info=None):
 937         print "  ---> Mise en oeuvre de l'observer"
 938         print "       var  =",var.valueserie(-1)
 939         print "       info =",info
 940     ADD.setDataObserver( 'Analysis', HookFunction=obs, Scheduler = [2, 4], HookParameters = "Second observer")
 941     # Attention, il faut décaler le stockage de 1 pour suivre le pas interne
 942     # car le pas 0 correspond à l'analyse ci-dessus.
 943     for i in range(1,6):
 944         print
 945         print "Action sur la variable observée, étape :",i
 946         ADD.get('Analysis').store( [i, i, i] )
 947     print
 948
 949     print "Mise en debug et hors debug"
 950     print "Nombre d'analyses  :", ADD.get("Analysis").stepnumber()
 951     ADD.setDebug()
 952     ADD.analyze()
 953     ADD.unsetDebug()
 954     print "Nombre d'analyses  :", ADD.get("Analysis").stepnumber()
 955     ADD.analyze()
 956     print "Nombre d'analyses  :", ADD.get("Analysis").stepnumber()
 957     print