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
38 logging.debug("Succeed initial import of scipy.optimize with Scipy %s", scipy.version.version)
40 logging.debug("Fail initial import of scipy.optimize")
42 from BasicObjects import Operator, Covariance
43 from PlatformInfo import uniq
45 # ==============================================================================
46 class AssimilationStudy:
48 Cette classe sert d'interface pour l'utilisation de l'assimilation de
49 données. Elle contient les méthodes ou accesseurs nécessaires à la
50 construction d'un calcul d'assimilation.
52 def __init__(self, name=""):
54 Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
55 élémentaire. Ces variables sont ensuite disponibles pour implémenter un
56 algorithme élémentaire particulier.
58 - Background : vecteur Xb
59 - Observation : vecteur Y (potentiellement temporel) d'observations
60 - State : vecteur d'état dont une partie est le vecteur de contrôle.
61 Cette information n'est utile que si l'on veut faire des calculs sur
62 l'état complet, mais elle n'est pas indispensable pour l'assimilation.
63 - Control : vecteur X contenant toutes les variables de contrôle, i.e.
64 les paramètres ou l'état dont on veut estimer la valeur pour obtenir
66 - ObservationOperator...: opérateur d'observation H
68 Les observations présentent une erreur dont la matrice de covariance est
69 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
72 self.__name = str(name)
83 self.__X = Persistence.OneVector()
84 self.__Parameters = {}
85 self.__StoredDiagnostics = {}
86 self.__StoredInputs = {}
88 # Variables temporaires
90 self.__algorithmFile = None
91 self.__algorithmName = None
92 self.__diagnosticFile = None
94 # Récupère le chemin du répertoire parent et l'ajoute au path
95 # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
96 # qui est activée dans Persistence)
97 self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
98 sys.path.insert(0, self.__parent)
99 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
101 # ---------------------------------------------------------
102 def setBackground(self,
104 asPersistentVector = None,
109 Permet de définir l'estimation a priori :
110 - asVector : entrée des données, comme un vecteur compatible avec le
111 constructeur de numpy.matrix
112 - asPersistentVector : entrée des données, comme un vecteur de type
113 persistent contruit avec la classe ad-hoc "Persistence"
114 - Scheduler est le contrôle temporel des données
115 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
116 être rendue disponible au même titre que les variables de calcul
118 if asVector is not None:
119 if type( asVector ) is type( numpy.matrix([]) ):
120 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
122 self.__Xb = numpy.matrix( asVector, numpy.float ).T
123 elif asPersistentVector is not None:
124 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
125 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
126 for member in asPersistentVector:
127 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
129 self.__Xb = asPersistentVector
131 raise ValueError("Error: improperly defined background, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
133 self.__StoredInputs["Background"] = self.__Xb
136 def setBackgroundError(self,
138 asEyeByScalar = None,
139 asEyeByVector = None,
143 Permet de définir la covariance des erreurs d'ébauche :
144 - asCovariance : entrée des données, comme une matrice compatible avec
145 le constructeur de numpy.matrix
146 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
147 multiplicatif d'une matrice de corrélation identité, aucune matrice
148 n'étant donc explicitement à donner
149 - asEyeByVector : entrée des données comme un seul vecteur de variance,
150 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
151 n'étant donc explicitement à donner
152 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
153 être rendue disponible au même titre que les variables de calcul
155 self.__B = Covariance(
156 name = "BackgroundError",
157 asCovariance = asCovariance,
158 asEyeByScalar = asEyeByScalar,
159 asEyeByVector = asEyeByVector,
162 self.__StoredInputs["BackgroundError"] = self.__B
165 # -----------------------------------------------------------
166 def setObservation(self,
168 asPersistentVector = None,
173 Permet de définir les observations :
174 - asVector : entrée des données, comme un vecteur compatible avec le
175 constructeur de numpy.matrix
176 - asPersistentVector : entrée des données, comme un vecteur de type
177 persistent contruit avec la classe ad-hoc "Persistence"
178 - Scheduler est le contrôle temporel des données disponibles
179 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
180 être rendue disponible au même titre que les variables de calcul
182 if asVector is not None:
183 if type( asVector ) is type( numpy.matrix([]) ):
184 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
186 self.__Y = numpy.matrix( asVector, numpy.float ).T
187 elif asPersistentVector is not None:
188 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
189 self.__Y = Persistence.OneVector("Observation", basetype=numpy.matrix)
190 for member in asPersistentVector:
191 self.__Y.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
193 self.__Y = asPersistentVector
195 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
197 self.__StoredInputs["Observation"] = self.__Y
200 def setObservationError(self,
202 asEyeByScalar = None,
203 asEyeByVector = None,
207 Permet de définir la covariance des erreurs d'observations :
208 - asCovariance : entrée des données, comme une matrice compatible avec
209 le constructeur de numpy.matrix
210 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
211 multiplicatif d'une matrice de corrélation identité, aucune matrice
212 n'étant donc explicitement à donner
213 - asEyeByVector : entrée des données comme un seul vecteur de variance,
214 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
215 n'étant donc explicitement à donner
216 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
217 être rendue disponible au même titre que les variables de calcul
219 self.__R = Covariance(
220 name = "ObservationError",
221 asCovariance = asCovariance,
222 asEyeByScalar = asEyeByScalar,
223 asEyeByVector = asEyeByVector,
226 self.__StoredInputs["ObservationError"] = self.__R
229 def setObservationOperator(self,
236 Permet de définir un opérateur d'observation H. L'ordre de priorité des
237 définitions et leur sens sont les suivants :
238 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
239 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
240 "Direct" n'est pas définie, on prend la fonction "Tangent".
241 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
242 des opérateurs tangents et adjoints. On utilise par défaut des
243 différences finies non centrées ou centrées (si "withCenteredDF" est
244 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
245 du point courant ou sur le point fixe "withdX".
246 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
247 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
248 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
249 matrice fournie doit être sous une forme compatible avec le
250 constructeur de numpy.matrix.
251 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
252 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
253 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
254 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
255 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
256 être rendue disponible au même titre que les variables de calcul
257 L'argument "asFunction" peut prendre la forme complète suivante, avec
258 les valeurs par défaut standards :
259 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
260 "useApproximatedDerivatives":False,
261 "withCenteredDF" :False,
262 "withIncrement" :0.01,
264 "withAvoidingRedundancy" :True,
265 "withToleranceInRedundancy" :1.e-18,
266 "withLenghtOfRedundancy" :-1,
267 "withmpEnabled" :False,
268 "withmpWorkers" :None,
271 if (type(asFunction) is type({})) and \
272 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
273 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
274 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
275 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
276 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
277 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
278 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
279 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
280 if not asFunction.has_key("withmpEnabled"): asFunction["withmpEnabled"] = False
281 if not asFunction.has_key("withmpWorkers"): asFunction["withmpWorkers"] = None
282 from daNumerics.ApproximatedDerivatives import FDApproximation
283 FDA = FDApproximation(
284 Function = asFunction["Direct"],
285 centeredDF = asFunction["withCenteredDF"],
286 increment = asFunction["withIncrement"],
287 dX = asFunction["withdX"],
288 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
289 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
290 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
291 mpEnabled = asFunction["withmpEnabled"],
292 mpWorkers = asFunction["withmpWorkers"],
294 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator )
295 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
296 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
297 elif (type(asFunction) is type({})) and \
298 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
299 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
300 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
301 self.__HO["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
303 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"] )
304 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
305 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
306 elif asMatrix is not None:
307 matrice = numpy.matrix( asMatrix, numpy.float )
308 self.__HO["Direct"] = Operator( fromMatrix = matrice )
309 self.__HO["Tangent"] = Operator( fromMatrix = matrice )
310 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T )
313 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
315 if appliedToX is not None:
316 self.__HO["AppliedToX"] = {}
317 if type(appliedToX) is not dict:
318 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
319 for key in appliedToX.keys():
320 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
321 # Pour le cas où l'on a une vraie matrice
322 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
323 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
324 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
325 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
327 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
329 self.__HO["AppliedToX"] = None
332 self.__StoredInputs["ObservationOperator"] = self.__HO
335 # -----------------------------------------------------------
336 def setEvolutionModel(self,
343 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
344 définitions et leur sens sont les suivants :
345 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
346 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
347 "Direct" n'est pas définie, on prend la fonction "Tangent".
348 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
349 des opérateurs tangents et adjoints. On utilise par défaut des
350 différences finies non centrées ou centrées (si "withCenteredDF" est
351 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
352 du point courant ou sur le point fixe "withdX".
353 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
354 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
355 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
356 matrice fournie doit être sous une forme compatible avec le
357 constructeur de numpy.matrix.
358 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
359 être rendue disponible au même titre que les variables de calcul
360 L'argument "asFunction" peut prendre la forme complète suivante, avec
361 les valeurs par défaut standards :
362 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
363 "useApproximatedDerivatives":False,
364 "withCenteredDF" :False,
365 "withIncrement" :0.01,
367 "withAvoidingRedundancy" :True,
368 "withToleranceInRedundancy" :1.e-18,
369 "withLenghtOfRedundancy" :-1,
370 "withmpEnabled" :False,
371 "withmpWorkers" :None,
374 if (type(asFunction) is type({})) and \
375 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
376 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
377 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
378 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
379 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
380 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
381 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
382 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
383 if not asFunction.has_key("withmpEnabled"): asFunction["withmpEnabled"] = False
384 if not asFunction.has_key("withmpWorkers"): asFunction["withmpWorkers"] = None
385 from daNumerics.ApproximatedDerivatives import FDApproximation
386 FDA = FDApproximation(
387 Function = asFunction["Direct"],
388 centeredDF = asFunction["withCenteredDF"],
389 increment = asFunction["withIncrement"],
390 dX = asFunction["withdX"],
391 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
392 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
393 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
394 mpEnabled = asFunction["withmpEnabled"],
395 mpWorkers = asFunction["withmpWorkers"],
397 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
398 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
399 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
400 elif (type(asFunction) is type({})) and \
401 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
402 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
403 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
404 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
406 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
407 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
408 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
409 elif asMatrix is not None:
410 matrice = numpy.matrix( asMatrix, numpy.float )
411 self.__EM["Direct"] = Operator( fromMatrix = matrice )
412 self.__EM["Tangent"] = Operator( fromMatrix = matrice )
413 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T )
416 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
419 self.__StoredInputs["EvolutionModel"] = self.__EM
422 def setEvolutionError(self,
424 asEyeByScalar = None,
425 asEyeByVector = None,
429 Permet de définir la covariance des erreurs de modèle :
430 - asCovariance : entrée des données, comme une matrice compatible avec
431 le constructeur de numpy.matrix
432 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
433 multiplicatif d'une matrice de corrélation identité, aucune matrice
434 n'étant donc explicitement à donner
435 - asEyeByVector : entrée des données comme un seul vecteur de variance,
436 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
437 n'étant donc explicitement à donner
438 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
439 être rendue disponible au même titre que les variables de calcul
441 self.__Q = Covariance(
442 name = "ObservationError",
443 asCovariance = asCovariance,
444 asEyeByScalar = asEyeByScalar,
445 asEyeByVector = asEyeByVector,
448 self.__StoredInputs["EvolutionError"] = self.__Q
451 # -----------------------------------------------------------
452 def setControlModel(self,
453 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
454 "useApproximatedDerivatives":False,
455 "withCenteredDF" :False,
456 "withIncrement" :0.01,
464 Permet de définir un opérateur de controle C. L'ordre de priorité des
465 définitions et leur sens sont les suivants :
466 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
467 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
468 "Direct" n'est pas définie, on prend la fonction "Tangent".
469 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
470 des opérateurs tangents et adjoints. On utilise par défaut des
471 différences finies non centrées ou centrées (si "withCenteredDF" est
472 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
473 du point courant ou sur le point fixe "withdX".
474 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
475 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
476 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
477 matrice fournie doit être sous une forme compatible avec le
478 constructeur de numpy.matrix.
479 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
480 être rendue disponible au même titre que les variables de calcul
482 if (type(asFunction) is type({})) and \
483 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
484 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
485 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
486 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
487 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
488 from daNumerics.ApproximatedDerivatives import FDApproximation
489 FDA = FDApproximation(
490 Function = asFunction["Direct"],
491 centeredDF = asFunction["withCenteredDF"],
492 increment = asFunction["withIncrement"],
493 dX = asFunction["withdX"] )
494 self.__CM["Direct"] = Operator( fromMethod = FDA.DirectOperator )
495 self.__CM["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
496 self.__CM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
497 elif (type(asFunction) is type({})) and \
498 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
499 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
500 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
501 self.__CM["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
503 self.__CM["Direct"] = Operator( fromMethod = asFunction["Direct"] )
504 self.__CM["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
505 self.__CM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
506 elif asMatrix is not None:
507 matrice = numpy.matrix( asMatrix, numpy.float )
508 self.__CM["Direct"] = Operator( fromMatrix = matrice )
509 self.__CM["Tangent"] = Operator( fromMatrix = matrice )
510 self.__CM["Adjoint"] = Operator( fromMatrix = matrice.T )
513 raise ValueError("Improperly defined input control model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
516 self.__StoredInputs["ControlModel"] = self.__CM
519 def setControlInput(self,
521 asPersistentVector = None,
526 Permet de définir le controle en entree :
527 - asVector : entrée des données, comme un vecteur compatible avec le
528 constructeur de numpy.matrix
529 - asPersistentVector : entrée des données, comme un vecteur de type
530 persistent contruit avec la classe ad-hoc "Persistence"
531 - Scheduler est le contrôle temporel des données disponibles
532 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
533 être rendue disponible au même titre que les variables de calcul
535 if asVector is not None:
536 if isinstance(asVector,numpy.matrix):
537 self.__U = numpy.matrix( asVector.A1, numpy.float ).T
539 self.__U = numpy.matrix( asVector, numpy.float ).T
540 elif asPersistentVector is not None:
541 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
542 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.matrix)
543 for member in asPersistentVector:
544 self.__U.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
546 self.__U = asPersistentVector
548 raise ValueError("Error: improperly defined control input, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
550 self.__StoredInputs["ControlInput"] = self.__U
553 # -----------------------------------------------------------
554 def setControls (self,
559 Permet de définir la valeur initiale du vecteur X contenant toutes les
560 variables de contrôle, i.e. les paramètres ou l'état dont on veut
561 estimer la valeur pour obtenir les observations. C'est utile pour un
562 algorithme itératif/incrémental.
563 - asVector : entrée des données, comme un vecteur compatible avec le
564 constructeur de numpy.matrix.
565 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
566 être rendue disponible au même titre que les variables de calcul
568 if asVector is not None:
569 self.__X.store( asVector )
571 self.__StoredInputs["Controls"] = self.__X
574 # -----------------------------------------------------------
575 def setAlgorithm(self, choice = None ):
577 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
578 d'assimilation. L'argument est un champ caractère se rapportant au nom
579 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
580 d'assimilation sur les arguments fixes.
583 raise ValueError("Error: algorithm choice has to be given")
584 if self.__algorithmName is not None:
585 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
586 daDirectory = "daAlgorithms"
588 # Recherche explicitement le fichier complet
589 # ------------------------------------------
591 for directory in sys.path:
592 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
593 module_path = os.path.abspath(os.path.join(directory, daDirectory))
594 if module_path is None:
595 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
597 # Importe le fichier complet comme un module
598 # ------------------------------------------
600 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
601 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
602 self.__algorithmName = str(choice)
603 sys.path = sys_path_tmp ; del sys_path_tmp
604 except ImportError, e:
605 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
607 # Instancie un objet du type élémentaire du fichier
608 # -------------------------------------------------
609 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
610 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
613 def setAlgorithmParameters(self, asDico=None):
615 Permet de définir les paramètres de l'algorithme, sous la forme d'un
618 if asDico is not None:
619 self.__Parameters.update( dict( asDico ) )
621 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
624 def getAlgorithmParameters(self, noDetails=True):
626 Renvoie la liste des paramètres requis selon l'algorithme
628 return self.__algorithm.getRequiredParameters(noDetails)
630 # -----------------------------------------------------------
631 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
633 Permet de sélectionner un diagnostic a effectuer.
636 raise ValueError("Error: diagnostic choice has to be given")
637 daDirectory = "daDiagnostics"
639 # Recherche explicitement le fichier complet
640 # ------------------------------------------
642 for directory in sys.path:
643 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
644 module_path = os.path.abspath(os.path.join(directory, daDirectory))
645 if module_path is None:
646 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
648 # Importe le fichier complet comme un module
649 # ------------------------------------------
651 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
652 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
653 sys.path = sys_path_tmp ; del sys_path_tmp
654 except ImportError, e:
655 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
657 # Instancie un objet du type élémentaire du fichier
658 # -------------------------------------------------
659 if self.__StoredInputs.has_key(name):
660 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
661 elif self.__StoredDiagnostics.has_key(name):
662 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
664 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
668 parameters = parameters )
671 # -----------------------------------------------------------
672 def shape_validate(self):
674 Validation de la correspondance correcte des tailles des variables et
675 des matrices s'il y en a.
677 if self.__Xb is None: __Xb_shape = (0,)
678 elif hasattr(self.__Xb,"size"): __Xb_shape = (self.__Xb.size,)
679 elif hasattr(self.__Xb,"shape"):
680 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
681 else: __Xb_shape = self.__Xb.shape()
682 else: raise TypeError("The background (Xb) has no attribute of shape: problem !")
684 if self.__Y is None: __Y_shape = (0,)
685 elif hasattr(self.__Y,"size"): __Y_shape = (self.__Y.size,)
686 elif hasattr(self.__Y,"shape"):
687 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
688 else: __Y_shape = self.__Y.shape()
689 else: raise TypeError("The observation (Y) has no attribute of shape: problem !")
691 if self.__U is None: __U_shape = (0,)
692 elif hasattr(self.__U,"size"): __U_shape = (self.__U.size,)
693 elif hasattr(self.__U,"shape"):
694 if type(self.__U.shape) is tuple: __U_shape = self.__U.shape
695 else: __U_shape = self.__U.shape()
696 else: raise TypeError("The control (U) has no attribute of shape: problem !")
698 if self.__B is None: __B_shape = (0,0)
699 elif hasattr(self.__B,"shape"):
700 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
701 else: __B_shape = self.__B.shape()
702 else: raise TypeError("The a priori errors covariance matrix (B) has no attribute of shape: problem !")
704 if self.__R is None: __R_shape = (0,0)
705 elif hasattr(self.__R,"shape"):
706 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
707 else: __R_shape = self.__R.shape()
708 else: raise TypeError("The observation errors covariance matrix (R) has no attribute of shape: problem !")
710 if self.__Q is None: __Q_shape = (0,0)
711 elif hasattr(self.__Q,"shape"):
712 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
713 else: __Q_shape = self.__Q.shape()
714 else: raise TypeError("The evolution errors covariance matrix (Q) has no attribute of shape: problem !")
716 if len(self.__HO) == 0: __HO_shape = (0,0)
717 elif type(self.__HO) is type({}): __HO_shape = (0,0)
718 elif hasattr(self.__HO["Direct"],"shape"):
719 if type(self.__HO["Direct"].shape) is tuple: __HO_shape = self.__HO["Direct"].shape
720 else: __HO_shape = self.__HO["Direct"].shape()
721 else: raise TypeError("The observation operator (H) has no attribute of shape: problem !")
723 if len(self.__EM) == 0: __EM_shape = (0,0)
724 elif type(self.__EM) is type({}): __EM_shape = (0,0)
725 elif hasattr(self.__EM["Direct"],"shape"):
726 if type(self.__EM["Direct"].shape) is tuple: __EM_shape = self.__EM["Direct"].shape
727 else: __EM_shape = self.__EM["Direct"].shape()
728 else: raise TypeError("The evolution model (EM) has no attribute of shape: problem !")
730 if len(self.__CM) == 0: __CM_shape = (0,0)
731 elif type(self.__CM) is type({}): __CM_shape = (0,0)
732 elif hasattr(self.__CM["Direct"],"shape"):
733 if type(self.__CM["Direct"].shape) is tuple: __CM_shape = self.__CM["Direct"].shape
734 else: __CM_shape = self.__CM["Direct"].shape()
735 else: raise TypeError("The control model (CM) has no attribute of shape: problem !")
737 # Vérification des conditions
738 # ---------------------------
739 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
740 raise ValueError("Shape characteristic of background (Xb) is incorrect: \"%s\"."%(__Xb_shape,))
741 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
742 raise ValueError("Shape characteristic of observation (Y) is incorrect: \"%s\"."%(__Y_shape,))
744 if not( min(__B_shape) == max(__B_shape) ):
745 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) is incorrect: \"%s\"."%(__B_shape,))
746 if not( min(__R_shape) == max(__R_shape) ):
747 raise ValueError("Shape characteristic of observation errors covariance matrix (R) is incorrect: \"%s\"."%(__R_shape,))
748 if not( min(__Q_shape) == max(__Q_shape) ):
749 raise ValueError("Shape characteristic of evolution errors covariance matrix (Q) is incorrect: \"%s\"."%(__Q_shape,))
750 if not( min(__EM_shape) == max(__EM_shape) ):
751 raise ValueError("Shape characteristic of evolution operator (EM) is incorrect: \"%s\"."%(__EM_shape,))
753 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
754 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and state (X) \"%s\" are incompatible."%(__HO_shape,__Xb_shape))
755 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
756 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation (Y) \"%s\" are incompatible."%(__HO_shape,__Y_shape))
757 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] ):
758 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and a priori errors covariance matrix (B) \"%s\" are incompatible."%(__HO_shape,__B_shape))
759 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] ):
760 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation errors covariance matrix (R) \"%s\" are incompatible."%(__HO_shape,__R_shape))
762 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
763 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
764 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
765 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
766 for member in asPersistentVector:
767 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
768 __Xb_shape = min(__B_shape)
770 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) \"%s\" and background (Xb) \"%s\" are incompatible."%(__B_shape,__Xb_shape))
772 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
773 raise ValueError("Shape characteristic of observation errors covariance matrix (R) \"%s\" and observation (Y) \"%s\" are incompatible."%(__R_shape,__Y_shape))
775 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) ):
776 raise ValueError("Shape characteristic of evolution model (EM) \"%s\" and state (X) \"%s\" are incompatible."%(__EM_shape,__Xb_shape))
778 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) ):
779 raise ValueError("Shape characteristic of control model (CM) \"%s\" and control (U) \"%s\" are incompatible."%(__CM_shape,__U_shape))
781 if self.__StoredInputs.has_key("AlgorithmParameters") \
782 and self.__StoredInputs["AlgorithmParameters"].has_key("Bounds") \
783 and (type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type([]) or type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type(())) \
784 and (len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) != max(__Xb_shape)):
785 raise ValueError("The number \"%s\" of bound pairs for the state (X) components is different of the size \"%s\" of the state itself." \
786 %(len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]),max(__Xb_shape)))
790 # -----------------------------------------------------------
793 Permet de lancer le calcul d'assimilation.
795 Le nom de la méthode à activer est toujours "run". Les paramètres en
796 arguments de la méthode sont fixés. En sortie, on obtient les résultats
797 dans la variable de type dictionnaire "StoredVariables", qui contient en
798 particulier des objets de Persistence pour les analyses, OMA...
800 Operator.CM.clearCache()
801 self.shape_validate()
803 self.__algorithm.run(
813 Parameters = self.__Parameters,
817 # -----------------------------------------------------------
818 def get(self, key=None):
820 Renvoie les résultats disponibles après l'exécution de la méthode
821 d'assimilation, ou les diagnostics disponibles. Attention, quand un
822 diagnostic porte le même nom qu'une variable stockée, c'est la variable
823 stockée qui est renvoyée, et le diagnostic est inatteignable.
826 if self.__algorithm.has_key(key):
827 return self.__algorithm.get( key )
828 elif self.__StoredInputs.has_key(key):
829 return self.__StoredInputs[key]
830 elif self.__StoredDiagnostics.has_key(key):
831 return self.__StoredDiagnostics[key]
833 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
835 allvariables = self.__algorithm.get()
836 allvariables.update( self.__StoredDiagnostics )
837 allvariables.update( self.__StoredInputs )
840 def get_available_variables(self):
842 Renvoie les variables potentiellement utilisables pour l'étude,
843 initialement stockées comme données d'entrées ou dans les algorithmes,
844 identifiés par les chaînes de caractères. L'algorithme doit avoir été
845 préalablement choisi sinon la méthode renvoie "None".
847 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
851 if len( self.__algorithm.keys()) > 0:
852 variables.extend( self.__algorithm.get().keys() )
853 if len( self.__StoredInputs.keys() ) > 0:
854 variables.extend( self.__StoredInputs.keys() )
858 def get_available_algorithms(self):
860 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
861 par les chaînes de caractères.
864 for directory in sys.path:
865 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
866 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
867 root, ext = os.path.splitext(fname)
868 if ext == '.py' and root != '__init__':
873 def get_available_diagnostics(self):
875 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
876 par les chaînes de caractères.
879 for directory in sys.path:
880 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
881 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
882 root, ext = os.path.splitext(fname)
883 if ext == '.py' and root != '__init__':
888 # -----------------------------------------------------------
889 def get_algorithms_main_path(self):
891 Renvoie le chemin pour le répertoire principal contenant les algorithmes
892 dans un sous-répertoire "daAlgorithms"
896 def add_algorithms_path(self, asPath=None):
898 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
899 se trouve un sous-répertoire "daAlgorithms"
901 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
902 pas indispensable de le rajouter ici.
904 if not os.path.isdir(asPath):
905 raise ValueError("The given "+asPath+" argument must exist as a directory")
906 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
907 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
908 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
909 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
910 sys.path.insert(0, os.path.abspath(asPath))
911 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
914 def get_diagnostics_main_path(self):
916 Renvoie le chemin pour le répertoire principal contenant les diagnostics
917 dans un sous-répertoire "daDiagnostics"
921 def add_diagnostics_path(self, asPath=None):
923 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
924 se trouve un sous-répertoire "daDiagnostics"
926 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
927 pas indispensable de le rajouter ici.
929 if not os.path.isdir(asPath):
930 raise ValueError("The given "+asPath+" argument must exist as a directory")
931 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
932 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
933 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
934 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
935 sys.path.insert(0, os.path.abspath(asPath))
936 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
939 # -----------------------------------------------------------
940 def setDataObserver(self,
943 HookParameters = None,
947 Permet d'associer un observer à une ou des variables nommées gérées en
948 interne, activable selon des règles définies dans le Scheduler. A chaque
949 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
950 les arguments (variable persistante VariableName, paramètres HookParameters).
953 if type( self.__algorithm ) is dict:
954 raise ValueError("No observer can be build before choosing an algorithm.")
956 # Vérification du nom de variable et typage
957 # -----------------------------------------
958 if type( VariableName ) is str:
959 VariableNames = [VariableName,]
960 elif type( VariableName ) is list:
961 VariableNames = map( str, VariableName )
963 raise ValueError("The observer requires a name or a list of names of variables.")
965 # Association interne de l'observer à la variable
966 # -----------------------------------------------
967 for n in VariableNames:
968 if not self.__algorithm.has_key( n ):
969 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
971 self.__algorithm.StoredVariables[ n ].setDataObserver(
972 Scheduler = Scheduler,
973 HookFunction = HookFunction,
974 HookParameters = HookParameters,
977 def removeDataObserver(self,
982 Permet de retirer un observer à une ou des variable nommée.
985 if type( self.__algorithm ) is dict:
986 raise ValueError("No observer can be removed before choosing an algorithm.")
988 # Vérification du nom de variable et typage
989 # -----------------------------------------
990 if type( VariableName ) is str:
991 VariableNames = [VariableName,]
992 elif type( VariableName ) is list:
993 VariableNames = map( str, VariableName )
995 raise ValueError("The observer requires a name or a list of names of variables.")
997 # Association interne de l'observer à la variable
998 # -----------------------------------------------
999 for n in VariableNames:
1000 if not self.__algorithm.has_key( n ):
1001 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
1003 self.__algorithm.StoredVariables[ n ].removeDataObserver(
1004 HookFunction = HookFunction,
1007 # -----------------------------------------------------------
1008 def setDebug(self, level=10):
1010 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
1011 appel pour changer le niveau de verbosité, avec :
1012 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
1014 log = logging.getLogger()
1015 log.setLevel( level )
1017 def unsetDebug(self):
1019 Remet le logger au niveau par défaut
1021 log = logging.getLogger()
1022 log.setLevel( logging.WARNING )
1024 def prepare_to_pickle(self):
1026 Retire les variables non pickelisables
1028 self.__algorithmFile = None
1029 self.__diagnosticFile = None
1034 # ==============================================================================
1035 if __name__ == "__main__":
1036 print '\n AUTODIAGNOSTIC \n'