src/daComposant/daCore/AssimilationStudy.py

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