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) d'observations
55 - State : vecteur d'état dont une partie est le vecteur de contrôle.
56 Cette information n'est utile que si l'on veut faire des calculs sur
57 l'état complet, mais elle n'est pas indispensable pour l'assimilation.
58 - Control : vecteur X contenant toutes les variables de contrôle, i.e.
59 les paramètres ou l'état dont on veut estimer la valeur pour obtenir
61 - ObservationOperator...: opérateur d'observation H
63 Les observations présentent une erreur dont la matrice de covariance est
64 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
67 self.__name = str(name)
78 self.__X = Persistence.OneVector()
79 self.__Parameters = {}
80 self.__StoredDiagnostics = {}
81 self.__StoredInputs = {}
83 self.__B_scalar = None
84 self.__R_scalar = None
85 self.__Q_scalar = None
86 self.__Parameters["B_scalar"] = None
87 self.__Parameters["R_scalar"] = None
88 self.__Parameters["Q_scalar"] = None
90 # Variables temporaires
92 self.__algorithmFile = None
93 self.__algorithmName = None
94 self.__diagnosticFile = None
96 # Récupère le chemin du répertoire parent et l'ajoute au path
97 # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
98 # qui est activée dans Persistence)
99 self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
100 sys.path.insert(0, self.__parent)
101 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
103 # ---------------------------------------------------------
104 def setBackground(self,
106 asPersistentVector = None,
111 Permet de définir l'estimation a priori :
112 - asVector : entrée des données, comme un vecteur compatible avec le
113 constructeur de numpy.matrix
114 - asPersistentVector : entrée des données, comme un vecteur de type
115 persistent contruit avec la classe ad-hoc "Persistence"
116 - Scheduler est le contrôle temporel des données
117 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
118 être rendue disponible au même titre que les variables de calcul
120 if asVector is not None:
121 if type( asVector ) is type( numpy.matrix([]) ):
122 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
124 self.__Xb = numpy.matrix( asVector, numpy.float ).T
125 elif asPersistentVector is not None:
126 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
127 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
128 for member in asPersistentVector:
129 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
131 self.__Xb = asPersistentVector
133 raise ValueError("Error: improperly defined background, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
135 self.__StoredInputs["Background"] = self.__Xb
138 def setBackgroundError(self,
140 asEyeByScalar = None,
141 asEyeByVector = None,
145 Permet de définir la covariance des erreurs d'ébauche :
146 - asCovariance : entrée des données, comme une matrice compatible avec
147 le constructeur de numpy.matrix
148 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
149 multiplicatif d'une matrice de corrélation identité, aucune matrice
150 n'étant donc explicitement à donner
151 - asEyeByVector : entrée des données comme un seul vecteur de variance,
152 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
153 n'étant donc explicitement à donner
154 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
155 être rendue disponible au même titre que les variables de calcul
157 if asEyeByScalar is not None:
158 self.__B_scalar = float(asEyeByScalar)
160 elif asEyeByVector is not None:
161 self.__B_scalar = None
162 self.__B = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
163 elif asCovariance is not None:
164 self.__B_scalar = None
165 self.__B = numpy.matrix( asCovariance, numpy.float )
167 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")
169 self.__Parameters["B_scalar"] = self.__B_scalar
171 self.__StoredInputs["BackgroundError"] = self.__B
174 # -----------------------------------------------------------
175 def setObservation(self,
177 asPersistentVector = None,
182 Permet de définir les observations :
183 - asVector : entrée des données, comme un vecteur compatible avec le
184 constructeur de numpy.matrix
185 - asPersistentVector : entrée des données, comme un vecteur de type
186 persistent contruit avec la classe ad-hoc "Persistence"
187 - Scheduler est le contrôle temporel des données disponibles
188 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
189 être rendue disponible au même titre que les variables de calcul
191 if asVector is not None:
192 if type( asVector ) is type( numpy.matrix([]) ):
193 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
195 self.__Y = numpy.matrix( asVector, numpy.float ).T
196 elif asPersistentVector is not None:
197 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
198 self.__Y = Persistence.OneVector("Observation", basetype=numpy.matrix)
199 for member in asPersistentVector:
200 self.__Y.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
202 self.__Y = asPersistentVector
204 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
206 self.__StoredInputs["Observation"] = self.__Y
209 def setObservationError(self,
211 asEyeByScalar = None,
212 asEyeByVector = None,
216 Permet de définir la covariance des erreurs d'observations :
217 - asCovariance : entrée des données, comme une matrice compatible avec
218 le constructeur de numpy.matrix
219 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
220 multiplicatif d'une matrice de corrélation identité, aucune matrice
221 n'étant donc explicitement à donner
222 - asEyeByVector : entrée des données comme un seul vecteur de variance,
223 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
224 n'étant donc explicitement à donner
225 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
226 être rendue disponible au même titre que les variables de calcul
228 if asEyeByScalar is not None:
229 self.__R_scalar = float(asEyeByScalar)
231 elif asEyeByVector is not None:
232 self.__R_scalar = None
233 self.__R = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
234 elif asCovariance is not None:
235 self.__R_scalar = None
236 self.__R = numpy.matrix( asCovariance, numpy.float )
238 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")
240 self.__Parameters["R_scalar"] = self.__R_scalar
242 self.__StoredInputs["ObservationError"] = self.__R
245 def setObservationOperator(self,
252 Permet de définir un opérateur d'observation H. L'ordre de priorité des
253 définitions et leur sens sont les suivants :
254 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
255 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
256 "Direct" n'est pas définie, on prend la fonction "Tangent".
257 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
258 des opérateurs tangents et adjoints. On utilise par défaut des
259 différences finies non centrées ou centrées (si "withCenteredDF" est
260 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
261 du point courant ou sur le point fixe "withdX".
262 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
263 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
264 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
265 matrice fournie doit être sous une forme compatible avec le
266 constructeur de numpy.matrix.
267 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
268 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
269 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
270 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
271 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
272 être rendue disponible au même titre que les variables de calcul
273 L'argument "asFunction" peut prendre la forme complète suivante, avec
274 les valeurs par défaut standards :
275 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
276 "useApproximatedDerivatives":False,
277 "withCenteredDF" :False,
278 "withIncrement" :0.01,
280 "withAvoidingRedundancy" :True,
281 "withToleranceInRedundancy" :1.e-18,
282 "withLenghtOfRedundancy" :-1,
285 if (type(asFunction) is type({})) and \
286 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
287 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
288 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
289 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
290 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
291 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
292 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
293 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
294 from daNumerics.ApproximatedDerivatives import FDApproximation
295 FDA = FDApproximation(
296 Function = asFunction["Direct"],
297 centeredDF = asFunction["withCenteredDF"],
298 increment = asFunction["withIncrement"],
299 dX = asFunction["withdX"],
300 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
301 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
302 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
304 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator )
305 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
306 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
307 elif (type(asFunction) is type({})) and \
308 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
309 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
310 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
311 self.__HO["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
313 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"] )
314 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
315 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
316 elif asMatrix is not None:
317 matrice = numpy.matrix( asMatrix, numpy.float )
318 self.__HO["Direct"] = Operator( fromMatrix = matrice )
319 self.__HO["Tangent"] = Operator( fromMatrix = matrice )
320 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T )
323 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
325 if appliedToX is not None:
326 self.__HO["AppliedToX"] = {}
327 if type(appliedToX) is not dict:
328 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
329 for key in appliedToX.keys():
330 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
331 # Pour le cas où l'on a une vraie matrice
332 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
333 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
334 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
335 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
337 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
339 self.__HO["AppliedToX"] = None
342 self.__StoredInputs["ObservationOperator"] = self.__HO
345 # -----------------------------------------------------------
346 def setEvolutionModel(self,
353 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
354 définitions et leur sens sont les suivants :
355 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
356 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
357 "Direct" n'est pas définie, on prend la fonction "Tangent".
358 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
359 des opérateurs tangents et adjoints. On utilise par défaut des
360 différences finies non centrées ou centrées (si "withCenteredDF" est
361 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
362 du point courant ou sur le point fixe "withdX".
363 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
364 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
365 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
366 matrice fournie doit être sous une forme compatible avec le
367 constructeur de numpy.matrix.
368 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
369 être rendue disponible au même titre que les variables de calcul
370 L'argument "asFunction" peut prendre la forme complète suivante, avec
371 les valeurs par défaut standards :
372 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
373 "useApproximatedDerivatives":False,
374 "withCenteredDF" :False,
375 "withIncrement" :0.01,
377 "withAvoidingRedundancy" :True,
378 "withToleranceInRedundancy" :1.e-18,
379 "withLenghtOfRedundancy" :-1,
382 if (type(asFunction) is type({})) and \
383 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
384 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
385 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
386 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
387 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
388 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
389 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
390 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
391 from daNumerics.ApproximatedDerivatives import FDApproximation
392 FDA = FDApproximation(
393 Function = asFunction["Direct"],
394 centeredDF = asFunction["withCenteredDF"],
395 increment = asFunction["withIncrement"],
396 dX = asFunction["withdX"],
397 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
398 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
399 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
401 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
402 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
403 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
404 elif (type(asFunction) is type({})) and \
405 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
406 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
407 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
408 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
410 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
411 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
412 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
413 elif asMatrix is not None:
414 matrice = numpy.matrix( asMatrix, numpy.float )
415 self.__EM["Direct"] = Operator( fromMatrix = matrice )
416 self.__EM["Tangent"] = Operator( fromMatrix = matrice )
417 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T )
420 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
423 self.__StoredInputs["EvolutionModel"] = self.__EM
426 def setEvolutionError(self,
428 asEyeByScalar = None,
429 asEyeByVector = None,
433 Permet de définir la covariance des erreurs de modèle :
434 - asCovariance : entrée des données, comme une matrice compatible avec
435 le constructeur de numpy.matrix
436 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
437 multiplicatif d'une matrice de corrélation identité, aucune matrice
438 n'étant donc explicitement à donner
439 - asEyeByVector : entrée des données comme un seul vecteur de variance,
440 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
441 n'étant donc explicitement à donner
442 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
443 être rendue disponible au même titre que les variables de calcul
445 if asEyeByScalar is not None:
446 self.__Q_scalar = float(asEyeByScalar)
448 elif asEyeByVector is not None:
449 self.__Q_scalar = None
450 self.__Q = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
451 elif asCovariance is not None:
452 self.__Q_scalar = None
453 self.__Q = numpy.matrix( asCovariance, numpy.float )
455 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")
457 self.__Parameters["Q_scalar"] = self.__Q_scalar
459 self.__StoredInputs["EvolutionError"] = self.__Q
462 # -----------------------------------------------------------
463 def setControlModel(self,
464 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
465 "useApproximatedDerivatives":False,
466 "withCenteredDF" :False,
467 "withIncrement" :0.01,
475 Permet de définir un opérateur de controle C. L'ordre de priorité des
476 définitions et leur sens sont les suivants :
477 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
478 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
479 "Direct" n'est pas définie, on prend la fonction "Tangent".
480 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
481 des opérateurs tangents et adjoints. On utilise par défaut des
482 différences finies non centrées ou centrées (si "withCenteredDF" est
483 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
484 du point courant ou sur le point fixe "withdX".
485 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
486 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
487 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
488 matrice fournie doit être sous une forme compatible avec le
489 constructeur de numpy.matrix.
490 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
491 être rendue disponible au même titre que les variables de calcul
493 if (type(asFunction) is type({})) and \
494 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
495 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
496 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
497 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
498 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
499 from daNumerics.ApproximatedDerivatives import FDApproximation
500 FDA = FDApproximation(
501 Function = asFunction["Direct"],
502 centeredDF = asFunction["withCenteredDF"],
503 increment = asFunction["withIncrement"],
504 dX = asFunction["withdX"] )
505 self.__CM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
506 self.__CM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
507 self.__CM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
508 elif (type(asFunction) is type({})) and \
509 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
510 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
511 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
512 self.__CM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
514 self.__CM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
515 self.__CM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
516 self.__CM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
517 elif asMatrix is not None:
518 matrice = numpy.matrix( asMatrix, numpy.float )
519 self.__CM["Direct"] = Operator( fromMatrix = matrice )
520 self.__CM["Tangent"] = Operator( fromMatrix = matrice )
521 self.__CM["Adjoint"] = Operator( fromMatrix = matrice.T )
524 raise ValueError("Improperly defined input control model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
527 self.__StoredInputs["ControlModel"] = self.__CM
530 def setControlInput(self,
532 asPersistentVector = None,
537 Permet de définir le controle en entree :
538 - asVector : entrée des données, comme un vecteur compatible avec le
539 constructeur de numpy.matrix
540 - asPersistentVector : entrée des données, comme un vecteur de type
541 persistent contruit avec la classe ad-hoc "Persistence"
542 - Scheduler est le contrôle temporel des données disponibles
543 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
544 être rendue disponible au même titre que les variables de calcul
546 if asVector is not None:
547 if isinstance(asVector,numpy.matrix):
548 self.__U = numpy.matrix( asVector.A1, numpy.float ).T
550 self.__U = numpy.matrix( asVector, numpy.float ).T
551 elif asPersistentVector is not None:
552 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
553 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.matrix)
554 for member in asPersistentVector:
555 self.__U.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
557 self.__U = asPersistentVector
559 raise ValueError("Error: improperly defined control input, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
561 self.__StoredInputs["ControlInput"] = self.__U
564 # -----------------------------------------------------------
565 def setControls (self,
570 Permet de définir la valeur initiale du vecteur X contenant toutes les
571 variables de contrôle, i.e. les paramètres ou l'état dont on veut
572 estimer la valeur pour obtenir les observations. C'est utile pour un
573 algorithme itératif/incrémental.
574 - asVector : entrée des données, comme un vecteur compatible avec le
575 constructeur de numpy.matrix.
576 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
577 être rendue disponible au même titre que les variables de calcul
579 if asVector is not None:
580 self.__X.store( asVector )
582 self.__StoredInputs["Controls"] = self.__X
585 # -----------------------------------------------------------
586 def setAlgorithm(self, choice = None ):
588 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
589 d'assimilation. L'argument est un champ caractère se rapportant au nom
590 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
591 d'assimilation sur les arguments fixes.
594 raise ValueError("Error: algorithm choice has to be given")
595 if self.__algorithmName is not None:
596 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
597 daDirectory = "daAlgorithms"
599 # Recherche explicitement le fichier complet
600 # ------------------------------------------
602 for directory in sys.path:
603 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
604 module_path = os.path.abspath(os.path.join(directory, daDirectory))
605 if module_path is None:
606 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
608 # Importe le fichier complet comme un module
609 # ------------------------------------------
611 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
612 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
613 self.__algorithmName = str(choice)
614 sys.path = sys_path_tmp ; del sys_path_tmp
615 except ImportError, e:
616 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
618 # Instancie un objet du type élémentaire du fichier
619 # -------------------------------------------------
620 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
621 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
624 def setAlgorithmParameters(self, asDico=None):
626 Permet de définir les paramètres de l'algorithme, sous la forme d'un
629 if asDico is not None:
630 self.__Parameters.update( dict( asDico ) )
632 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
635 def getAlgorithmParameters(self, noDetails=True):
637 Renvoie la liste des paramètres requis selon l'algorithme
639 return self.__algorithm.getRequiredParameters(noDetails)
641 # -----------------------------------------------------------
642 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
644 Permet de sélectionner un diagnostic a effectuer.
647 raise ValueError("Error: diagnostic choice has to be given")
648 daDirectory = "daDiagnostics"
650 # Recherche explicitement le fichier complet
651 # ------------------------------------------
653 for directory in sys.path:
654 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
655 module_path = os.path.abspath(os.path.join(directory, daDirectory))
656 if module_path is None:
657 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
659 # Importe le fichier complet comme un module
660 # ------------------------------------------
662 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
663 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
664 sys.path = sys_path_tmp ; del sys_path_tmp
665 except ImportError, e:
666 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
668 # Instancie un objet du type élémentaire du fichier
669 # -------------------------------------------------
670 if self.__StoredInputs.has_key(name):
671 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
672 elif self.__StoredDiagnostics.has_key(name):
673 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
675 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
679 parameters = parameters )
682 # -----------------------------------------------------------
683 def shape_validate(self):
685 Validation de la correspondance correcte des tailles des variables et
686 des matrices s'il y en a.
688 if self.__Xb is None: __Xb_shape = (0,)
689 elif hasattr(self.__Xb,"shape"):
690 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
691 else: __Xb_shape = self.__Xb.shape()
692 else: raise TypeError("Xb has no attribute of shape: problem !")
694 if self.__Y is None: __Y_shape = (0,)
695 elif hasattr(self.__Y,"shape"):
696 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
697 else: __Y_shape = self.__Y.shape()
698 else: raise TypeError("Y has no attribute of shape: problem !")
700 if self.__U is None: __U_shape = (0,)
701 elif hasattr(self.__U,"shape"):
702 if type(self.__U.shape) is tuple: __U_shape = self.__U.shape
703 else: __U_shape = self.__U.shape()
704 else: raise TypeError("U has no attribute of shape: problem !")
706 if self.__B is None and self.__B_scalar is None:
708 elif self.__B is None and self.__B_scalar is not None:
709 __B_shape = (max(__Xb_shape),max(__Xb_shape))
710 elif hasattr(self.__B,"shape"):
711 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
712 else: __B_shape = self.__B.shape()
713 else: raise TypeError("B has no attribute of shape: problem !")
715 if self.__R is None and self.__R_scalar is None:
717 elif self.__R is None and self.__R_scalar is not None:
718 __R_shape = (max(__Y_shape),max(__Y_shape))
719 elif hasattr(self.__R,"shape"):
720 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
721 else: __R_shape = self.__R.shape()
722 else: raise TypeError("R has no attribute of shape: problem !")
724 if self.__Q is None: __Q_shape = (0,0)
725 elif hasattr(self.__Q,"shape"):
726 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
727 else: __Q_shape = self.__Q.shape()
728 else: raise TypeError("Q has no attribute of shape: problem !")
730 if len(self.__HO) == 0: __HO_shape = (0,0)
731 elif type(self.__HO) is type({}): __HO_shape = (0,0)
732 elif hasattr(self.__HO["Direct"],"shape"):
733 if type(self.__HO["Direct"].shape) is tuple: __HO_shape = self.__HO["Direct"].shape
734 else: __HO_shape = self.__HO["Direct"].shape()
735 else: raise TypeError("H has no attribute of shape: problem !")
737 if len(self.__EM) == 0: __EM_shape = (0,0)
738 elif type(self.__EM) is type({}): __EM_shape = (0,0)
739 elif hasattr(self.__EM["Direct"],"shape"):
740 if type(self.__EM["Direct"].shape) is tuple: __EM_shape = self.__EM["Direct"].shape
741 else: __EM_shape = self.__EM["Direct"].shape()
742 else: raise TypeError("EM has no attribute of shape: problem !")
744 if len(self.__CM) == 0: __CM_shape = (0,0)
745 elif type(self.__CM) is type({}): __CM_shape = (0,0)
746 elif hasattr(self.__CM["Direct"],"shape"):
747 if type(self.__CM["Direct"].shape) is tuple: __CM_shape = self.__CM["Direct"].shape
748 else: __CM_shape = self.__CM["Direct"].shape()
749 else: raise TypeError("CM has no attribute of shape: problem !")
751 # Vérification des conditions
752 # ---------------------------
753 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
754 raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
755 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
756 raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
758 if not( min(__B_shape) == max(__B_shape) ):
759 raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
760 if not( min(__R_shape) == max(__R_shape) ):
761 raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
762 if not( min(__Q_shape) == max(__Q_shape) ):
763 raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
764 if not( min(__EM_shape) == max(__EM_shape) ):
765 raise ValueError("Shape characteristic of EM is incorrect: \"%s\""%(__EM_shape,))
767 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
768 raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__HO_shape,__Xb_shape))
769 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
770 raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__HO_shape,__Y_shape))
771 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] ):
772 raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__HO_shape,__B_shape))
773 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] ):
774 raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__HO_shape,__R_shape))
776 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
777 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
778 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
779 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
780 for member in asPersistentVector:
781 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
782 __Xb_shape = min(__B_shape)
784 raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
786 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
787 raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
789 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) ):
790 raise ValueError("Shape characteristic of EM \"%s\" and X \"%s\" are incompatible"%(__EM_shape,__Xb_shape))
792 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) ):
793 raise ValueError("Shape characteristic of CM \"%s\" and U \"%s\" are incompatible"%(__CM_shape,__U_shape))
797 # -----------------------------------------------------------
800 Permet de lancer le calcul d'assimilation.
802 Le nom de la méthode à activer est toujours "run". Les paramètres en
803 arguments de la méthode sont fixés. En sortie, on obtient les résultats
804 dans la variable de type dictionnaire "StoredVariables", qui contient en
805 particulier des objets de Persistence pour les analyses, OMA...
807 self.shape_validate()
809 self.__algorithm.run(
819 Parameters = self.__Parameters,
823 # -----------------------------------------------------------
824 def get(self, key=None):
826 Renvoie les résultats disponibles après l'exécution de la méthode
827 d'assimilation, ou les diagnostics disponibles. Attention, quand un
828 diagnostic porte le même nom qu'une variable stockée, c'est la variable
829 stockée qui est renvoyée, et le diagnostic est inatteignable.
832 if self.__algorithm.has_key(key):
833 return self.__algorithm.get( key )
834 elif self.__StoredInputs.has_key(key):
835 return self.__StoredInputs[key]
836 elif self.__StoredDiagnostics.has_key(key):
837 return self.__StoredDiagnostics[key]
839 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
841 allvariables = self.__algorithm.get()
842 allvariables.update( self.__StoredDiagnostics )
843 allvariables.update( self.__StoredInputs )
846 def get_available_variables(self):
848 Renvoie les variables potentiellement utilisables pour l'étude,
849 initialement stockées comme données d'entrées ou dans les algorithmes,
850 identifiés par les chaînes de caractères. L'algorithme doit avoir été
851 préalablement choisi sinon la méthode renvoie "None".
853 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
857 if len( self.__algorithm.keys()) > 0:
858 variables.extend( self.__algorithm.get().keys() )
859 if len( self.__StoredInputs.keys() ) > 0:
860 variables.extend( self.__StoredInputs.keys() )
864 def get_available_algorithms(self):
866 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
867 par les chaînes de caractères.
870 for directory in sys.path:
871 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
872 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
873 root, ext = os.path.splitext(fname)
874 if ext == '.py' and root != '__init__':
879 def get_available_diagnostics(self):
881 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
882 par les chaînes de caractères.
885 for directory in sys.path:
886 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
887 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
888 root, ext = os.path.splitext(fname)
889 if ext == '.py' and root != '__init__':
894 # -----------------------------------------------------------
895 def get_algorithms_main_path(self):
897 Renvoie le chemin pour le répertoire principal contenant les algorithmes
898 dans un sous-répertoire "daAlgorithms"
902 def add_algorithms_path(self, asPath=None):
904 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
905 se trouve un sous-répertoire "daAlgorithms"
907 Remarque : si le chemin a déjà été ajouté pour les diagnostics, 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,"daAlgorithms")):
913 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
914 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
915 raise ValueError("The given \""+asPath+"/daAlgorithms\" 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 def get_diagnostics_main_path(self):
922 Renvoie le chemin pour le répertoire principal contenant les diagnostics
923 dans un sous-répertoire "daDiagnostics"
927 def add_diagnostics_path(self, asPath=None):
929 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
930 se trouve un sous-répertoire "daDiagnostics"
932 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
933 pas indispensable de le rajouter ici.
935 if not os.path.isdir(asPath):
936 raise ValueError("The given "+asPath+" argument must exist as a directory")
937 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
938 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
939 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
940 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
941 sys.path.insert(0, os.path.abspath(asPath))
942 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
945 # -----------------------------------------------------------
946 def setDataObserver(self,
949 HookParameters = None,
953 Permet d'associer un observer à une ou des variables nommées gérées en
954 interne, activable selon des règles définies dans le Scheduler. A chaque
955 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
956 les arguments (variable persistante VariableName, paramètres HookParameters).
959 if type( self.__algorithm ) is dict:
960 raise ValueError("No observer can be build before choosing an algorithm.")
962 # Vérification du nom de variable et typage
963 # -----------------------------------------
964 if type( VariableName ) is str:
965 VariableNames = [VariableName,]
966 elif type( VariableName ) is list:
967 VariableNames = map( str, VariableName )
969 raise ValueError("The observer requires a name or a list of names of variables.")
971 # Association interne de l'observer à la variable
972 # -----------------------------------------------
973 for n in VariableNames:
974 if not self.__algorithm.has_key( n ):
975 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
977 self.__algorithm.StoredVariables[ n ].setDataObserver(
978 Scheduler = Scheduler,
979 HookFunction = HookFunction,
980 HookParameters = HookParameters,
983 def removeDataObserver(self,
988 Permet de retirer un observer à une ou des variable nommée.
991 if type( self.__algorithm ) is dict:
992 raise ValueError("No observer can be removed before choosing an algorithm.")
994 # Vérification du nom de variable et typage
995 # -----------------------------------------
996 if type( VariableName ) is str:
997 VariableNames = [VariableName,]
998 elif type( VariableName ) is list:
999 VariableNames = map( str, VariableName )
1001 raise ValueError("The observer requires a name or a list of names of variables.")
1003 # Association interne de l'observer à la variable
1004 # -----------------------------------------------
1005 for n in VariableNames:
1006 if not self.__algorithm.has_key( n ):
1007 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
1009 self.__algorithm.StoredVariables[ n ].removeDataObserver(
1010 HookFunction = HookFunction,
1013 # -----------------------------------------------------------
1014 def setDebug(self, level=10):
1016 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
1017 appel pour changer le niveau de verbosité, avec :
1018 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
1021 log = logging.getLogger()
1022 log.setLevel( level )
1024 def unsetDebug(self):
1026 Remet le logger au niveau par défaut
1029 log = logging.getLogger()
1030 log.setLevel( logging.WARNING )
1032 def prepare_to_pickle(self):
1033 self.__algorithmFile = None
1034 self.__diagnosticFile = None
1039 # ==============================================================================
1040 if __name__ == "__main__":
1041 print '\n AUTODIAGNOSTIC \n'
1043 ADD = AssimilationStudy("Ma premiere etude BLUE")
1045 ADD.setBackground (asVector = [0, 1, 2])
1046 ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
1047 ADD.setObservation (asVector = [0.5, 1.5, 2.5])
1048 ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
1049 ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
1051 ADD.setAlgorithm(choice="Blue")
1055 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1056 print "Ebauche :", [0, 1, 2]
1057 print "Observation :", [0.5, 1.5, 2.5]
1058 print "Demi-somme :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
1059 print " qui doit être identique à :"
1060 print "Analyse résultante :", ADD.get("Analysis")[0]
1063 print "Algorithmes disponibles.......................:", ADD.get_available_algorithms()
1064 # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path()
1065 print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics()
1066 # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path()
1067 print "Variables disponibles.........................:", ADD.get_available_variables()
1070 print "Paramètres requis par algorithme :"
1071 for algo in ADD.get_available_algorithms():
1072 tmpADD = AssimilationStudy("Un algorithme")
1073 tmpADD.setAlgorithm(choice=algo)
1074 print " %25s : %s"%(algo,tmpADD.getAlgorithmParameters())
1078 ADD.setDiagnostic("RMS", "Ma RMS")
1080 liste = ADD.get().keys()
1082 print "Variables et diagnostics nommés disponibles...:", liste
1085 print "Exemple de mise en place d'un observeur :"
1086 def obs(var=None,info=None):
1087 print " ---> Mise en oeuvre de l'observer"
1088 print " var =",var[-1]
1089 print " info =",info
1090 ADD.setDataObserver( 'Analysis', HookFunction=obs, Scheduler = [2, 4], HookParameters = "Second observer")
1091 # Attention, il faut décaler le stockage de 1 pour suivre le pas interne
1092 # car le pas 0 correspond à l'analyse ci-dessus.
1093 for i in range(1,6):
1095 print "Action sur la variable observée, étape :",i
1096 ADD.get('Analysis').store( [i, i, i] )
1099 print "Mise en debug et hors debug"
1100 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1104 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
1106 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()