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