1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2013 EDF R&D
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.
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.
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
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 Classe principale pour la préparation, la réalisation et la restitution de
25 calculs d'assimilation de données.
27 Ce module est destiné à être appelé par AssimilationStudy pour constituer
28 les objets élémentaires de l'étude.
30 __author__ = "Jean-Philippe ARGAUD"
34 import Logging ; Logging.Logging() # A importer en premier
37 from BasicObjects import Operator
38 from PlatformInfo import uniq
40 # ==============================================================================
41 class AssimilationStudy:
43 Cette classe sert d'interface pour l'utilisation de l'assimilation de
44 données. Elle contient les méthodes ou accesseurs nécessaires à la
45 construction d'un calcul d'assimilation.
47 def __init__(self, name=""):
49 Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
50 élémentaire. Ces variables sont ensuite disponibles pour implémenter un
51 algorithme élémentaire particulier.
53 Background............: vecteur Xb
54 Observation...........: vecteur Y (potentiellement temporel)
56 State.................: vecteur d'état dont une partie est le vecteur de
57 contrôle. Cette information n'est utile que si l'on veut faire des
58 calculs sur l'état complet, mais elle n'est pas indispensable pour
60 Control...............: vecteur X contenant toutes les variables de
61 contrôle, i.e. les paramètres ou l'état dont on veut estimer la
62 valeur pour obtenir les observations
63 ObservationOperator...: opérateur d'observation H
65 Les observations présentent une erreur dont la matrice de covariance est
66 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
69 self.__name = str(name)
80 self.__X = Persistence.OneVector()
81 self.__Parameters = {}
82 self.__StoredDiagnostics = {}
83 self.__StoredInputs = {}
85 self.__B_scalar = None
86 self.__R_scalar = None
87 self.__Q_scalar = None
88 self.__Parameters["B_scalar"] = None
89 self.__Parameters["R_scalar"] = None
90 self.__Parameters["Q_scalar"] = None
92 # Variables temporaires
94 self.__algorithmFile = None
95 self.__algorithmName = None
96 self.__diagnosticFile = None
98 # Récupère le chemin du répertoire parent et l'ajoute au path
99 # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
100 # qui est activée dans Persistence)
101 self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
102 sys.path.insert(0, self.__parent)
103 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
105 # ---------------------------------------------------------
106 def setBackground(self,
108 asPersistentVector = None,
113 Permet de définir l'estimation a priori :
114 - asVector : entrée des données, comme un vecteur compatible avec le
115 constructeur de numpy.matrix
116 - asPersistentVector : entrée des données, comme un vecteur de type
117 persistent contruit avec la classe ad-hoc "Persistence"
118 - Scheduler est le contrôle temporel des données
119 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
120 être rendue disponible au même titre que les variables de calcul
122 if asVector is not None:
123 if type( asVector ) is type( numpy.matrix([]) ):
124 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
126 self.__Xb = numpy.matrix( asVector, numpy.float ).T
127 elif asPersistentVector is not None:
128 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
129 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
130 for member in asPersistentVector:
131 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
133 self.__Xb = asPersistentVector
135 raise ValueError("Error: improperly defined background, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
137 self.__StoredInputs["Background"] = self.__Xb
140 def setBackgroundError(self,
142 asEyeByScalar = None,
143 asEyeByVector = None,
147 Permet de définir la covariance des erreurs d'ébauche :
148 - asCovariance : entrée des données, comme une matrice compatible avec
149 le constructeur de numpy.matrix
150 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
151 multiplicatif d'une matrice de corrélation identité, aucune matrice
152 n'étant donc explicitement à donner
153 - asEyeByVector : entrée des données comme un seul vecteur de variance,
154 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
155 n'étant donc explicitement à donner
156 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
157 être rendue disponible au même titre que les variables de calcul
159 if asEyeByScalar is not None:
160 self.__B_scalar = float(asEyeByScalar)
162 elif asEyeByVector is not None:
163 self.__B_scalar = None
164 self.__B = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
165 elif asCovariance is not None:
166 self.__B_scalar = None
167 self.__B = numpy.matrix( asCovariance, numpy.float )
169 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")
171 self.__Parameters["B_scalar"] = self.__B_scalar
173 self.__StoredInputs["BackgroundError"] = self.__B
176 # -----------------------------------------------------------
177 def setObservation(self,
179 asPersistentVector = None,
184 Permet de définir les observations :
185 - asVector : entrée des données, comme un vecteur compatible avec le
186 constructeur de numpy.matrix
187 - asPersistentVector : entrée des données, comme un vecteur de type
188 persistent contruit avec la classe ad-hoc "Persistence"
189 - Scheduler est le contrôle temporel des données disponibles
190 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
191 être rendue disponible au même titre que les variables de calcul
193 if asVector is not None:
194 if type( asVector ) is type( numpy.matrix([]) ):
195 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
197 self.__Y = numpy.matrix( asVector, numpy.float ).T
198 elif asPersistentVector is not None:
199 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
200 self.__Y = Persistence.OneVector("Observation", basetype=numpy.matrix)
201 for member in asPersistentVector:
202 self.__Y.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
204 self.__Y = asPersistentVector
206 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
208 self.__StoredInputs["Observation"] = self.__Y
211 def setObservationError(self,
213 asEyeByScalar = None,
214 asEyeByVector = None,
218 Permet de définir la covariance des erreurs d'observations :
219 - asCovariance : entrée des données, comme une matrice compatible avec
220 le constructeur de numpy.matrix
221 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
222 multiplicatif d'une matrice de corrélation identité, aucune matrice
223 n'étant donc explicitement à donner
224 - asEyeByVector : entrée des données comme un seul vecteur de variance,
225 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
226 n'étant donc explicitement à donner
227 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
228 être rendue disponible au même titre que les variables de calcul
230 if asEyeByScalar is not None:
231 self.__R_scalar = float(asEyeByScalar)
233 elif asEyeByVector is not None:
234 self.__R_scalar = None
235 self.__R = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
236 elif asCovariance is not None:
237 self.__R_scalar = None
238 self.__R = numpy.matrix( asCovariance, numpy.float )
240 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")
242 self.__Parameters["R_scalar"] = self.__R_scalar
244 self.__StoredInputs["ObservationError"] = self.__R
247 def setObservationOperator(self,
248 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
249 "useApproximatedDerivatives":False,
250 "withCenteredDF" :False,
251 "withIncrement" :0.01,
259 Permet de définir un opérateur d'observation H. L'ordre de priorité des
260 définitions et leur sens sont les suivants :
261 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
262 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
263 "Direct" n'est pas définie, on prend la fonction "Tangent".
264 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
265 des opérateurs tangents et adjoints. On utilise par défaut des
266 différences finies non centrées ou centrées (si "withCenteredDF" est
267 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
268 du point courant ou sur le point fixe "withdX".
269 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
270 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
271 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
272 matrice fournie doit être sous une forme compatible avec le
273 constructeur de numpy.matrix.
274 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
275 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
276 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
277 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
278 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
279 être rendue disponible au même titre que les variables de calcul
281 if (type(asFunction) is type({})) and \
282 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
283 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
284 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
285 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
286 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
287 from daNumerics.ApproximatedDerivatives import FDApproximation
288 FDA = FDApproximation(
289 Function = asFunction["Direct"],
290 centeredDF = asFunction["withCenteredDF"],
291 increment = asFunction["withIncrement"],
292 dX = asFunction["withdX"] )
293 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator )
294 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
295 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
296 elif (type(asFunction) is type({})) and \
297 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
298 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
299 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
300 self.__HO["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
302 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"] )
303 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
304 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
305 elif asMatrix is not None:
306 matrice = numpy.matrix( asMatrix, numpy.float )
307 self.__HO["Direct"] = Operator( fromMatrix = matrice )
308 self.__HO["Tangent"] = Operator( fromMatrix = matrice )
309 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T )
312 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
314 if appliedToX is not None:
315 self.__HO["AppliedToX"] = {}
316 if type(appliedToX) is not dict:
317 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
318 for key in appliedToX.keys():
319 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
320 # Pour le cas où l'on a une vraie matrice
321 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
322 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
323 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
324 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
326 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
328 self.__HO["AppliedToX"] = None
331 self.__StoredInputs["ObservationOperator"] = self.__HO
334 # -----------------------------------------------------------
335 def setEvolutionModel(self,
336 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
337 "useApproximatedDerivatives":False,
338 "withCenteredDF" :False,
339 "withIncrement" :0.01,
347 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
348 définitions et leur sens sont les suivants :
349 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
350 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
351 "Direct" n'est pas définie, on prend la fonction "Tangent".
352 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
353 des opérateurs tangents et adjoints. On utilise par défaut des
354 différences finies non centrées ou centrées (si "withCenteredDF" est
355 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
356 du point courant ou sur le point fixe "withdX".
357 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
358 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
359 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
360 matrice fournie doit être sous une forme compatible avec le
361 constructeur de numpy.matrix.
362 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
363 être rendue disponible au même titre que les variables de calcul
365 if (type(asFunction) is type({})) and \
366 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
367 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
368 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
369 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
370 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
371 from daNumerics.ApproximatedDerivatives import FDApproximation
372 FDA = FDApproximation(
373 Function = asFunction["Direct"],
374 centeredDF = asFunction["withCenteredDF"],
375 increment = asFunction["withIncrement"],
376 dX = asFunction["withdX"] )
377 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
378 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
379 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
380 elif (type(asFunction) is type({})) and \
381 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
382 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
383 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
384 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
386 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
387 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
388 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
389 elif asMatrix is not None:
390 matrice = numpy.matrix( asMatrix, numpy.float )
391 self.__EM["Direct"] = Operator( fromMatrix = matrice )
392 self.__EM["Tangent"] = Operator( fromMatrix = matrice )
393 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T )
396 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
399 self.__StoredInputs["EvolutionModel"] = self.__EM
402 def setEvolutionError(self,
404 asEyeByScalar = None,
405 asEyeByVector = None,
409 Permet de définir la covariance des erreurs de modèle :
410 - asCovariance : entrée des données, comme une matrice compatible avec
411 le constructeur de numpy.matrix
412 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
413 multiplicatif d'une matrice de corrélation identité, aucune matrice
414 n'étant donc explicitement à donner
415 - asEyeByVector : entrée des données comme un seul vecteur de variance,
416 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
417 n'étant donc explicitement à donner
418 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
419 être rendue disponible au même titre que les variables de calcul
421 if asEyeByScalar is not None:
422 self.__Q_scalar = float(asEyeByScalar)
424 elif asEyeByVector is not None:
425 self.__Q_scalar = None
426 self.__Q = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
427 elif asCovariance is not None:
428 self.__Q_scalar = None
429 self.__Q = numpy.matrix( asCovariance, numpy.float )
431 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")
433 self.__Parameters["Q_scalar"] = self.__Q_scalar
435 self.__StoredInputs["EvolutionError"] = self.__Q
438 # -----------------------------------------------------------
439 def setControlModel(self,
440 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
441 "useApproximatedDerivatives":False,
442 "withCenteredDF" :False,
443 "withIncrement" :0.01,
451 Permet de définir un opérateur de controle C. L'ordre de priorité des
452 définitions et leur sens sont les suivants :
453 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
454 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
455 "Direct" n'est pas définie, on prend la fonction "Tangent".
456 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
457 des opérateurs tangents et adjoints. On utilise par défaut des
458 différences finies non centrées ou centrées (si "withCenteredDF" est
459 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
460 du point courant ou sur le point fixe "withdX".
461 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
462 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
463 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
464 matrice fournie doit être sous une forme compatible avec le
465 constructeur de numpy.matrix.
466 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
467 être rendue disponible au même titre que les variables de calcul
469 if (type(asFunction) is type({})) and \
470 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
471 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
472 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
473 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
474 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
475 from daNumerics.ApproximatedDerivatives import FDApproximation
476 FDA = FDApproximation(
477 Function = asFunction["Direct"],
478 centeredDF = asFunction["withCenteredDF"],
479 increment = asFunction["withIncrement"],
480 dX = asFunction["withdX"] )
481 self.__CM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
482 self.__CM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
483 self.__CM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
484 elif (type(asFunction) is type({})) and \
485 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
486 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
487 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
488 self.__CM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
490 self.__CM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
491 self.__CM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
492 self.__CM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
493 elif asMatrix is not None:
494 matrice = numpy.matrix( asMatrix, numpy.float )
495 self.__CM["Direct"] = Operator( fromMatrix = matrice )
496 self.__CM["Tangent"] = Operator( fromMatrix = matrice )
497 self.__CM["Adjoint"] = Operator( fromMatrix = matrice.T )
500 raise ValueError("Improperly defined input control model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
503 self.__StoredInputs["ControlModel"] = self.__CM
506 def setControlInput(self,
508 asPersistentVector = None,
513 Permet de définir le controle en entree :
514 - asVector : entrée des données, comme un vecteur compatible avec le
515 constructeur de numpy.matrix
516 - asPersistentVector : entrée des données, comme un vecteur de type
517 persistent contruit avec la classe ad-hoc "Persistence"
518 - Scheduler est le contrôle temporel des données disponibles
519 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
520 être rendue disponible au même titre que les variables de calcul
522 if asVector is not None:
523 if isinstance(asVector,numpy.matrix):
524 self.__U = numpy.matrix( asVector.A1, numpy.float ).T
526 self.__U = numpy.matrix( asVector, numpy.float ).T
527 elif asPersistentVector is not None:
528 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
529 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.matrix)
530 for member in asPersistentVector:
531 self.__U.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
533 self.__U = asPersistentVector
535 raise ValueError("Error: improperly defined control input, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
537 self.__StoredInputs["ControlInput"] = self.__U
540 # -----------------------------------------------------------
541 def setControls (self,
546 Permet de définir la valeur initiale du vecteur X contenant toutes les
547 variables de contrôle, i.e. les paramètres ou l'état dont on veut
548 estimer la valeur pour obtenir les observations. C'est utile pour un
549 algorithme itératif/incrémental.
550 - asVector : entrée des données, comme un vecteur compatible avec le
551 constructeur de numpy.matrix.
552 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
553 être rendue disponible au même titre que les variables de calcul
555 if asVector is not None:
556 self.__X.store( asVector )
558 self.__StoredInputs["Controls"] = self.__X
561 # -----------------------------------------------------------
562 def setAlgorithm(self, choice = None ):
564 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
565 d'assimilation. L'argument est un champ caractère se rapportant au nom
566 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
567 d'assimilation sur les arguments fixes.
570 raise ValueError("Error: algorithm choice has to be given")
571 if self.__algorithmName is not None:
572 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
573 daDirectory = "daAlgorithms"
575 # Recherche explicitement le fichier complet
576 # ------------------------------------------
578 for directory in sys.path:
579 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
580 module_path = os.path.abspath(os.path.join(directory, daDirectory))
581 if module_path is None:
582 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
584 # Importe le fichier complet comme un module
585 # ------------------------------------------
587 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
588 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
589 self.__algorithmName = str(choice)
590 sys.path = sys_path_tmp ; del sys_path_tmp
591 except ImportError, e:
592 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
594 # Instancie un objet du type élémentaire du fichier
595 # -------------------------------------------------
596 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
597 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
600 def setAlgorithmParameters(self, asDico=None):
602 Permet de définir les paramètres de l'algorithme, sous la forme d'un
605 if asDico is not None:
606 self.__Parameters.update( dict( asDico ) )
608 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
611 def getAlgorithmParameters(self, noDetails=True):
613 Renvoie la liste des paramètres requis selon l'algorithme
615 return self.__algorithm.getRequiredParameters(noDetails)
617 # -----------------------------------------------------------
618 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
620 Permet de sélectionner un diagnostic a effectuer.
623 raise ValueError("Error: diagnostic choice has to be given")
624 daDirectory = "daDiagnostics"
626 # Recherche explicitement le fichier complet
627 # ------------------------------------------
629 for directory in sys.path:
630 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
631 module_path = os.path.abspath(os.path.join(directory, daDirectory))
632 if module_path is None:
633 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
635 # Importe le fichier complet comme un module
636 # ------------------------------------------
638 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
639 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
640 sys.path = sys_path_tmp ; del sys_path_tmp
641 except ImportError, e:
642 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
644 # Instancie un objet du type élémentaire du fichier
645 # -------------------------------------------------
646 if self.__StoredInputs.has_key(name):
647 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
648 elif self.__StoredDiagnostics.has_key(name):
649 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
651 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
655 parameters = parameters )
658 # -----------------------------------------------------------
659 def shape_validate(self):
661 Validation de la correspondance correcte des tailles des variables et
662 des matrices s'il y en a.
664 if self.__Xb is None: __Xb_shape = (0,)
665 elif hasattr(self.__Xb,"shape"):
666 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
667 else: __Xb_shape = self.__Xb.shape()
668 else: raise TypeError("Xb has no attribute of shape: problem !")
670 if self.__Y is None: __Y_shape = (0,)
671 elif hasattr(self.__Y,"shape"):
672 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
673 else: __Y_shape = self.__Y.shape()
674 else: raise TypeError("Y has no attribute of shape: problem !")
676 if self.__U is None: __U_shape = (0,)
677 elif hasattr(self.__U,"shape"):
678 if type(self.__U.shape) is tuple: __U_shape = self.__U.shape
679 else: __U_shape = self.__U.shape()
680 else: raise TypeError("U has no attribute of shape: problem !")
682 if self.__B is None and self.__B_scalar is None:
684 elif self.__B is None and self.__B_scalar is not None:
685 __B_shape = (max(__Xb_shape),max(__Xb_shape))
686 elif hasattr(self.__B,"shape"):
687 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
688 else: __B_shape = self.__B.shape()
689 else: raise TypeError("B has no attribute of shape: problem !")
691 if self.__R is None and self.__R_scalar is None:
693 elif self.__R is None and self.__R_scalar is not None:
694 __R_shape = (max(__Y_shape),max(__Y_shape))
695 elif hasattr(self.__R,"shape"):
696 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
697 else: __R_shape = self.__R.shape()
698 else: raise TypeError("R has no attribute of shape: problem !")
700 if self.__Q is None: __Q_shape = (0,0)
701 elif hasattr(self.__Q,"shape"):
702 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
703 else: __Q_shape = self.__Q.shape()
704 else: raise TypeError("Q has no attribute of shape: problem !")
706 if len(self.__HO) == 0: __HO_shape = (0,0)
707 elif type(self.__HO) is type({}): __HO_shape = (0,0)
708 elif hasattr(self.__HO["Direct"],"shape"):
709 if type(self.__HO["Direct"].shape) is tuple: __HO_shape = self.__HO["Direct"].shape
710 else: __HO_shape = self.__HO["Direct"].shape()
711 else: raise TypeError("H has no attribute of shape: problem !")
713 if len(self.__EM) == 0: __EM_shape = (0,0)
714 elif type(self.__EM) is type({}): __EM_shape = (0,0)
715 elif hasattr(self.__EM["Direct"],"shape"):
716 if type(self.__EM["Direct"].shape) is tuple: __EM_shape = self.__EM["Direct"].shape
717 else: __EM_shape = self.__EM["Direct"].shape()
718 else: raise TypeError("EM has no attribute of shape: problem !")
720 if len(self.__CM) == 0: __CM_shape = (0,0)
721 elif type(self.__CM) is type({}): __CM_shape = (0,0)
722 elif hasattr(self.__CM["Direct"],"shape"):
723 if type(self.__CM["Direct"].shape) is tuple: __CM_shape = self.__CM["Direct"].shape
724 else: __CM_shape = self.__CM["Direct"].shape()
725 else: raise TypeError("CM has no attribute of shape: problem !")
727 # Vérification des conditions
728 # ---------------------------
729 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
730 raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
731 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
732 raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
734 if not( min(__B_shape) == max(__B_shape) ):
735 raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
736 if not( min(__R_shape) == max(__R_shape) ):
737 raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
738 if not( min(__Q_shape) == max(__Q_shape) ):
739 raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
740 if not( min(__EM_shape) == max(__EM_shape) ):
741 raise ValueError("Shape characteristic of EM is incorrect: \"%s\""%(__EM_shape,))
743 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
744 raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__HO_shape,__Xb_shape))
745 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
746 raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__HO_shape,__Y_shape))
747 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and len(self.__B) > 0 and not( __HO_shape[1] == __B_shape[0] ):
748 raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__HO_shape,__B_shape))
749 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and len(self.__R) > 0 and not( __HO_shape[0] == __R_shape[1] ):
750 raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__HO_shape,__R_shape))
752 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
753 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
754 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
755 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
756 for member in asPersistentVector:
757 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
758 __Xb_shape = min(__B_shape)
760 raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
762 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
763 raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
765 if self.__EM is not None and len(self.__EM) > 0 and not(type(self.__EM) is type({})) and not( __EM_shape[1] == max(__Xb_shape) ):
766 raise ValueError("Shape characteristic of EM \"%s\" and X \"%s\" are incompatible"%(__EM_shape,__Xb_shape))
768 if self.__CM is not None and len(self.__CM) > 0 and not(type(self.__CM) is type({})) and not( __CM_shape[1] == max(__U_shape) ):
769 raise ValueError("Shape characteristic of CM \"%s\" and U \"%s\" are incompatible"%(__CM_shape,__U_shape))
773 # -----------------------------------------------------------
776 Permet de lancer le calcul d'assimilation.
778 Le nom de la méthode à activer est toujours "run". Les paramètres en
779 arguments de la méthode sont fixés. En sortie, on obtient les résultats
780 dans la variable de type dictionnaire "StoredVariables", qui contient en
781 particulier des objets de Persistence pour les analyses, OMA...
783 self.shape_validate()
785 self.__algorithm.run(
795 Parameters = self.__Parameters,
799 # -----------------------------------------------------------
800 def get(self, key=None):
802 Renvoie les résultats disponibles après l'exécution de la méthode
803 d'assimilation, ou les diagnostics disponibles. Attention, quand un
804 diagnostic porte le même nom qu'une variable stockée, c'est la variable
805 stockée qui est renvoyée, et le diagnostic est inatteignable.
808 if self.__algorithm.has_key(key):
809 return self.__algorithm.get( key )
810 elif self.__StoredInputs.has_key(key):
811 return self.__StoredInputs[key]
812 elif self.__StoredDiagnostics.has_key(key):
813 return self.__StoredDiagnostics[key]
815 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
817 allvariables = self.__algorithm.get()
818 allvariables.update( self.__StoredDiagnostics )
819 allvariables.update( self.__StoredInputs )
822 def get_available_variables(self):
824 Renvoie les variables potentiellement utilisables pour l'étude,
825 initialement stockées comme données d'entrées ou dans les algorithmes,
826 identifiés par les chaînes de caractères. L'algorithme doit avoir été
827 préalablement choisi sinon la méthode renvoie "None".
829 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
833 if len( self.__algorithm.keys()) > 0:
834 variables.extend( self.__algorithm.get().keys() )
835 if len( self.__StoredInputs.keys() ) > 0:
836 variables.extend( self.__StoredInputs.keys() )
840 def get_available_algorithms(self):
842 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
843 par les chaînes de caractères.
846 for directory in sys.path:
847 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
848 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
849 root, ext = os.path.splitext(fname)
850 if ext == '.py' and root != '__init__':
855 def get_available_diagnostics(self):
857 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
858 par les chaînes de caractères.
861 for directory in sys.path:
862 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
863 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
864 root, ext = os.path.splitext(fname)
865 if ext == '.py' and root != '__init__':
870 # -----------------------------------------------------------
871 def get_algorithms_main_path(self):
873 Renvoie le chemin pour le répertoire principal contenant les algorithmes
874 dans un sous-répertoire "daAlgorithms"
878 def add_algorithms_path(self, asPath=None):
880 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
881 se trouve un sous-répertoire "daAlgorithms"
883 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
884 pas indispensable de le rajouter ici.
886 if not os.path.isdir(asPath):
887 raise ValueError("The given "+asPath+" argument must exist as a directory")
888 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
889 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
890 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
891 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
892 sys.path.insert(0, os.path.abspath(asPath))
893 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
896 def get_diagnostics_main_path(self):
898 Renvoie le chemin pour le répertoire principal contenant les diagnostics
899 dans un sous-répertoire "daDiagnostics"
903 def add_diagnostics_path(self, asPath=None):
905 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
906 se trouve un sous-répertoire "daDiagnostics"
908 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
909 pas indispensable de le rajouter ici.
911 if not os.path.isdir(asPath):
912 raise ValueError("The given "+asPath+" argument must exist as a directory")
913 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
914 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
915 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
916 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
917 sys.path.insert(0, os.path.abspath(asPath))
918 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
921 # -----------------------------------------------------------
922 def setDataObserver(self,
925 HookParameters = None,
929 Permet d'associer un observer à une ou des variables nommées gérées en
930 interne, activable selon des règles définies dans le Scheduler. A chaque
931 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
932 les arguments (variable persistante VariableName, paramètres HookParameters).
935 if type( self.__algorithm ) is dict:
936 raise ValueError("No observer can be build before choosing an algorithm.")
938 # Vérification du nom de variable et typage
939 # -----------------------------------------
940 if type( VariableName ) is str:
941 VariableNames = [VariableName,]
942 elif type( VariableName ) is list:
943 VariableNames = map( str, VariableName )
945 raise ValueError("The observer requires a name or a list of names of variables.")
947 # Association interne de l'observer à la variable
948 # -----------------------------------------------
949 for n in VariableNames:
950 if not self.__algorithm.has_key( n ):
951 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
953 self.__algorithm.StoredVariables[ n ].setDataObserver(
954 Scheduler = Scheduler,
955 HookFunction = HookFunction,
956 HookParameters = HookParameters,
959 def removeDataObserver(self,
964 Permet de retirer un observer à une ou des variable nommée.
967 if type( self.__algorithm ) is dict:
968 raise ValueError("No observer can be removed before choosing an algorithm.")
970 # Vérification du nom de variable et typage
971 # -----------------------------------------
972 if type( VariableName ) is str:
973 VariableNames = [VariableName,]
974 elif type( VariableName ) is list:
975 VariableNames = map( str, VariableName )
977 raise ValueError("The observer requires a name or a list of names of variables.")
979 # Association interne de l'observer à la variable
980 # -----------------------------------------------
981 for n in VariableNames:
982 if not self.__algorithm.has_key( n ):
983 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
985 self.__algorithm.StoredVariables[ n ].removeDataObserver(
986 HookFunction = HookFunction,
989 # -----------------------------------------------------------
990 def setDebug(self, level=10):
992 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
993 appel pour changer le niveau de verbosité, avec :
994 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
997 log = logging.getLogger()
998 log.setLevel( level )
1000 def unsetDebug(self):
1002 Remet le logger au niveau par défaut
1005 log = logging.getLogger()
1006 log.setLevel( logging.WARNING )
1008 def prepare_to_pickle(self):
1009 self.__algorithmFile = None
1010 self.__diagnosticFile = None
1015 # ==============================================================================
1016 if __name__ == "__main__":
1017 print '\n AUTODIAGNOSTIC \n'
1019 ADD = AssimilationStudy("Ma premiere etude BLUE")
1021 ADD.setBackground (asVector = [0, 1, 2])
1022 ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
1023 ADD.setObservation (asVector = [0.5, 1.5, 2.5])
1024 ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
1025 ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
1027 ADD.setAlgorithm(choice="Blue")
1031 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1032 print "Ebauche :", [0, 1, 2]
1033 print "Observation :", [0.5, 1.5, 2.5]
1034 print "Demi-somme :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
1035 print " qui doit être identique à :"
1036 print "Analyse résultante :", ADD.get("Analysis")[0]
1039 print "Algorithmes disponibles.......................:", ADD.get_available_algorithms()
1040 # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path()
1041 print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics()
1042 # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path()
1043 print "Variables disponibles.........................:", ADD.get_available_variables()
1046 print "Paramètres requis par algorithme :"
1047 for algo in ADD.get_available_algorithms():
1048 tmpADD = AssimilationStudy("Un algorithme")
1049 tmpADD.setAlgorithm(choice=algo)
1050 print " %25s : %s"%(algo,tmpADD.getAlgorithmParameters())
1054 ADD.setDiagnostic("RMS", "Ma RMS")
1056 liste = ADD.get().keys()
1058 print "Variables et diagnostics nommés disponibles...:", liste
1061 print "Exemple de mise en place d'un observeur :"
1062 def obs(var=None,info=None):
1063 print " ---> Mise en oeuvre de l'observer"
1064 print " var =",var[-1]
1065 print " info =",info
1066 ADD.setDataObserver( 'Analysis', HookFunction=obs, Scheduler = [2, 4], HookParameters = "Second observer")
1067 # Attention, il faut décaler le stockage de 1 pour suivre le pas interne
1068 # car le pas 0 correspond à l'analyse ci-dessus.
1069 for i in range(1,6):
1071 print "Action sur la variable observée, étape :",i
1072 ADD.get('Analysis').store( [i, i, i] )
1075 print "Mise en debug et hors debug"
1076 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1080 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1082 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()