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