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
36 from BasicObjects import Operator
37 from PlatformInfo import uniq
39 # ==============================================================================
40 class AssimilationStudy:
42 Cette classe sert d'interface pour l'utilisation de l'assimilation de
43 données. Elle contient les méthodes ou accesseurs nécessaires à la
44 construction d'un calcul d'assimilation.
46 def __init__(self, name=""):
48 Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
49 élémentaire. Ces variables sont ensuite disponibles pour implémenter un
50 algorithme élémentaire particulier.
52 Background............: vecteur Xb
53 Observation...........: vecteur Y (potentiellement temporel)
55 State.................: vecteur d'état dont une partie est le vecteur de
56 contrôle. Cette information n'est utile que si l'on veut faire des
57 calculs sur l'état complet, mais elle n'est pas indispensable pour
59 Control...............: vecteur X contenant toutes les variables de
60 contrôle, i.e. les paramètres ou l'état dont on veut estimer la
61 valeur pour obtenir les observations
62 ObservationOperator...: opérateur d'observation H
64 Les observations présentent une erreur dont la matrice de covariance est
65 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
68 self.__name = str(name)
79 self.__X = Persistence.OneVector()
80 self.__Parameters = {}
81 self.__StoredDiagnostics = {}
82 self.__StoredInputs = {}
84 self.__B_scalar = None
85 self.__R_scalar = None
86 self.__Q_scalar = None
87 self.__Parameters["B_scalar"] = None
88 self.__Parameters["R_scalar"] = None
89 self.__Parameters["Q_scalar"] = None
91 # Variables temporaires
93 self.__algorithmFile = None
94 self.__algorithmName = None
95 self.__diagnosticFile = None
97 # Récupère le chemin du répertoire parent et l'ajoute au path
98 # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
99 # qui est activée dans Persistence)
100 self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
101 sys.path.insert(0, self.__parent)
102 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
104 # ---------------------------------------------------------
105 def setBackground(self,
107 asPersistentVector = None,
112 Permet de définir l'estimation a priori :
113 - asVector : entrée des données, comme un vecteur compatible avec le
114 constructeur de numpy.matrix
115 - asPersistentVector : entrée des données, comme un vecteur de type
116 persistent contruit avec la classe ad-hoc "Persistence"
117 - Scheduler est le contrôle temporel des données
118 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
119 être rendue disponible au même titre que les variables de calcul
121 if asVector is not None:
122 if type( asVector ) is type( numpy.matrix([]) ):
123 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
125 self.__Xb = numpy.matrix( asVector, numpy.float ).T
126 elif asPersistentVector is not None:
127 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
128 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
129 for member in asPersistentVector:
130 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
132 self.__Xb = asPersistentVector
134 raise ValueError("Error: improperly defined background, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
136 self.__StoredInputs["Background"] = self.__Xb
139 def setBackgroundError(self,
141 asEyeByScalar = None,
142 asEyeByVector = None,
146 Permet de définir la covariance des erreurs d'ébauche :
147 - asCovariance : entrée des données, comme une matrice compatible avec
148 le constructeur de numpy.matrix
149 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
150 multiplicatif d'une matrice de corrélation identité, aucune matrice
151 n'étant donc explicitement à donner
152 - asEyeByVector : entrée des données comme un seul vecteur de variance,
153 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
154 n'étant donc explicitement à donner
155 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
156 être rendue disponible au même titre que les variables de calcul
158 if asEyeByScalar is not None:
159 self.__B_scalar = float(asEyeByScalar)
161 elif asEyeByVector is not None:
162 self.__B_scalar = None
163 self.__B = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
164 elif asCovariance is not None:
165 self.__B_scalar = None
166 self.__B = numpy.matrix( asCovariance, numpy.float )
168 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")
170 self.__Parameters["B_scalar"] = self.__B_scalar
172 self.__StoredInputs["BackgroundError"] = self.__B
175 # -----------------------------------------------------------
176 def setObservation(self,
178 asPersistentVector = None,
183 Permet de définir les observations :
184 - asVector : entrée des données, comme un vecteur compatible avec le
185 constructeur de numpy.matrix
186 - asPersistentVector : entrée des données, comme un vecteur de type
187 persistent contruit avec la classe ad-hoc "Persistence"
188 - Scheduler est le contrôle temporel des données disponibles
189 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
190 être rendue disponible au même titre que les variables de calcul
192 if asVector is not None:
193 if type( asVector ) is type( numpy.matrix([]) ):
194 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
196 self.__Y = numpy.matrix( asVector, numpy.float ).T
197 elif asPersistentVector is not None:
198 if type( asPersistentVector ) is list or type( asPersistentVector ) is tuple:
199 self.__Y = Persistence.OneVector("Observation", basetype=numpy.array)
200 for y in asPersistentVector:
203 self.__Y = asPersistentVector
205 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
207 self.__StoredInputs["Observation"] = self.__Y
210 def setObservationError(self,
212 asEyeByScalar = None,
213 asEyeByVector = None,
217 Permet de définir la covariance des erreurs d'observations :
218 - asCovariance : entrée des données, comme une matrice compatible avec
219 le constructeur de numpy.matrix
220 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
221 multiplicatif d'une matrice de corrélation identité, aucune matrice
222 n'étant donc explicitement à donner
223 - asEyeByVector : entrée des données comme un seul vecteur de variance,
224 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
225 n'étant donc explicitement à donner
226 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
227 être rendue disponible au même titre que les variables de calcul
229 if asEyeByScalar is not None:
230 self.__R_scalar = float(asEyeByScalar)
232 elif asEyeByVector is not None:
233 self.__R_scalar = None
234 self.__R = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
235 elif asCovariance is not None:
236 self.__R_scalar = None
237 self.__R = numpy.matrix( asCovariance, numpy.float )
239 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")
241 self.__Parameters["R_scalar"] = self.__R_scalar
243 self.__StoredInputs["ObservationError"] = self.__R
246 def setObservationOperator(self,
247 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
248 "useApproximatedDerivatives":False,
249 "withCenteredDF" :False,
250 "withIncrement" :0.01,
258 Permet de définir un opérateur d'observation H. L'ordre de priorité des
259 définitions et leur sens sont les suivants :
260 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
261 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
262 "Direct" n'est pas définie, on prend la fonction "Tangent".
263 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
264 des opérateurs tangents et adjoints. On utilise par défaut des
265 différences finies non centrées ou centrées (si "withCenteredDF" est
266 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
267 du point courant ou sur le point fixe "withdX".
268 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
269 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
270 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
271 matrice fournie doit être sous une forme compatible avec le
272 constructeur de numpy.matrix.
273 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
274 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
275 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
276 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
277 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
278 être rendue disponible au même titre que les variables de calcul
280 if (type(asFunction) is type({})) and \
281 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
282 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
283 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
284 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
285 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
286 from daNumerics.ApproximatedDerivatives import FDApproximation
287 FDA = FDApproximation(
288 Function = asFunction["Direct"],
289 centeredDF = asFunction["withCenteredDF"],
290 increment = asFunction["withIncrement"],
291 dX = asFunction["withdX"] )
292 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator )
293 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
294 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
295 elif (type(asFunction) is type({})) and \
296 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
297 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
298 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
299 self.__HO["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
301 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"] )
302 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
303 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
304 elif asMatrix is not None:
305 matrice = numpy.matrix( asMatrix, numpy.float )
306 self.__HO["Direct"] = Operator( fromMatrix = matrice )
307 self.__HO["Tangent"] = Operator( fromMatrix = matrice )
308 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T )
311 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
313 if appliedToX is not None:
314 self.__HO["AppliedToX"] = {}
315 if type(appliedToX) is not dict:
316 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
317 for key in appliedToX.keys():
318 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
319 # Pour le cas où l'on a une vraie matrice
320 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
321 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
322 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
323 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
325 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
327 self.__HO["AppliedToX"] = None
330 self.__StoredInputs["ObservationOperator"] = self.__HO
333 # -----------------------------------------------------------
334 def setEvolutionModel(self,
335 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
336 "useApproximatedDerivatives":False,
337 "withCenteredDF" :False,
338 "withIncrement" :0.01,
346 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
347 définitions et leur sens sont les suivants :
348 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
349 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
350 "Direct" n'est pas définie, on prend la fonction "Tangent".
351 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
352 des opérateurs tangents et adjoints. On utilise par défaut des
353 différences finies non centrées ou centrées (si "withCenteredDF" est
354 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
355 du point courant ou sur le point fixe "withdX".
356 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
357 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
358 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
359 matrice fournie doit être sous une forme compatible avec le
360 constructeur de numpy.matrix.
361 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
362 être rendue disponible au même titre que les variables de calcul
364 if (type(asFunction) is type({})) and \
365 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
366 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
367 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
368 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
369 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
370 from daNumerics.ApproximatedDerivatives import FDApproximation
371 FDA = FDApproximation(
372 Function = asFunction["Direct"],
373 centeredDF = asFunction["withCenteredDF"],
374 increment = asFunction["withIncrement"],
375 dX = asFunction["withdX"] )
376 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
377 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
378 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
379 elif (type(asFunction) is type({})) and \
380 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
381 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
382 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
383 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
385 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
386 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
387 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
388 elif asMatrix is not None:
389 matrice = numpy.matrix( asMatrix, numpy.float )
390 self.__EM["Direct"] = Operator( fromMatrix = matrice )
391 self.__EM["Tangent"] = Operator( fromMatrix = matrice )
392 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T )
395 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
398 self.__StoredInputs["EvolutionModel"] = self.__EM
401 def setEvolutionError(self,
403 asEyeByScalar = None,
404 asEyeByVector = None,
408 Permet de définir la covariance des erreurs de modèle :
409 - asCovariance : entrée des données, comme une matrice compatible avec
410 le constructeur de numpy.matrix
411 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
412 multiplicatif d'une matrice de corrélation identité, aucune matrice
413 n'étant donc explicitement à donner
414 - asEyeByVector : entrée des données comme un seul vecteur de variance,
415 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
416 n'étant donc explicitement à donner
417 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
418 être rendue disponible au même titre que les variables de calcul
420 if asEyeByScalar is not None:
421 self.__Q_scalar = float(asEyeByScalar)
423 elif asEyeByVector is not None:
424 self.__Q_scalar = None
425 self.__Q = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
426 elif asCovariance is not None:
427 self.__Q_scalar = None
428 self.__Q = numpy.matrix( asCovariance, numpy.float )
430 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")
432 self.__Parameters["Q_scalar"] = self.__Q_scalar
434 self.__StoredInputs["EvolutionError"] = self.__Q
437 # -----------------------------------------------------------
438 def setControlModel(self,
439 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
440 "useApproximatedDerivatives":False,
441 "withCenteredDF" :False,
442 "withIncrement" :0.01,
450 Permet de définir un opérateur de controle C. L'ordre de priorité des
451 définitions et leur sens sont les suivants :
452 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
453 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
454 "Direct" n'est pas définie, on prend la fonction "Tangent".
455 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
456 des opérateurs tangents et adjoints. On utilise par défaut des
457 différences finies non centrées ou centrées (si "withCenteredDF" est
458 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
459 du point courant ou sur le point fixe "withdX".
460 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
461 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
462 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
463 matrice fournie doit être sous une forme compatible avec le
464 constructeur de numpy.matrix.
465 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
466 être rendue disponible au même titre que les variables de calcul
468 if (type(asFunction) is type({})) and \
469 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
470 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
471 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
472 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
473 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
474 from daNumerics.ApproximatedDerivatives import FDApproximation
475 FDA = FDApproximation(
476 Function = asFunction["Direct"],
477 centeredDF = asFunction["withCenteredDF"],
478 increment = asFunction["withIncrement"],
479 dX = asFunction["withdX"] )
480 self.__CM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
481 self.__CM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
482 self.__CM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
483 elif (type(asFunction) is type({})) and \
484 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
485 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
486 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
487 self.__CM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
489 self.__CM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
490 self.__CM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
491 self.__CM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
492 elif asMatrix is not None:
493 matrice = numpy.matrix( asMatrix, numpy.float )
494 self.__CM["Direct"] = Operator( fromMatrix = matrice )
495 self.__CM["Tangent"] = Operator( fromMatrix = matrice )
496 self.__CM["Adjoint"] = Operator( fromMatrix = matrice.T )
499 raise ValueError("Improperly defined input control model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
502 self.__StoredInputs["ControlModel"] = self.__CM
505 def setControlInput(self,
507 asPersistentVector = None,
512 Permet de définir le controle en entree :
513 - asVector : entrée des données, comme un vecteur compatible avec le
514 constructeur de numpy.matrix
515 - asPersistentVector : entrée des données, comme un vecteur de type
516 persistent contruit avec la classe ad-hoc "Persistence"
517 - Scheduler est le contrôle temporel des données disponibles
518 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
519 être rendue disponible au même titre que les variables de calcul
521 if asVector is not None:
522 if isinstance(asVector,numpy.matrix):
523 self.__U = numpy.matrix( asVector.A1, numpy.float ).T
525 self.__U = numpy.matrix( asVector, numpy.float ).T
526 elif asPersistentVector is not None:
527 if isinstance(asPersistentVector,list) or isinstance( asPersistentVector,tuple):
528 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.array)
529 for y in asPersistentVector:
532 self.__U = asPersistentVector
534 raise ValueError("Error: improperly defined control input, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
536 self.__StoredInputs["ControlInput"] = self.__U
539 # -----------------------------------------------------------
540 def setControls (self,
545 Permet de définir la valeur initiale du vecteur X contenant toutes les
546 variables de contrôle, i.e. les paramètres ou l'état dont on veut
547 estimer la valeur pour obtenir les observations. C'est utile pour un
548 algorithme itératif/incrémental.
549 - asVector : entrée des données, comme un vecteur compatible avec le
550 constructeur de numpy.matrix.
551 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
552 être rendue disponible au même titre que les variables de calcul
554 if asVector is not None:
555 self.__X.store( asVector )
557 self.__StoredInputs["Controls"] = self.__X
560 # -----------------------------------------------------------
561 def setAlgorithm(self, choice = None ):
563 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
564 d'assimilation. L'argument est un champ caractère se rapportant au nom
565 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
566 d'assimilation sur les arguments fixes.
569 raise ValueError("Error: algorithm choice has to be given")
570 if self.__algorithmName is not None:
571 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
572 daDirectory = "daAlgorithms"
574 # Recherche explicitement le fichier complet
575 # ------------------------------------------
577 for directory in sys.path:
578 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
579 module_path = os.path.abspath(os.path.join(directory, daDirectory))
580 if module_path is None:
581 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
583 # Importe le fichier complet comme un module
584 # ------------------------------------------
586 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
587 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
588 self.__algorithmName = str(choice)
589 sys.path = sys_path_tmp ; del sys_path_tmp
590 except ImportError, e:
591 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
593 # Instancie un objet du type élémentaire du fichier
594 # -------------------------------------------------
595 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
596 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
599 def setAlgorithmParameters(self, asDico=None):
601 Permet de définir les paramètres de l'algorithme, sous la forme d'un
604 if asDico is not None:
605 self.__Parameters.update( dict( asDico ) )
607 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
610 def getAlgorithmParameters(self, noDetails=True):
612 Renvoie la liste des paramètres requis selon l'algorithme
614 return self.__algorithm.getRequiredParameters(noDetails)
616 # -----------------------------------------------------------
617 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
619 Permet de sélectionner un diagnostic a effectuer.
622 raise ValueError("Error: diagnostic choice has to be given")
623 daDirectory = "daDiagnostics"
625 # Recherche explicitement le fichier complet
626 # ------------------------------------------
628 for directory in sys.path:
629 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
630 module_path = os.path.abspath(os.path.join(directory, daDirectory))
631 if module_path is None:
632 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
634 # Importe le fichier complet comme un module
635 # ------------------------------------------
637 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
638 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
639 sys.path = sys_path_tmp ; del sys_path_tmp
640 except ImportError, e:
641 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
643 # Instancie un objet du type élémentaire du fichier
644 # -------------------------------------------------
645 if self.__StoredInputs.has_key(name):
646 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
647 elif self.__StoredDiagnostics.has_key(name):
648 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
650 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
654 parameters = parameters )
657 # -----------------------------------------------------------
658 def shape_validate(self):
660 Validation de la correspondance correcte des tailles des variables et
661 des matrices s'il y en a.
663 if self.__Xb is None: __Xb_shape = (0,)
664 elif hasattr(self.__Xb,"shape"):
665 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
666 else: __Xb_shape = self.__Xb.shape()
667 else: raise TypeError("Xb has no attribute of shape: problem !")
669 if self.__Y is None: __Y_shape = (0,)
670 elif hasattr(self.__Y,"shape"):
671 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
672 else: __Y_shape = self.__Y.shape()
673 else: raise TypeError("Y has no attribute of shape: problem !")
675 if self.__U is None: __U_shape = (0,)
676 elif hasattr(self.__U,"shape"):
677 if type(self.__U.shape) is tuple: __U_shape = self.__U.shape
678 else: __U_shape = self.__U.shape()
679 else: raise TypeError("U has no attribute of shape: problem !")
681 if self.__B is None and self.__B_scalar is None:
683 elif self.__B is None and self.__B_scalar is not None:
684 __B_shape = (max(__Xb_shape),max(__Xb_shape))
685 elif hasattr(self.__B,"shape"):
686 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
687 else: __B_shape = self.__B.shape()
688 else: raise TypeError("B has no attribute of shape: problem !")
690 if self.__R is None and self.__R_scalar is None:
692 elif self.__R is None and self.__R_scalar is not None:
693 __R_shape = (max(__Y_shape),max(__Y_shape))
694 elif hasattr(self.__R,"shape"):
695 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
696 else: __R_shape = self.__R.shape()
697 else: raise TypeError("R has no attribute of shape: problem !")
699 if self.__Q is None: __Q_shape = (0,0)
700 elif hasattr(self.__Q,"shape"):
701 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
702 else: __Q_shape = self.__Q.shape()
703 else: raise TypeError("Q has no attribute of shape: problem !")
705 if len(self.__HO) == 0: __HO_shape = (0,0)
706 elif type(self.__HO) is type({}): __HO_shape = (0,0)
707 elif hasattr(self.__HO["Direct"],"shape"):
708 if type(self.__HO["Direct"].shape) is tuple: __HO_shape = self.__HO["Direct"].shape
709 else: __HO_shape = self.__HO["Direct"].shape()
710 else: raise TypeError("H has no attribute of shape: problem !")
712 if len(self.__EM) == 0: __EM_shape = (0,0)
713 elif type(self.__EM) is type({}): __EM_shape = (0,0)
714 elif hasattr(self.__EM["Direct"],"shape"):
715 if type(self.__EM["Direct"].shape) is tuple: __EM_shape = self.__EM["Direct"].shape
716 else: __EM_shape = self.__EM["Direct"].shape()
717 else: raise TypeError("EM has no attribute of shape: problem !")
719 if len(self.__CM) == 0: __CM_shape = (0,0)
720 elif type(self.__CM) is type({}): __CM_shape = (0,0)
721 elif hasattr(self.__CM["Direct"],"shape"):
722 if type(self.__CM["Direct"].shape) is tuple: __CM_shape = self.__CM["Direct"].shape
723 else: __CM_shape = self.__CM["Direct"].shape()
724 else: raise TypeError("CM has no attribute of shape: problem !")
726 # Vérification des conditions
727 # ---------------------------
728 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
729 raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
730 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
731 raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
733 if not( min(__B_shape) == max(__B_shape) ):
734 raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
735 if not( min(__R_shape) == max(__R_shape) ):
736 raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
737 if not( min(__Q_shape) == max(__Q_shape) ):
738 raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
739 if not( min(__EM_shape) == max(__EM_shape) ):
740 raise ValueError("Shape characteristic of EM is incorrect: \"%s\""%(__EM_shape,))
742 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
743 raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__HO_shape,__Xb_shape))
744 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
745 raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__HO_shape,__Y_shape))
746 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] ):
747 raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__HO_shape,__B_shape))
748 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] ):
749 raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__HO_shape,__R_shape))
751 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
752 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
753 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
754 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
755 for member in asPersistentVector:
756 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
757 __Xb_shape = min(__B_shape)
759 raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
761 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
762 raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
764 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) ):
765 raise ValueError("Shape characteristic of EM \"%s\" and X \"%s\" are incompatible"%(__EM_shape,__Xb_shape))
767 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) ):
768 raise ValueError("Shape characteristic of CM \"%s\" and U \"%s\" are incompatible"%(__CM_shape,__U_shape))
772 # -----------------------------------------------------------
775 Permet de lancer le calcul d'assimilation.
777 Le nom de la méthode à activer est toujours "run". Les paramètres en
778 arguments de la méthode sont fixés. En sortie, on obtient les résultats
779 dans la variable de type dictionnaire "StoredVariables", qui contient en
780 particulier des objets de Persistence pour les analyses, OMA...
782 self.shape_validate()
784 self.__algorithm.run(
794 Parameters = self.__Parameters,
798 # -----------------------------------------------------------
799 def get(self, key=None):
801 Renvoie les résultats disponibles après l'exécution de la méthode
802 d'assimilation, ou les diagnostics disponibles. Attention, quand un
803 diagnostic porte le même nom qu'une variable stockée, c'est la variable
804 stockée qui est renvoyée, et le diagnostic est inatteignable.
807 if self.__algorithm.has_key(key):
808 return self.__algorithm.get( key )
809 elif self.__StoredInputs.has_key(key):
810 return self.__StoredInputs[key]
811 elif self.__StoredDiagnostics.has_key(key):
812 return self.__StoredDiagnostics[key]
814 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
816 allvariables = self.__algorithm.get()
817 allvariables.update( self.__StoredDiagnostics )
818 allvariables.update( self.__StoredInputs )
821 def get_available_variables(self):
823 Renvoie les variables potentiellement utilisables pour l'étude,
824 initialement stockées comme données d'entrées ou dans les algorithmes,
825 identifiés par les chaînes de caractères. L'algorithme doit avoir été
826 préalablement choisi sinon la méthode renvoie "None".
828 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
832 if len( self.__algorithm.keys()) > 0:
833 variables.extend( self.__algorithm.get().keys() )
834 if len( self.__StoredInputs.keys() ) > 0:
835 variables.extend( self.__StoredInputs.keys() )
839 def get_available_algorithms(self):
841 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
842 par les chaînes de caractères.
845 for directory in sys.path:
846 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
847 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
848 root, ext = os.path.splitext(fname)
849 if ext == '.py' and root != '__init__':
854 def get_available_diagnostics(self):
856 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
857 par les chaînes de caractères.
860 for directory in sys.path:
861 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
862 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
863 root, ext = os.path.splitext(fname)
864 if ext == '.py' and root != '__init__':
869 # -----------------------------------------------------------
870 def get_algorithms_main_path(self):
872 Renvoie le chemin pour le répertoire principal contenant les algorithmes
873 dans un sous-répertoire "daAlgorithms"
877 def add_algorithms_path(self, asPath=None):
879 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
880 se trouve un sous-répertoire "daAlgorithms"
882 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
883 pas indispensable de le rajouter ici.
885 if not os.path.isdir(asPath):
886 raise ValueError("The given "+asPath+" argument must exist as a directory")
887 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
888 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
889 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
890 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
891 sys.path.insert(0, os.path.abspath(asPath))
892 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
895 def get_diagnostics_main_path(self):
897 Renvoie le chemin pour le répertoire principal contenant les diagnostics
898 dans un sous-répertoire "daDiagnostics"
902 def add_diagnostics_path(self, asPath=None):
904 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
905 se trouve un sous-répertoire "daDiagnostics"
907 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
908 pas indispensable de le rajouter ici.
910 if not os.path.isdir(asPath):
911 raise ValueError("The given "+asPath+" argument must exist as a directory")
912 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
913 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
914 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
915 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
916 sys.path.insert(0, os.path.abspath(asPath))
917 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
920 # -----------------------------------------------------------
921 def setDataObserver(self,
924 HookParameters = None,
928 Permet d'associer un observer à une ou des variables nommées gérées en
929 interne, activable selon des règles définies dans le Scheduler. A chaque
930 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
931 les arguments (variable persistante VariableName, paramètres HookParameters).
934 if type( self.__algorithm ) is dict:
935 raise ValueError("No observer can be build before choosing an algorithm.")
937 # Vérification du nom de variable et typage
938 # -----------------------------------------
939 if type( VariableName ) is str:
940 VariableNames = [VariableName,]
941 elif type( VariableName ) is list:
942 VariableNames = map( str, VariableName )
944 raise ValueError("The observer requires a name or a list of names of variables.")
946 # Association interne de l'observer à la variable
947 # -----------------------------------------------
948 for n in VariableNames:
949 if not self.__algorithm.has_key( n ):
950 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
952 self.__algorithm.StoredVariables[ n ].setDataObserver(
953 Scheduler = Scheduler,
954 HookFunction = HookFunction,
955 HookParameters = HookParameters,
958 def removeDataObserver(self,
963 Permet de retirer un observer à une ou des variable nommée.
966 if type( self.__algorithm ) is dict:
967 raise ValueError("No observer can be removed before choosing an algorithm.")
969 # Vérification du nom de variable et typage
970 # -----------------------------------------
971 if type( VariableName ) is str:
972 VariableNames = [VariableName,]
973 elif type( VariableName ) is list:
974 VariableNames = map( str, VariableName )
976 raise ValueError("The observer requires a name or a list of names of variables.")
978 # Association interne de l'observer à la variable
979 # -----------------------------------------------
980 for n in VariableNames:
981 if not self.__algorithm.has_key( n ):
982 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
984 self.__algorithm.StoredVariables[ n ].removeDataObserver(
985 HookFunction = HookFunction,
988 # -----------------------------------------------------------
989 def setDebug(self, level=10):
991 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
992 appel pour changer le niveau de verbosité, avec :
993 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
996 log = logging.getLogger()
997 log.setLevel( level )
999 def unsetDebug(self):
1001 Remet le logger au niveau par défaut
1004 log = logging.getLogger()
1005 log.setLevel( logging.WARNING )
1007 def prepare_to_pickle(self):
1008 self.__algorithmFile = None
1009 self.__diagnosticFile = None
1014 # ==============================================================================
1015 if __name__ == "__main__":
1016 print '\n AUTODIAGNOSTIC \n'
1018 ADD = AssimilationStudy("Ma premiere etude BLUE")
1020 ADD.setBackground (asVector = [0, 1, 2])
1021 ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
1022 ADD.setObservation (asVector = [0.5, 1.5, 2.5])
1023 ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
1024 ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
1026 ADD.setAlgorithm(choice="Blue")
1030 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1031 print "Ebauche :", [0, 1, 2]
1032 print "Observation :", [0.5, 1.5, 2.5]
1033 print "Demi-somme :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
1034 print " qui doit être identique à :"
1035 print "Analyse résultante :", ADD.get("Analysis")[0]
1038 print "Algorithmes disponibles.......................:", ADD.get_available_algorithms()
1039 # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path()
1040 print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics()
1041 # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path()
1042 print "Variables disponibles.........................:", ADD.get_available_variables()
1045 print "Paramètres requis par algorithme :"
1046 for algo in ADD.get_available_algorithms():
1047 tmpADD = AssimilationStudy("Un algorithme")
1048 tmpADD.setAlgorithm(choice=algo)
1049 print " %25s : %s"%(algo,tmpADD.getAlgorithmParameters())
1053 ADD.setDiagnostic("RMS", "Ma RMS")
1055 liste = ADD.get().keys()
1057 print "Variables et diagnostics nommés disponibles...:", liste
1060 print "Exemple de mise en place d'un observeur :"
1061 def obs(var=None,info=None):
1062 print " ---> Mise en oeuvre de l'observer"
1063 print " var =",var[-1]
1064 print " info =",info
1065 ADD.setDataObserver( 'Analysis', HookFunction=obs, Scheduler = [2, 4], HookParameters = "Second observer")
1066 # Attention, il faut décaler le stockage de 1 pour suivre le pas interne
1067 # car le pas 0 correspond à l'analyse ci-dessus.
1068 for i in range(1,6):
1070 print "Action sur la variable observée, étape :",i
1071 ADD.get('Analysis').store( [i, i, i] )
1074 print "Mise en debug et hors debug"
1075 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1079 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1081 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()