1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2015 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"
31 __all__ = ["AssimilationStudy"]
35 import ExtendedLogging ; ExtendedLogging.ExtendedLogging() # A importer en premier
39 logging.debug("Succeed initial import of scipy.optimize with Scipy %s", scipy.version.version)
41 logging.debug("Fail initial import of scipy.optimize")
43 from BasicObjects import Operator, Covariance
44 from PlatformInfo import uniq
46 # ==============================================================================
47 class AssimilationStudy:
49 Cette classe sert d'interface pour l'utilisation de l'assimilation de
50 données. Elle contient les méthodes ou accesseurs nécessaires à la
51 construction d'un calcul d'assimilation.
53 def __init__(self, name=""):
55 Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
56 élémentaire. Ces variables sont ensuite disponibles pour implémenter un
57 algorithme élémentaire particulier.
59 - Background : vecteur Xb
60 - Observation : vecteur Y (potentiellement temporel) d'observations
61 - State : vecteur d'état dont une partie est le vecteur de contrôle.
62 Cette information n'est utile que si l'on veut faire des calculs sur
63 l'état complet, mais elle n'est pas indispensable pour l'assimilation.
64 - Control : vecteur X contenant toutes les variables de contrôle, i.e.
65 les paramètres ou l'état dont on veut estimer la valeur pour obtenir
67 - ObservationOperator...: opérateur d'observation H
69 Les observations présentent une erreur dont la matrice de covariance est
70 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
73 self.__name = str(name)
84 self.__X = Persistence.OneVector()
85 self.__Parameters = {}
86 self.__StoredDiagnostics = {}
87 self.__StoredInputs = {}
89 # Variables temporaires
91 self.__algorithmFile = None
92 self.__algorithmName = None
93 self.__diagnosticFile = None
95 # Récupère le chemin du répertoire parent et l'ajoute au path
96 # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
97 # qui est activée dans Persistence)
98 self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
99 sys.path.insert(0, self.__parent)
100 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
102 # ---------------------------------------------------------
103 def setBackground(self,
105 asPersistentVector = None,
110 Permet de définir l'estimation a priori :
111 - asVector : entrée des données, comme un vecteur compatible avec le
112 constructeur de numpy.matrix
113 - asPersistentVector : entrée des données, comme un vecteur de type
114 persistent contruit avec la classe ad-hoc "Persistence"
115 - Scheduler est le contrôle temporel des données
116 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
117 être rendue disponible au même titre que les variables de calcul
119 if asVector is not None:
120 if type( asVector ) is type( numpy.matrix([]) ):
121 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
123 self.__Xb = numpy.matrix( asVector, numpy.float ).T
124 elif asPersistentVector is not None:
125 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
126 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
127 for member in asPersistentVector:
128 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
130 self.__Xb = asPersistentVector
132 raise ValueError("Error: improperly defined background, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
134 self.__StoredInputs["Background"] = self.__Xb
137 def setBackgroundError(self,
139 asEyeByScalar = None,
140 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 - asCovObject : entrée des données comme un objet ayant des méthodes
155 particulieres de type matriciel
156 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
157 être rendue disponible au même titre que les variables de calcul
159 self.__B = Covariance(
160 name = "BackgroundError",
161 asCovariance = asCovariance,
162 asEyeByScalar = asEyeByScalar,
163 asEyeByVector = asEyeByVector,
164 asCovObject = asCovObject,
167 self.__StoredInputs["BackgroundError"] = self.__B
170 # -----------------------------------------------------------
171 def setObservation(self,
173 asPersistentVector = None,
178 Permet de définir les observations :
179 - asVector : entrée des données, comme un vecteur compatible avec le
180 constructeur de numpy.matrix
181 - asPersistentVector : entrée des données, comme un vecteur de type
182 persistent contruit avec la classe ad-hoc "Persistence"
183 - Scheduler est le contrôle temporel des données disponibles
184 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
185 être rendue disponible au même titre que les variables de calcul
187 if asVector is not None:
188 if type( asVector ) is type( numpy.matrix([]) ):
189 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
191 self.__Y = numpy.matrix( asVector, numpy.float ).T
192 elif asPersistentVector is not None:
193 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
194 self.__Y = Persistence.OneVector("Observation", basetype=numpy.matrix)
195 for member in asPersistentVector:
196 self.__Y.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
198 self.__Y = asPersistentVector
200 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
202 self.__StoredInputs["Observation"] = self.__Y
205 def setObservationError(self,
207 asEyeByScalar = None,
208 asEyeByVector = None,
213 Permet de définir la covariance des erreurs d'observations :
214 - asCovariance : entrée des données, comme une matrice compatible avec
215 le constructeur de numpy.matrix
216 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
217 multiplicatif d'une matrice de corrélation identité, aucune matrice
218 n'étant donc explicitement à donner
219 - asEyeByVector : entrée des données comme un seul vecteur de variance,
220 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
221 n'étant donc explicitement à donner
222 - asCovObject : entrée des données comme un objet ayant des méthodes
223 particulieres de type matriciel
224 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
225 être rendue disponible au même titre que les variables de calcul
227 self.__R = Covariance(
228 name = "ObservationError",
229 asCovariance = asCovariance,
230 asEyeByScalar = asEyeByScalar,
231 asEyeByVector = asEyeByVector,
232 asCovObject = asCovObject,
235 self.__StoredInputs["ObservationError"] = self.__R
238 def setObservationOperator(self,
246 Permet de définir un opérateur d'observation H. L'ordre de priorité des
247 définitions et leur sens sont les suivants :
248 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
249 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
250 "Direct" n'est pas définie, on prend la fonction "Tangent".
251 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
252 des opérateurs tangents et adjoints. On utilise par défaut des
253 différences finies non centrées ou centrées (si "withCenteredDF" est
254 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
255 du point courant ou sur le point fixe "withdX".
256 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
257 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
258 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
259 matrice fournie doit être sous une forme compatible avec le
260 constructeur de numpy.matrix.
261 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
262 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
263 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
264 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
265 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
266 être rendue disponible au même titre que les variables de calcul
267 L'argument "asFunction" peut prendre la forme complète suivante, avec
268 les valeurs par défaut standards :
269 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
270 "useApproximatedDerivatives":False,
271 "withCenteredDF" :False,
272 "withIncrement" :0.01,
274 "withAvoidingRedundancy" :True,
275 "withToleranceInRedundancy" :1.e-18,
276 "withLenghtOfRedundancy" :-1,
277 "withmpEnabled" :False,
278 "withmpWorkers" :None,
281 if (type(asFunction) is type({})) and \
282 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
283 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
284 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
285 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
286 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
287 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
288 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
289 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
290 if not asFunction.has_key("withmpEnabled"): asFunction["withmpEnabled"] = False
291 if not asFunction.has_key("withmpWorkers"): asFunction["withmpWorkers"] = None
292 from daNumerics.ApproximatedDerivatives import FDApproximation
293 FDA = FDApproximation(
294 Function = asFunction["Direct"],
295 centeredDF = asFunction["withCenteredDF"],
296 increment = asFunction["withIncrement"],
297 dX = asFunction["withdX"],
298 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
299 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
300 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
301 mpEnabled = asFunction["withmpEnabled"],
302 mpWorkers = asFunction["withmpWorkers"],
304 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
305 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
306 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
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"], avoidingRedundancy = avoidRC )
313 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
314 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
315 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
316 elif asMatrix is not None:
317 matrice = numpy.matrix( asMatrix, numpy.float )
318 self.__HO["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
319 self.__HO["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
320 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
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,
354 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
355 définitions et leur sens sont les suivants :
356 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
357 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
358 "Direct" n'est pas définie, on prend la fonction "Tangent".
359 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
360 des opérateurs tangents et adjoints. On utilise par défaut des
361 différences finies non centrées ou centrées (si "withCenteredDF" est
362 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
363 du point courant ou sur le point fixe "withdX".
364 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
365 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
366 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
367 matrice fournie doit être sous une forme compatible avec le
368 constructeur de numpy.matrix.
369 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
370 être rendue disponible au même titre que les variables de calcul
371 L'argument "asFunction" peut prendre la forme complète suivante, avec
372 les valeurs par défaut standards :
373 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
374 "useApproximatedDerivatives":False,
375 "withCenteredDF" :False,
376 "withIncrement" :0.01,
378 "withAvoidingRedundancy" :True,
379 "withToleranceInRedundancy" :1.e-18,
380 "withLenghtOfRedundancy" :-1,
381 "withmpEnabled" :False,
382 "withmpWorkers" :None,
385 if (type(asFunction) is type({})) and \
386 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
387 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
388 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
389 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
390 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
391 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
392 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
393 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
394 if not asFunction.has_key("withmpEnabled"): asFunction["withmpEnabled"] = False
395 if not asFunction.has_key("withmpWorkers"): asFunction["withmpWorkers"] = None
396 from daNumerics.ApproximatedDerivatives import FDApproximation
397 FDA = FDApproximation(
398 Function = asFunction["Direct"],
399 centeredDF = asFunction["withCenteredDF"],
400 increment = asFunction["withIncrement"],
401 dX = asFunction["withdX"],
402 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
403 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
404 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
405 mpEnabled = asFunction["withmpEnabled"],
406 mpWorkers = asFunction["withmpWorkers"],
408 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
409 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
410 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
411 elif (type(asFunction) is type({})) and \
412 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
413 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
414 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
415 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
417 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
418 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
419 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
420 elif asMatrix is not None:
421 matrice = numpy.matrix( asMatrix, numpy.float )
422 self.__EM["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
423 self.__EM["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
424 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
427 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
430 self.__StoredInputs["EvolutionModel"] = self.__EM
433 def setEvolutionError(self,
435 asEyeByScalar = None,
436 asEyeByVector = None,
441 Permet de définir la covariance des erreurs de modèle :
442 - asCovariance : entrée des données, comme une matrice compatible avec
443 le constructeur de numpy.matrix
444 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
445 multiplicatif d'une matrice de corrélation identité, aucune matrice
446 n'étant donc explicitement à donner
447 - asEyeByVector : entrée des données comme un seul vecteur de variance,
448 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
449 n'étant donc explicitement à donner
450 - asCovObject : entrée des données comme un objet ayant des méthodes
451 particulieres de type matriciel
452 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
453 être rendue disponible au même titre que les variables de calcul
455 self.__Q = Covariance(
456 name = "EvolutionError",
457 asCovariance = asCovariance,
458 asEyeByScalar = asEyeByScalar,
459 asEyeByVector = asEyeByVector,
460 asCovObject = asCovObject,
463 self.__StoredInputs["EvolutionError"] = self.__Q
466 # -----------------------------------------------------------
467 def setControlModel(self,
468 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
469 "useApproximatedDerivatives":False,
470 "withCenteredDF" :False,
471 "withIncrement" :0.01,
480 Permet de définir un opérateur de controle C. L'ordre de priorité des
481 définitions et leur sens sont les suivants :
482 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
483 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
484 "Direct" n'est pas définie, on prend la fonction "Tangent".
485 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
486 des opérateurs tangents et adjoints. On utilise par défaut des
487 différences finies non centrées ou centrées (si "withCenteredDF" est
488 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
489 du point courant ou sur le point fixe "withdX".
490 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
491 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
492 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
493 matrice fournie doit être sous une forme compatible avec le
494 constructeur de numpy.matrix.
495 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
496 être rendue disponible au même titre que les variables de calcul
498 if (type(asFunction) is type({})) and \
499 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
500 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
501 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
502 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
503 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
504 from daNumerics.ApproximatedDerivatives import FDApproximation
505 FDA = FDApproximation(
506 Function = asFunction["Direct"],
507 centeredDF = asFunction["withCenteredDF"],
508 increment = asFunction["withIncrement"],
509 dX = asFunction["withdX"] )
510 self.__CM["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
511 self.__CM["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
512 self.__CM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
513 elif (type(asFunction) is type({})) and \
514 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
515 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
516 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
517 self.__CM["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
519 self.__CM["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
520 self.__CM["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
521 self.__CM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
522 elif asMatrix is not None:
523 matrice = numpy.matrix( asMatrix, numpy.float )
524 self.__CM["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
525 self.__CM["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
526 self.__CM["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
529 raise ValueError("Improperly defined input control model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
532 self.__StoredInputs["ControlModel"] = self.__CM
535 def setControlInput(self,
537 asPersistentVector = None,
542 Permet de définir le controle en entree :
543 - asVector : entrée des données, comme un vecteur compatible avec le
544 constructeur de numpy.matrix
545 - asPersistentVector : entrée des données, comme un vecteur de type
546 persistent contruit avec la classe ad-hoc "Persistence"
547 - Scheduler est le contrôle temporel des données disponibles
548 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
549 être rendue disponible au même titre que les variables de calcul
551 if asVector is not None:
552 if isinstance(asVector,numpy.matrix):
553 self.__U = numpy.matrix( asVector.A1, numpy.float ).T
555 self.__U = numpy.matrix( asVector, numpy.float ).T
556 elif asPersistentVector is not None:
557 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
558 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.matrix)
559 for member in asPersistentVector:
560 self.__U.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
562 self.__U = asPersistentVector
564 raise ValueError("Error: improperly defined control input, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
566 self.__StoredInputs["ControlInput"] = self.__U
569 # -----------------------------------------------------------
570 def setControls (self,
575 Permet de définir la valeur initiale du vecteur X contenant toutes les
576 variables de contrôle, i.e. les paramètres ou l'état dont on veut
577 estimer la valeur pour obtenir les observations. C'est utile pour un
578 algorithme itératif/incrémental.
579 - asVector : entrée des données, comme un vecteur compatible avec le
580 constructeur de numpy.matrix.
581 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
582 être rendue disponible au même titre que les variables de calcul
584 if asVector is not None:
585 self.__X.store( asVector )
587 self.__StoredInputs["Controls"] = self.__X
590 # -----------------------------------------------------------
591 def setAlgorithm(self, choice = None ):
593 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
594 d'assimilation. L'argument est un champ caractère se rapportant au nom
595 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
596 d'assimilation sur les arguments fixes.
599 raise ValueError("Error: algorithm choice has to be given")
600 if self.__algorithmName is not None:
601 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
602 daDirectory = "daAlgorithms"
604 # Recherche explicitement le fichier complet
605 # ------------------------------------------
607 for directory in sys.path:
608 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
609 module_path = os.path.abspath(os.path.join(directory, daDirectory))
610 if module_path is None:
611 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
613 # Importe le fichier complet comme un module
614 # ------------------------------------------
616 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
617 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
618 self.__algorithmName = str(choice)
619 sys.path = sys_path_tmp ; del sys_path_tmp
620 except ImportError, e:
621 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
623 # Instancie un objet du type élémentaire du fichier
624 # -------------------------------------------------
625 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
626 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
629 def setAlgorithmParameters(self, asDico=None):
631 Permet de définir les paramètres de l'algorithme, sous la forme d'un
634 if asDico is not None:
635 self.__Parameters.update( dict( asDico ) )
637 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
640 def getAlgorithmParameters(self, noDetails=True):
642 Renvoie la liste des paramètres requis selon l'algorithme
644 return self.__algorithm.getRequiredParameters(noDetails)
646 # -----------------------------------------------------------
647 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
649 Permet de sélectionner un diagnostic a effectuer.
652 raise ValueError("Error: diagnostic choice has to be given")
653 daDirectory = "daDiagnostics"
655 # Recherche explicitement le fichier complet
656 # ------------------------------------------
658 for directory in sys.path:
659 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
660 module_path = os.path.abspath(os.path.join(directory, daDirectory))
661 if module_path is None:
662 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
664 # Importe le fichier complet comme un module
665 # ------------------------------------------
667 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
668 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
669 sys.path = sys_path_tmp ; del sys_path_tmp
670 except ImportError, e:
671 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
673 # Instancie un objet du type élémentaire du fichier
674 # -------------------------------------------------
675 if self.__StoredInputs.has_key(name):
676 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
677 elif self.__StoredDiagnostics.has_key(name):
678 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
680 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
684 parameters = parameters )
687 # -----------------------------------------------------------
688 def shape_validate(self):
690 Validation de la correspondance correcte des tailles des variables et
691 des matrices s'il y en a.
693 if self.__Xb is None: __Xb_shape = (0,)
694 elif hasattr(self.__Xb,"size"): __Xb_shape = (self.__Xb.size,)
695 elif hasattr(self.__Xb,"shape"):
696 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
697 else: __Xb_shape = self.__Xb.shape()
698 else: raise TypeError("The background (Xb) has no attribute of shape: problem !")
700 if self.__Y is None: __Y_shape = (0,)
701 elif hasattr(self.__Y,"size"): __Y_shape = (self.__Y.size,)
702 elif hasattr(self.__Y,"shape"):
703 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
704 else: __Y_shape = self.__Y.shape()
705 else: raise TypeError("The observation (Y) has no attribute of shape: problem !")
707 if self.__U is None: __U_shape = (0,)
708 elif hasattr(self.__U,"size"): __U_shape = (self.__U.size,)
709 elif hasattr(self.__U,"shape"):
710 if type(self.__U.shape) is tuple: __U_shape = self.__U.shape
711 else: __U_shape = self.__U.shape()
712 else: raise TypeError("The control (U) has no attribute of shape: problem !")
714 if self.__B is None: __B_shape = (0,0)
715 elif hasattr(self.__B,"shape"):
716 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
717 else: __B_shape = self.__B.shape()
718 else: raise TypeError("The a priori errors covariance matrix (B) has no attribute of shape: problem !")
720 if self.__R is None: __R_shape = (0,0)
721 elif hasattr(self.__R,"shape"):
722 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
723 else: __R_shape = self.__R.shape()
724 else: raise TypeError("The observation errors covariance matrix (R) has no attribute of shape: problem !")
726 if self.__Q is None: __Q_shape = (0,0)
727 elif hasattr(self.__Q,"shape"):
728 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
729 else: __Q_shape = self.__Q.shape()
730 else: raise TypeError("The evolution errors covariance matrix (Q) has no attribute of shape: problem !")
732 if len(self.__HO) == 0: __HO_shape = (0,0)
733 elif type(self.__HO) is type({}): __HO_shape = (0,0)
734 elif hasattr(self.__HO["Direct"],"shape"):
735 if type(self.__HO["Direct"].shape) is tuple: __HO_shape = self.__HO["Direct"].shape
736 else: __HO_shape = self.__HO["Direct"].shape()
737 else: raise TypeError("The observation operator (H) has no attribute of shape: problem !")
739 if len(self.__EM) == 0: __EM_shape = (0,0)
740 elif type(self.__EM) is type({}): __EM_shape = (0,0)
741 elif hasattr(self.__EM["Direct"],"shape"):
742 if type(self.__EM["Direct"].shape) is tuple: __EM_shape = self.__EM["Direct"].shape
743 else: __EM_shape = self.__EM["Direct"].shape()
744 else: raise TypeError("The evolution model (EM) has no attribute of shape: problem !")
746 if len(self.__CM) == 0: __CM_shape = (0,0)
747 elif type(self.__CM) is type({}): __CM_shape = (0,0)
748 elif hasattr(self.__CM["Direct"],"shape"):
749 if type(self.__CM["Direct"].shape) is tuple: __CM_shape = self.__CM["Direct"].shape
750 else: __CM_shape = self.__CM["Direct"].shape()
751 else: raise TypeError("The control model (CM) has no attribute of shape: problem !")
753 # Vérification des conditions
754 # ---------------------------
755 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
756 raise ValueError("Shape characteristic of background (Xb) is incorrect: \"%s\"."%(__Xb_shape,))
757 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
758 raise ValueError("Shape characteristic of observation (Y) is incorrect: \"%s\"."%(__Y_shape,))
760 if not( min(__B_shape) == max(__B_shape) ):
761 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) is incorrect: \"%s\"."%(__B_shape,))
762 if not( min(__R_shape) == max(__R_shape) ):
763 raise ValueError("Shape characteristic of observation errors covariance matrix (R) is incorrect: \"%s\"."%(__R_shape,))
764 if not( min(__Q_shape) == max(__Q_shape) ):
765 raise ValueError("Shape characteristic of evolution errors covariance matrix (Q) is incorrect: \"%s\"."%(__Q_shape,))
766 if not( min(__EM_shape) == max(__EM_shape) ):
767 raise ValueError("Shape characteristic of evolution operator (EM) is incorrect: \"%s\"."%(__EM_shape,))
769 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
770 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and state (X) \"%s\" are incompatible."%(__HO_shape,__Xb_shape))
771 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
772 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation (Y) \"%s\" are incompatible."%(__HO_shape,__Y_shape))
773 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] ):
774 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and a priori errors covariance matrix (B) \"%s\" are incompatible."%(__HO_shape,__B_shape))
775 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] ):
776 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation errors covariance matrix (R) \"%s\" are incompatible."%(__HO_shape,__R_shape))
778 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
779 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
780 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
781 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
782 for member in asPersistentVector:
783 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
784 __Xb_shape = min(__B_shape)
786 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) \"%s\" and background (Xb) \"%s\" are incompatible."%(__B_shape,__Xb_shape))
788 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
789 raise ValueError("Shape characteristic of observation errors covariance matrix (R) \"%s\" and observation (Y) \"%s\" are incompatible."%(__R_shape,__Y_shape))
791 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) ):
792 raise ValueError("Shape characteristic of evolution model (EM) \"%s\" and state (X) \"%s\" are incompatible."%(__EM_shape,__Xb_shape))
794 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) ):
795 raise ValueError("Shape characteristic of control model (CM) \"%s\" and control (U) \"%s\" are incompatible."%(__CM_shape,__U_shape))
797 if self.__StoredInputs.has_key("AlgorithmParameters") \
798 and self.__StoredInputs["AlgorithmParameters"].has_key("Bounds") \
799 and (type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type([]) or type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type(())) \
800 and (len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) != max(__Xb_shape)):
801 raise ValueError("The number \"%s\" of bound pairs for the state (X) components is different of the size \"%s\" of the state itself." \
802 %(len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]),max(__Xb_shape)))
806 # -----------------------------------------------------------
809 Permet de lancer le calcul d'assimilation.
811 Le nom de la méthode à activer est toujours "run". Les paramètres en
812 arguments de la méthode sont fixés. En sortie, on obtient les résultats
813 dans la variable de type dictionnaire "StoredVariables", qui contient en
814 particulier des objets de Persistence pour les analyses, OMA...
816 Operator.CM.clearCache()
817 self.shape_validate()
819 self.__algorithm.run(
829 Parameters = self.__Parameters,
833 # -----------------------------------------------------------
834 def get(self, key=None):
836 Renvoie les résultats disponibles après l'exécution de la méthode
837 d'assimilation, ou les diagnostics disponibles. Attention, quand un
838 diagnostic porte le même nom qu'une variable stockée, c'est la variable
839 stockée qui est renvoyée, et le diagnostic est inatteignable.
842 if self.__algorithm.has_key(key):
843 return self.__algorithm.get( key )
844 elif self.__StoredInputs.has_key(key):
845 return self.__StoredInputs[key]
846 elif self.__StoredDiagnostics.has_key(key):
847 return self.__StoredDiagnostics[key]
849 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
851 allvariables = self.__algorithm.get()
852 allvariables.update( self.__StoredDiagnostics )
853 allvariables.update( self.__StoredInputs )
856 def get_available_variables(self):
858 Renvoie les variables potentiellement utilisables pour l'étude,
859 initialement stockées comme données d'entrées ou dans les algorithmes,
860 identifiés par les chaînes de caractères. L'algorithme doit avoir été
861 préalablement choisi sinon la méthode renvoie "None".
863 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
867 if len( self.__algorithm.keys()) > 0:
868 variables.extend( self.__algorithm.get().keys() )
869 if len( self.__StoredInputs.keys() ) > 0:
870 variables.extend( self.__StoredInputs.keys() )
874 def get_available_algorithms(self):
876 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
877 par les chaînes de caractères.
880 for directory in sys.path:
881 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
882 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
883 root, ext = os.path.splitext(fname)
884 if ext == '.py' and root != '__init__':
889 def get_available_diagnostics(self):
891 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
892 par les chaînes de caractères.
895 for directory in sys.path:
896 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
897 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
898 root, ext = os.path.splitext(fname)
899 if ext == '.py' and root != '__init__':
904 # -----------------------------------------------------------
905 def get_algorithms_main_path(self):
907 Renvoie le chemin pour le répertoire principal contenant les algorithmes
908 dans un sous-répertoire "daAlgorithms"
912 def add_algorithms_path(self, asPath=None):
914 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
915 se trouve un sous-répertoire "daAlgorithms"
917 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
918 pas indispensable de le rajouter ici.
920 if not os.path.isdir(asPath):
921 raise ValueError("The given "+asPath+" argument must exist as a directory")
922 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
923 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
924 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
925 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
926 sys.path.insert(0, os.path.abspath(asPath))
927 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
930 def get_diagnostics_main_path(self):
932 Renvoie le chemin pour le répertoire principal contenant les diagnostics
933 dans un sous-répertoire "daDiagnostics"
937 def add_diagnostics_path(self, asPath=None):
939 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
940 se trouve un sous-répertoire "daDiagnostics"
942 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
943 pas indispensable de le rajouter ici.
945 if not os.path.isdir(asPath):
946 raise ValueError("The given "+asPath+" argument must exist as a directory")
947 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
948 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
949 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
950 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
951 sys.path.insert(0, os.path.abspath(asPath))
952 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
955 # -----------------------------------------------------------
956 def setDataObserver(self,
959 HookParameters = None,
963 Permet d'associer un observer à une ou des variables nommées gérées en
964 interne, activable selon des règles définies dans le Scheduler. A chaque
965 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
966 les arguments (variable persistante VariableName, paramètres HookParameters).
969 if type( self.__algorithm ) is dict:
970 raise ValueError("No observer can be build 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 set on a variable named %s which does not exist."%n)
987 self.__algorithm.StoredVariables[ n ].setDataObserver(
988 Scheduler = Scheduler,
989 HookFunction = HookFunction,
990 HookParameters = HookParameters,
993 def removeDataObserver(self,
998 Permet de retirer un observer à une ou des variable nommée.
1001 if type( self.__algorithm ) is dict:
1002 raise ValueError("No observer can be removed before choosing an algorithm.")
1004 # Vérification du nom de variable et typage
1005 # -----------------------------------------
1006 if type( VariableName ) is str:
1007 VariableNames = [VariableName,]
1008 elif type( VariableName ) is list:
1009 VariableNames = map( str, VariableName )
1011 raise ValueError("The observer requires a name or a list of names of variables.")
1013 # Association interne de l'observer à la variable
1014 # -----------------------------------------------
1015 for n in VariableNames:
1016 if not self.__algorithm.has_key( n ):
1017 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
1019 self.__algorithm.StoredVariables[ n ].removeDataObserver(
1020 HookFunction = HookFunction,
1023 # -----------------------------------------------------------
1024 def setDebug(self, level=10):
1026 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
1027 appel pour changer le niveau de verbosité, avec :
1028 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
1030 log = logging.getLogger()
1031 log.setLevel( level )
1033 def unsetDebug(self):
1035 Remet le logger au niveau par défaut
1037 log = logging.getLogger()
1038 log.setLevel( logging.WARNING )
1041 # return set(self.__dict__.keys() + dir(self.__class__))
1042 return ['get', '__doc__', '__init__', '__module__']
1044 def prepare_to_pickle(self):
1046 Retire les variables non pickelisables
1049 del self.__CM # non pickelisable
1050 del self.__EM # non pickelisable
1051 del self.__HO # non pickelisable
1052 del self.__Parameters
1059 del self.__algorithmFile # non pickelisable
1060 del self.__algorithmName
1061 del self.__diagnosticFile # non pickelisable
1062 self.__class__.__doc__ = ""
1064 # ==============================================================================
1065 if __name__ == "__main__":
1066 print '\n AUTODIAGNOSTIC \n'