1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2014 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 ExtendedLogging ; ExtendedLogging.ExtendedLogging() # A importer en premier
40 from BasicObjects import Operator, Covariance
41 from PlatformInfo import uniq
43 # ==============================================================================
44 class AssimilationStudy:
46 Cette classe sert d'interface pour l'utilisation de l'assimilation de
47 données. Elle contient les méthodes ou accesseurs nécessaires à la
48 construction d'un calcul d'assimilation.
50 def __init__(self, name=""):
52 Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
53 élémentaire. Ces variables sont ensuite disponibles pour implémenter un
54 algorithme élémentaire particulier.
56 - Background : vecteur Xb
57 - Observation : vecteur Y (potentiellement temporel) d'observations
58 - State : vecteur d'état dont une partie est le vecteur de contrôle.
59 Cette information n'est utile que si l'on veut faire des calculs sur
60 l'état complet, mais elle n'est pas indispensable pour l'assimilation.
61 - Control : vecteur X contenant toutes les variables de contrôle, i.e.
62 les paramètres ou l'état dont on veut estimer la valeur pour obtenir
64 - ObservationOperator...: opérateur d'observation H
66 Les observations présentent une erreur dont la matrice de covariance est
67 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
70 self.__name = str(name)
81 self.__X = Persistence.OneVector()
82 self.__Parameters = {}
83 self.__StoredDiagnostics = {}
84 self.__StoredInputs = {}
86 # Variables temporaires
88 self.__algorithmFile = None
89 self.__algorithmName = None
90 self.__diagnosticFile = None
92 # Récupère le chemin du répertoire parent et l'ajoute au path
93 # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
94 # qui est activée dans Persistence)
95 self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
96 sys.path.insert(0, self.__parent)
97 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
99 # ---------------------------------------------------------
100 def setBackground(self,
102 asPersistentVector = None,
107 Permet de définir l'estimation a priori :
108 - asVector : entrée des données, comme un vecteur compatible avec le
109 constructeur de numpy.matrix
110 - asPersistentVector : entrée des données, comme un vecteur de type
111 persistent contruit avec la classe ad-hoc "Persistence"
112 - Scheduler est le contrôle temporel des données
113 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
114 être rendue disponible au même titre que les variables de calcul
116 if asVector is not None:
117 if type( asVector ) is type( numpy.matrix([]) ):
118 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
120 self.__Xb = numpy.matrix( asVector, numpy.float ).T
121 elif asPersistentVector is not None:
122 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
123 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
124 for member in asPersistentVector:
125 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
127 self.__Xb = asPersistentVector
129 raise ValueError("Error: improperly defined background, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
131 self.__StoredInputs["Background"] = self.__Xb
134 def setBackgroundError(self,
136 asEyeByScalar = None,
137 asEyeByVector = None,
141 Permet de définir la covariance des erreurs d'ébauche :
142 - asCovariance : entrée des données, comme une matrice compatible avec
143 le constructeur de numpy.matrix
144 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
145 multiplicatif d'une matrice de corrélation identité, aucune matrice
146 n'étant donc explicitement à donner
147 - asEyeByVector : entrée des données comme un seul vecteur de variance,
148 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
149 n'étant donc explicitement à donner
150 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
151 être rendue disponible au même titre que les variables de calcul
153 self.__B = Covariance(
154 name = "BackgroundError",
155 asCovariance = asCovariance,
156 asEyeByScalar = asEyeByScalar,
157 asEyeByVector = asEyeByVector,
160 self.__StoredInputs["BackgroundError"] = self.__B
163 # -----------------------------------------------------------
164 def setObservation(self,
166 asPersistentVector = None,
171 Permet de définir les observations :
172 - asVector : entrée des données, comme un vecteur compatible avec le
173 constructeur de numpy.matrix
174 - asPersistentVector : entrée des données, comme un vecteur de type
175 persistent contruit avec la classe ad-hoc "Persistence"
176 - Scheduler est le contrôle temporel des données disponibles
177 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
178 être rendue disponible au même titre que les variables de calcul
180 if asVector is not None:
181 if type( asVector ) is type( numpy.matrix([]) ):
182 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
184 self.__Y = numpy.matrix( asVector, numpy.float ).T
185 elif asPersistentVector is not None:
186 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
187 self.__Y = Persistence.OneVector("Observation", basetype=numpy.matrix)
188 for member in asPersistentVector:
189 self.__Y.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
191 self.__Y = asPersistentVector
193 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
195 self.__StoredInputs["Observation"] = self.__Y
198 def setObservationError(self,
200 asEyeByScalar = None,
201 asEyeByVector = None,
205 Permet de définir la covariance des erreurs d'observations :
206 - asCovariance : entrée des données, comme une matrice compatible avec
207 le constructeur de numpy.matrix
208 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
209 multiplicatif d'une matrice de corrélation identité, aucune matrice
210 n'étant donc explicitement à donner
211 - asEyeByVector : entrée des données comme un seul vecteur de variance,
212 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
213 n'étant donc explicitement à donner
214 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
215 être rendue disponible au même titre que les variables de calcul
217 self.__R = Covariance(
218 name = "ObservationError",
219 asCovariance = asCovariance,
220 asEyeByScalar = asEyeByScalar,
221 asEyeByVector = asEyeByVector,
224 self.__StoredInputs["ObservationError"] = self.__R
227 def setObservationOperator(self,
234 Permet de définir un opérateur d'observation H. L'ordre de priorité des
235 définitions et leur sens sont les suivants :
236 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
237 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
238 "Direct" n'est pas définie, on prend la fonction "Tangent".
239 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
240 des opérateurs tangents et adjoints. On utilise par défaut des
241 différences finies non centrées ou centrées (si "withCenteredDF" est
242 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
243 du point courant ou sur le point fixe "withdX".
244 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
245 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
246 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
247 matrice fournie doit être sous une forme compatible avec le
248 constructeur de numpy.matrix.
249 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
250 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
251 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
252 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
253 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
254 être rendue disponible au même titre que les variables de calcul
255 L'argument "asFunction" peut prendre la forme complète suivante, avec
256 les valeurs par défaut standards :
257 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
258 "useApproximatedDerivatives":False,
259 "withCenteredDF" :False,
260 "withIncrement" :0.01,
262 "withAvoidingRedundancy" :True,
263 "withToleranceInRedundancy" :1.e-18,
264 "withLenghtOfRedundancy" :-1,
267 if (type(asFunction) is type({})) and \
268 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
269 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
270 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
271 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
272 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
273 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
274 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
275 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
276 from daNumerics.ApproximatedDerivatives import FDApproximation
277 FDA = FDApproximation(
278 Function = asFunction["Direct"],
279 centeredDF = asFunction["withCenteredDF"],
280 increment = asFunction["withIncrement"],
281 dX = asFunction["withdX"],
282 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
283 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
284 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
286 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator )
287 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
288 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
289 elif (type(asFunction) is type({})) and \
290 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
291 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
292 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
293 self.__HO["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
295 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"] )
296 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
297 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
298 elif asMatrix is not None:
299 matrice = numpy.matrix( asMatrix, numpy.float )
300 self.__HO["Direct"] = Operator( fromMatrix = matrice )
301 self.__HO["Tangent"] = Operator( fromMatrix = matrice )
302 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T )
305 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
307 if appliedToX is not None:
308 self.__HO["AppliedToX"] = {}
309 if type(appliedToX) is not dict:
310 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
311 for key in appliedToX.keys():
312 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
313 # Pour le cas où l'on a une vraie matrice
314 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
315 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
316 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
317 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
319 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
321 self.__HO["AppliedToX"] = None
324 self.__StoredInputs["ObservationOperator"] = self.__HO
327 # -----------------------------------------------------------
328 def setEvolutionModel(self,
335 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
336 définitions et leur sens sont les suivants :
337 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
338 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
339 "Direct" n'est pas définie, on prend la fonction "Tangent".
340 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
341 des opérateurs tangents et adjoints. On utilise par défaut des
342 différences finies non centrées ou centrées (si "withCenteredDF" est
343 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
344 du point courant ou sur le point fixe "withdX".
345 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
346 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
347 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
348 matrice fournie doit être sous une forme compatible avec le
349 constructeur de numpy.matrix.
350 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
351 être rendue disponible au même titre que les variables de calcul
352 L'argument "asFunction" peut prendre la forme complète suivante, avec
353 les valeurs par défaut standards :
354 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
355 "useApproximatedDerivatives":False,
356 "withCenteredDF" :False,
357 "withIncrement" :0.01,
359 "withAvoidingRedundancy" :True,
360 "withToleranceInRedundancy" :1.e-18,
361 "withLenghtOfRedundancy" :-1,
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 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
371 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
372 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
373 from daNumerics.ApproximatedDerivatives import FDApproximation
374 FDA = FDApproximation(
375 Function = asFunction["Direct"],
376 centeredDF = asFunction["withCenteredDF"],
377 increment = asFunction["withIncrement"],
378 dX = asFunction["withdX"],
379 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
380 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
381 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
383 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
384 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
385 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
386 elif (type(asFunction) is type({})) and \
387 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
388 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
389 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
390 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
392 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
393 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
394 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
395 elif asMatrix is not None:
396 matrice = numpy.matrix( asMatrix, numpy.float )
397 self.__EM["Direct"] = Operator( fromMatrix = matrice )
398 self.__EM["Tangent"] = Operator( fromMatrix = matrice )
399 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T )
402 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
405 self.__StoredInputs["EvolutionModel"] = self.__EM
408 def setEvolutionError(self,
410 asEyeByScalar = None,
411 asEyeByVector = None,
415 Permet de définir la covariance des erreurs de modèle :
416 - asCovariance : entrée des données, comme une matrice compatible avec
417 le constructeur de numpy.matrix
418 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
419 multiplicatif d'une matrice de corrélation identité, aucune matrice
420 n'étant donc explicitement à donner
421 - asEyeByVector : entrée des données comme un seul vecteur de variance,
422 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
423 n'étant donc explicitement à donner
424 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
425 être rendue disponible au même titre que les variables de calcul
427 self.__Q = Covariance(
428 name = "ObservationError",
429 asCovariance = asCovariance,
430 asEyeByScalar = asEyeByScalar,
431 asEyeByVector = asEyeByVector,
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 type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
528 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.matrix)
529 for member in asPersistentVector:
530 self.__U.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
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,"size"): __Xb_shape = (self.__Xb.size,)
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("The background (Xb) has no attribute of shape: problem !")
670 if self.__Y is None: __Y_shape = (0,)
671 elif hasattr(self.__Y,"size"): __Y_shape = (self.__Y.size,)
672 elif hasattr(self.__Y,"shape"):
673 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
674 else: __Y_shape = self.__Y.shape()
675 else: raise TypeError("The observation (Y) has no attribute of shape: problem !")
677 if self.__U is None: __U_shape = (0,)
678 elif hasattr(self.__U,"size"): __U_shape = (self.__U.size,)
679 elif hasattr(self.__U,"shape"):
680 if type(self.__U.shape) is tuple: __U_shape = self.__U.shape
681 else: __U_shape = self.__U.shape()
682 else: raise TypeError("The control (U) has no attribute of shape: problem !")
684 if self.__B is None: __B_shape = (0,0)
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("The a priori errors covariance matrix (B) has no attribute of shape: problem !")
690 if self.__R is None: __R_shape = (0,0)
691 elif hasattr(self.__R,"shape"):
692 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
693 else: __R_shape = self.__R.shape()
694 else: raise TypeError("The observation errors covariance matrix (R) has no attribute of shape: problem !")
696 if self.__Q is None: __Q_shape = (0,0)
697 elif hasattr(self.__Q,"shape"):
698 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
699 else: __Q_shape = self.__Q.shape()
700 else: raise TypeError("The evolution errors covariance matrix (Q) has no attribute of shape: problem !")
702 if len(self.__HO) == 0: __HO_shape = (0,0)
703 elif type(self.__HO) is type({}): __HO_shape = (0,0)
704 elif hasattr(self.__HO["Direct"],"shape"):
705 if type(self.__HO["Direct"].shape) is tuple: __HO_shape = self.__HO["Direct"].shape
706 else: __HO_shape = self.__HO["Direct"].shape()
707 else: raise TypeError("The observation operator (H) has no attribute of shape: problem !")
709 if len(self.__EM) == 0: __EM_shape = (0,0)
710 elif type(self.__EM) is type({}): __EM_shape = (0,0)
711 elif hasattr(self.__EM["Direct"],"shape"):
712 if type(self.__EM["Direct"].shape) is tuple: __EM_shape = self.__EM["Direct"].shape
713 else: __EM_shape = self.__EM["Direct"].shape()
714 else: raise TypeError("The evolution model (EM) has no attribute of shape: problem !")
716 if len(self.__CM) == 0: __CM_shape = (0,0)
717 elif type(self.__CM) is type({}): __CM_shape = (0,0)
718 elif hasattr(self.__CM["Direct"],"shape"):
719 if type(self.__CM["Direct"].shape) is tuple: __CM_shape = self.__CM["Direct"].shape
720 else: __CM_shape = self.__CM["Direct"].shape()
721 else: raise TypeError("The control model (CM) has no attribute of shape: problem !")
723 # Vérification des conditions
724 # ---------------------------
725 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
726 raise ValueError("Shape characteristic of background (Xb) is incorrect: \"%s\"."%(__Xb_shape,))
727 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
728 raise ValueError("Shape characteristic of observation (Y) is incorrect: \"%s\"."%(__Y_shape,))
730 if not( min(__B_shape) == max(__B_shape) ):
731 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) is incorrect: \"%s\"."%(__B_shape,))
732 if not( min(__R_shape) == max(__R_shape) ):
733 raise ValueError("Shape characteristic of observation errors covariance matrix (R) is incorrect: \"%s\"."%(__R_shape,))
734 if not( min(__Q_shape) == max(__Q_shape) ):
735 raise ValueError("Shape characteristic of evolution errors covariance matrix (Q) is incorrect: \"%s\"."%(__Q_shape,))
736 if not( min(__EM_shape) == max(__EM_shape) ):
737 raise ValueError("Shape characteristic of evolution operator (EM) is incorrect: \"%s\"."%(__EM_shape,))
739 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
740 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and state (X) \"%s\" are incompatible."%(__HO_shape,__Xb_shape))
741 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
742 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation (Y) \"%s\" are incompatible."%(__HO_shape,__Y_shape))
743 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] ):
744 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and a priori errors covariance matrix (B) \"%s\" are incompatible."%(__HO_shape,__B_shape))
745 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] ):
746 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation errors covariance matrix (R) \"%s\" are incompatible."%(__HO_shape,__R_shape))
748 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
749 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
750 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
751 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
752 for member in asPersistentVector:
753 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
754 __Xb_shape = min(__B_shape)
756 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) \"%s\" and background (Xb) \"%s\" are incompatible."%(__B_shape,__Xb_shape))
758 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
759 raise ValueError("Shape characteristic of observation errors covariance matrix (R) \"%s\" and observation (Y) \"%s\" are incompatible."%(__R_shape,__Y_shape))
761 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) ):
762 raise ValueError("Shape characteristic of evolution model (EM) \"%s\" and state (X) \"%s\" are incompatible."%(__EM_shape,__Xb_shape))
764 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) ):
765 raise ValueError("Shape characteristic of control model (CM) \"%s\" and control (U) \"%s\" are incompatible."%(__CM_shape,__U_shape))
767 if self.__StoredInputs.has_key("AlgorithmParameters") \
768 and self.__StoredInputs["AlgorithmParameters"].has_key("Bounds") \
769 and (type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type([]) or type(self._parameters["Bounds"]) is type(())) \
770 and (len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) != max(__Xb_shape)):
771 raise ValueError("The number \"%s\" of bound pairs for the state (X) components is different of the size \"%s\" of the state itself."%(len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]),max(__Xb_shape)))
775 # -----------------------------------------------------------
778 Permet de lancer le calcul d'assimilation.
780 Le nom de la méthode à activer est toujours "run". Les paramètres en
781 arguments de la méthode sont fixés. En sortie, on obtient les résultats
782 dans la variable de type dictionnaire "StoredVariables", qui contient en
783 particulier des objets de Persistence pour les analyses, OMA...
785 self.shape_validate()
787 self.__algorithm.run(
797 Parameters = self.__Parameters,
801 # -----------------------------------------------------------
802 def get(self, key=None):
804 Renvoie les résultats disponibles après l'exécution de la méthode
805 d'assimilation, ou les diagnostics disponibles. Attention, quand un
806 diagnostic porte le même nom qu'une variable stockée, c'est la variable
807 stockée qui est renvoyée, et le diagnostic est inatteignable.
810 if self.__algorithm.has_key(key):
811 return self.__algorithm.get( key )
812 elif self.__StoredInputs.has_key(key):
813 return self.__StoredInputs[key]
814 elif self.__StoredDiagnostics.has_key(key):
815 return self.__StoredDiagnostics[key]
817 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
819 allvariables = self.__algorithm.get()
820 allvariables.update( self.__StoredDiagnostics )
821 allvariables.update( self.__StoredInputs )
824 def get_available_variables(self):
826 Renvoie les variables potentiellement utilisables pour l'étude,
827 initialement stockées comme données d'entrées ou dans les algorithmes,
828 identifiés par les chaînes de caractères. L'algorithme doit avoir été
829 préalablement choisi sinon la méthode renvoie "None".
831 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
835 if len( self.__algorithm.keys()) > 0:
836 variables.extend( self.__algorithm.get().keys() )
837 if len( self.__StoredInputs.keys() ) > 0:
838 variables.extend( self.__StoredInputs.keys() )
842 def get_available_algorithms(self):
844 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
845 par les chaînes de caractères.
848 for directory in sys.path:
849 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
850 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
851 root, ext = os.path.splitext(fname)
852 if ext == '.py' and root != '__init__':
857 def get_available_diagnostics(self):
859 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
860 par les chaînes de caractères.
863 for directory in sys.path:
864 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
865 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
866 root, ext = os.path.splitext(fname)
867 if ext == '.py' and root != '__init__':
872 # -----------------------------------------------------------
873 def get_algorithms_main_path(self):
875 Renvoie le chemin pour le répertoire principal contenant les algorithmes
876 dans un sous-répertoire "daAlgorithms"
880 def add_algorithms_path(self, asPath=None):
882 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
883 se trouve un sous-répertoire "daAlgorithms"
885 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
886 pas indispensable de le rajouter ici.
888 if not os.path.isdir(asPath):
889 raise ValueError("The given "+asPath+" argument must exist as a directory")
890 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
891 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
892 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
893 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
894 sys.path.insert(0, os.path.abspath(asPath))
895 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
898 def get_diagnostics_main_path(self):
900 Renvoie le chemin pour le répertoire principal contenant les diagnostics
901 dans un sous-répertoire "daDiagnostics"
905 def add_diagnostics_path(self, asPath=None):
907 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
908 se trouve un sous-répertoire "daDiagnostics"
910 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
911 pas indispensable de le rajouter ici.
913 if not os.path.isdir(asPath):
914 raise ValueError("The given "+asPath+" argument must exist as a directory")
915 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
916 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
917 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
918 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
919 sys.path.insert(0, os.path.abspath(asPath))
920 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
923 # -----------------------------------------------------------
924 def setDataObserver(self,
927 HookParameters = None,
931 Permet d'associer un observer à une ou des variables nommées gérées en
932 interne, activable selon des règles définies dans le Scheduler. A chaque
933 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
934 les arguments (variable persistante VariableName, paramètres HookParameters).
937 if type( self.__algorithm ) is dict:
938 raise ValueError("No observer can be build before choosing an algorithm.")
940 # Vérification du nom de variable et typage
941 # -----------------------------------------
942 if type( VariableName ) is str:
943 VariableNames = [VariableName,]
944 elif type( VariableName ) is list:
945 VariableNames = map( str, VariableName )
947 raise ValueError("The observer requires a name or a list of names of variables.")
949 # Association interne de l'observer à la variable
950 # -----------------------------------------------
951 for n in VariableNames:
952 if not self.__algorithm.has_key( n ):
953 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
955 self.__algorithm.StoredVariables[ n ].setDataObserver(
956 Scheduler = Scheduler,
957 HookFunction = HookFunction,
958 HookParameters = HookParameters,
961 def removeDataObserver(self,
966 Permet de retirer un observer à une ou des variable nommée.
969 if type( self.__algorithm ) is dict:
970 raise ValueError("No observer can be removed before choosing an algorithm.")
972 # Vérification du nom de variable et typage
973 # -----------------------------------------
974 if type( VariableName ) is str:
975 VariableNames = [VariableName,]
976 elif type( VariableName ) is list:
977 VariableNames = map( str, VariableName )
979 raise ValueError("The observer requires a name or a list of names of variables.")
981 # Association interne de l'observer à la variable
982 # -----------------------------------------------
983 for n in VariableNames:
984 if not self.__algorithm.has_key( n ):
985 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
987 self.__algorithm.StoredVariables[ n ].removeDataObserver(
988 HookFunction = HookFunction,
991 # -----------------------------------------------------------
992 def setDebug(self, level=10):
994 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
995 appel pour changer le niveau de verbosité, avec :
996 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
999 log = logging.getLogger()
1000 log.setLevel( level )
1002 def unsetDebug(self):
1004 Remet le logger au niveau par défaut
1007 log = logging.getLogger()
1008 log.setLevel( logging.WARNING )
1010 def prepare_to_pickle(self):
1011 self.__algorithmFile = None
1012 self.__diagnosticFile = None
1017 # ==============================================================================
1018 if __name__ == "__main__":
1019 print '\n AUTODIAGNOSTIC \n'