src/daComposant/daCore/AssimilationStudy.py

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