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