src/daComposant/daCore/AssimilationStudy.py

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