1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2016 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
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 = PlatformInfo.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 self.__Xb = numpy.matrix( numpy.ravel(numpy.matrix(asVector)), numpy.float ).T
121 elif asPersistentVector is not None:
122 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
123 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
124 for member in asPersistentVector:
125 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
127 self.__Xb = asPersistentVector
129 raise ValueError("Error: improperly defined background, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
131 self.__StoredInputs["Background"] = self.__Xb
134 def setBackgroundError(self,
136 asEyeByScalar = None,
137 asEyeByVector = None,
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 - asCovObject : entrée des données comme un objet ayant des méthodes
153 particulieres de type matriciel
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 self.__B = Covariance(
158 name = "BackgroundError",
159 asCovariance = asCovariance,
160 asEyeByScalar = asEyeByScalar,
161 asEyeByVector = asEyeByVector,
162 asCovObject = asCovObject,
163 toBeChecked = toBeChecked,
166 self.__StoredInputs["BackgroundError"] = self.__B
169 # -----------------------------------------------------------
170 def setObservation(self,
172 asPersistentVector = None,
177 Permet de définir les observations :
178 - asVector : entrée des données, comme un vecteur compatible avec le
179 constructeur de numpy.matrix
180 - asPersistentVector : entrée des données, comme un vecteur de type
181 persistent contruit avec la classe ad-hoc "Persistence"
182 - Scheduler est le contrôle temporel des données disponibles
183 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
184 être rendue disponible au même titre que les variables de calcul
186 if asVector is not None:
187 self.__Y = numpy.matrix( numpy.ravel(numpy.matrix(asVector)), numpy.float ).T
188 elif asPersistentVector is not None:
189 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
190 self.__Y = Persistence.OneVector("Observation", basetype=numpy.matrix)
191 for member in asPersistentVector:
192 self.__Y.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
194 self.__Y = asPersistentVector
196 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
198 self.__StoredInputs["Observation"] = self.__Y
201 def setObservationError(self,
203 asEyeByScalar = None,
204 asEyeByVector = None,
210 Permet de définir la covariance des erreurs d'observations :
211 - asCovariance : entrée des données, comme une matrice compatible avec
212 le constructeur de numpy.matrix
213 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
214 multiplicatif d'une matrice de corrélation identité, aucune matrice
215 n'étant donc explicitement à donner
216 - asEyeByVector : entrée des données comme un seul vecteur de variance,
217 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
218 n'étant donc explicitement à donner
219 - asCovObject : entrée des données comme un objet ayant des méthodes
220 particulieres de type matriciel
221 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
222 être rendue disponible au même titre que les variables de calcul
224 self.__R = Covariance(
225 name = "ObservationError",
226 asCovariance = asCovariance,
227 asEyeByScalar = asEyeByScalar,
228 asEyeByVector = asEyeByVector,
229 asCovObject = asCovObject,
230 toBeChecked = toBeChecked,
233 self.__StoredInputs["ObservationError"] = self.__R
236 def setObservationOperator(self,
244 Permet de définir un opérateur d'observation H. L'ordre de priorité des
245 définitions et leur sens sont les suivants :
246 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
247 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
248 "Direct" n'est pas définie, on prend la fonction "Tangent".
249 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
250 des opérateurs tangents et adjoints. On utilise par défaut des
251 différences finies non centrées ou centrées (si "withCenteredDF" est
252 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
253 du point courant ou sur le point fixe "withdX".
254 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
255 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
256 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
257 matrice fournie doit être sous une forme compatible avec le
258 constructeur de numpy.matrix.
259 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
260 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
261 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
262 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
263 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
264 être rendue disponible au même titre que les variables de calcul
265 L'argument "asFunction" peut prendre la forme complète suivante, avec
266 les valeurs par défaut standards :
267 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
268 "useApproximatedDerivatives":False,
269 "withCenteredDF" :False,
270 "withIncrement" :0.01,
272 "withAvoidingRedundancy" :True,
273 "withToleranceInRedundancy" :1.e-18,
274 "withLenghtOfRedundancy" :-1,
275 "withmpEnabled" :False,
276 "withmpWorkers" :None,
279 if isinstance(asFunction, dict) and \
280 ("useApproximatedDerivatives" in asFunction) and bool(asFunction["useApproximatedDerivatives"]) and \
281 ("Direct" in asFunction) and (asFunction["Direct"] is not None):
282 if "withCenteredDF" not in asFunction: asFunction["withCenteredDF"] = False
283 if "withIncrement" not in asFunction: asFunction["withIncrement"] = 0.01
284 if "withdX" not in asFunction: asFunction["withdX"] = None
285 if "withAvoidingRedundancy" not in asFunction: asFunction["withAvoidingRedundancy"] = True
286 if "withToleranceInRedundancy" not in asFunction: asFunction["withToleranceInRedundancy"] = 1.e-18
287 if "withLenghtOfRedundancy" not in asFunction: asFunction["withLenghtOfRedundancy"] = -1
288 if "withmpEnabled" not in asFunction: asFunction["withmpEnabled"] = False
289 if "withmpWorkers" not in asFunction: asFunction["withmpWorkers"] = None
290 from daNumerics.ApproximatedDerivatives import FDApproximation
291 FDA = FDApproximation(
292 Function = asFunction["Direct"],
293 centeredDF = asFunction["withCenteredDF"],
294 increment = asFunction["withIncrement"],
295 dX = asFunction["withdX"],
296 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
297 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
298 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
299 mpEnabled = asFunction["withmpEnabled"],
300 mpWorkers = asFunction["withmpWorkers"],
302 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
303 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
304 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
305 elif isinstance(asFunction, dict) and \
306 ("Tangent" in asFunction) and ("Adjoint" in asFunction) and \
307 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
308 if ("Direct" not in asFunction) or (asFunction["Direct"] is None):
309 self.__HO["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
311 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
312 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
313 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
314 elif asMatrix is not None:
315 matrice = numpy.matrix( asMatrix, numpy.float )
316 self.__HO["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
317 self.__HO["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
318 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
321 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
323 if appliedToX is not None:
324 self.__HO["AppliedToX"] = {}
325 if type(appliedToX) is not dict:
326 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
327 for key in appliedToX.keys():
328 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
329 # Pour le cas où l'on a une vraie matrice
330 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
331 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
332 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
333 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
335 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
337 self.__HO["AppliedToX"] = None
340 self.__StoredInputs["ObservationOperator"] = self.__HO
343 # -----------------------------------------------------------
344 def setEvolutionModel(self,
352 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
353 définitions et leur sens sont les suivants :
354 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
355 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
356 "Direct" n'est pas définie, on prend la fonction "Tangent".
357 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
358 des opérateurs tangents et adjoints. On utilise par défaut des
359 différences finies non centrées ou centrées (si "withCenteredDF" est
360 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
361 du point courant ou sur le point fixe "withdX".
362 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
363 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
364 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
365 matrice fournie doit être sous une forme compatible avec le
366 constructeur de numpy.matrix.
367 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
368 être rendue disponible au même titre que les variables de calcul
369 L'argument "asFunction" peut prendre la forme complète suivante, avec
370 les valeurs par défaut standards :
371 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
372 "useApproximatedDerivatives":False,
373 "withCenteredDF" :False,
374 "withIncrement" :0.01,
376 "withAvoidingRedundancy" :True,
377 "withToleranceInRedundancy" :1.e-18,
378 "withLenghtOfRedundancy" :-1,
379 "withmpEnabled" :False,
380 "withmpWorkers" :None,
383 if isinstance(asFunction, dict) and \
384 ("useApproximatedDerivatives" in asFunction) and bool(asFunction["useApproximatedDerivatives"]) and \
385 ("Direct" in asFunction) and (asFunction["Direct"] is not None):
386 if "withCenteredDF" not in asFunction: asFunction["withCenteredDF"] = False
387 if "withIncrement" not in asFunction: asFunction["withIncrement"] = 0.01
388 if "withdX" not in asFunction: asFunction["withdX"] = None
389 if "withAvoidingRedundancy" not in asFunction: asFunction["withAvoidingRedundancy"] = True
390 if "withToleranceInRedundancy" not in asFunction: asFunction["withToleranceInRedundancy"] = 1.e-18
391 if "withLenghtOfRedundancy" not in asFunction: asFunction["withLenghtOfRedundancy"] = -1
392 if "withmpEnabled" not in asFunction: asFunction["withmpEnabled"] = False
393 if "withmpWorkers" not in asFunction: asFunction["withmpWorkers"] = None
394 from daNumerics.ApproximatedDerivatives import FDApproximation
395 FDA = FDApproximation(
396 Function = asFunction["Direct"],
397 centeredDF = asFunction["withCenteredDF"],
398 increment = asFunction["withIncrement"],
399 dX = asFunction["withdX"],
400 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
401 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
402 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
403 mpEnabled = asFunction["withmpEnabled"],
404 mpWorkers = asFunction["withmpWorkers"],
406 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
407 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
408 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
409 elif isinstance(asFunction, dict) and \
410 ("Tangent" in asFunction) and ("Adjoint" in asFunction) and \
411 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
412 if ("Direct" not in asFunction) or (asFunction["Direct"] is None):
413 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
415 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
416 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
417 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
418 elif asMatrix is not None:
419 matrice = numpy.matrix( asMatrix, numpy.float )
420 self.__EM["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
421 self.__EM["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
422 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
425 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
428 self.__StoredInputs["EvolutionModel"] = self.__EM
431 def setEvolutionError(self,
433 asEyeByScalar = None,
434 asEyeByVector = None,
440 Permet de définir la covariance des erreurs de modèle :
441 - asCovariance : entrée des données, comme une matrice compatible avec
442 le constructeur de numpy.matrix
443 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
444 multiplicatif d'une matrice de corrélation identité, aucune matrice
445 n'étant donc explicitement à donner
446 - asEyeByVector : entrée des données comme un seul vecteur de variance,
447 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
448 n'étant donc explicitement à donner
449 - asCovObject : entrée des données comme un objet ayant des méthodes
450 particulieres de type matriciel
451 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
452 être rendue disponible au même titre que les variables de calcul
454 self.__Q = Covariance(
455 name = "EvolutionError",
456 asCovariance = asCovariance,
457 asEyeByScalar = asEyeByScalar,
458 asEyeByVector = asEyeByVector,
459 asCovObject = asCovObject,
460 toBeChecked = toBeChecked,
463 self.__StoredInputs["EvolutionError"] = self.__Q
466 # -----------------------------------------------------------
467 def setControlModel(self,
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
492 L'argument "asFunction" peut prendre la forme complète suivante, avec
493 les valeurs par défaut standards :
494 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
495 "useApproximatedDerivatives":False,
496 "withCenteredDF" :False,
497 "withIncrement" :0.01,
499 "withAvoidingRedundancy" :True,
500 "withToleranceInRedundancy" :1.e-18,
501 "withLenghtOfRedundancy" :-1,
502 "withmpEnabled" :False,
503 "withmpWorkers" :None,
506 if isinstance(asFunction, dict) and \
507 ("useApproximatedDerivatives" in asFunction) and bool(asFunction["useApproximatedDerivatives"]) and \
508 ("Direct" in asFunction) and (asFunction["Direct"] is not None):
509 if "withCenteredDF" not in asFunction: asFunction["withCenteredDF"] = False
510 if "withIncrement" not in asFunction: asFunction["withIncrement"] = 0.01
511 if "withdX" not in asFunction: asFunction["withdX"] = None
512 if "withAvoidingRedundancy" not in asFunction: asFunction["withAvoidingRedundancy"] = True
513 if "withToleranceInRedundancy" not in asFunction: asFunction["withToleranceInRedundancy"] = 1.e-18
514 if "withLenghtOfRedundancy" not in asFunction: asFunction["withLenghtOfRedundancy"] = -1
515 if "withmpEnabled" not in asFunction: asFunction["withmpEnabled"] = False
516 if "withmpWorkers" not in asFunction: asFunction["withmpWorkers"] = None
517 from daNumerics.ApproximatedDerivatives import FDApproximation
518 FDA = FDApproximation(
519 Function = asFunction["Direct"],
520 centeredDF = asFunction["withCenteredDF"],
521 increment = asFunction["withIncrement"],
522 dX = asFunction["withdX"],
523 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
524 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
525 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
526 mpEnabled = asFunction["withmpEnabled"],
527 mpWorkers = asFunction["withmpWorkers"],
529 self.__CM["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
530 self.__CM["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
531 self.__CM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
532 elif isinstance(asFunction, dict) and \
533 ("Tangent" in asFunction) and ("Adjoint" in asFunction) and \
534 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
535 if ("Direct" not in asFunction) or (asFunction["Direct"] is None):
536 self.__CM["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
538 self.__CM["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
539 self.__CM["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
540 self.__CM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
541 elif asMatrix is not None:
542 matrice = numpy.matrix( asMatrix, numpy.float )
543 self.__CM["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
544 self.__CM["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
545 self.__CM["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
548 raise ValueError("Improperly defined input control model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
551 self.__StoredInputs["ControlModel"] = self.__CM
554 def setControlInput(self,
556 asPersistentVector = None,
561 Permet de définir le controle en entree :
562 - asVector : entrée des données, comme un vecteur compatible avec le
563 constructeur de numpy.matrix
564 - asPersistentVector : entrée des données, comme un vecteur de type
565 persistent contruit avec la classe ad-hoc "Persistence"
566 - Scheduler est le contrôle temporel des données disponibles
567 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
568 être rendue disponible au même titre que les variables de calcul
570 if asVector is not None:
571 self.__U = numpy.matrix( numpy.ravel(numpy.matrix(asVector)), numpy.float ).T
572 elif asPersistentVector is not None:
573 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
574 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.matrix)
575 for member in asPersistentVector:
576 self.__U.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
578 self.__U = asPersistentVector
580 raise ValueError("Error: improperly defined control input, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
582 self.__StoredInputs["ControlInput"] = self.__U
585 # -----------------------------------------------------------
586 def setControls (self,
591 Permet de définir la valeur initiale du vecteur X contenant toutes les
592 variables de contrôle, i.e. les paramètres ou l'état dont on veut
593 estimer la valeur pour obtenir les observations. C'est utile pour un
594 algorithme itératif/incrémental.
595 - asVector : entrée des données, comme un vecteur compatible avec le
596 constructeur de numpy.matrix.
597 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
598 être rendue disponible au même titre que les variables de calcul
600 if asVector is not None:
601 self.__X.store( asVector )
603 self.__StoredInputs["Controls"] = self.__X
606 # -----------------------------------------------------------
607 def setAlgorithm(self, choice = None ):
609 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
610 d'assimilation. L'argument est un champ caractère se rapportant au nom
611 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
612 d'assimilation sur les arguments fixes.
615 raise ValueError("Error: algorithm choice has to be given")
616 if self.__algorithmName is not None:
617 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
618 daDirectory = "daAlgorithms"
620 # Recherche explicitement le fichier complet
621 # ------------------------------------------
623 for directory in sys.path:
624 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
625 module_path = os.path.abspath(os.path.join(directory, daDirectory))
626 if module_path is None:
627 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
629 # Importe le fichier complet comme un module
630 # ------------------------------------------
632 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
633 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
634 self.__algorithmName = str(choice)
635 sys.path = sys_path_tmp ; del sys_path_tmp
636 except ImportError, e:
637 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
639 # Instancie un objet du type élémentaire du fichier
640 # -------------------------------------------------
641 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
642 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
645 def setAlgorithmParameters(self, asDico=None):
647 Permet de définir les paramètres de l'algorithme, sous la forme d'un
650 if asDico is not None:
651 self.__Parameters.update( dict( asDico ) )
653 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
656 def getAlgorithmParameters(self, noDetails=True):
658 Renvoie la liste des paramètres requis selon l'algorithme
660 return self.__algorithm.getRequiredParameters(noDetails)
662 # -----------------------------------------------------------
663 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
665 Permet de sélectionner un diagnostic a effectuer.
668 raise ValueError("Error: diagnostic choice has to be given")
669 daDirectory = "daDiagnostics"
671 # Recherche explicitement le fichier complet
672 # ------------------------------------------
674 for directory in sys.path:
675 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
676 module_path = os.path.abspath(os.path.join(directory, daDirectory))
677 if module_path is None:
678 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
680 # Importe le fichier complet comme un module
681 # ------------------------------------------
683 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
684 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
685 sys.path = sys_path_tmp ; del sys_path_tmp
686 except ImportError, e:
687 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
689 # Instancie un objet du type élémentaire du fichier
690 # -------------------------------------------------
691 if name in self.__StoredInputs:
692 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
693 elif name in self.__StoredDiagnostics:
694 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
696 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
700 parameters = parameters )
703 # -----------------------------------------------------------
704 def shape_validate(self):
706 Validation de la correspondance correcte des tailles des variables et
707 des matrices s'il y en a.
709 if self.__Xb is None: __Xb_shape = (0,)
710 elif hasattr(self.__Xb,"size"): __Xb_shape = (self.__Xb.size,)
711 elif hasattr(self.__Xb,"shape"):
712 if isinstance(self.__Xb.shape, tuple): __Xb_shape = self.__Xb.shape
713 else: __Xb_shape = self.__Xb.shape()
714 else: raise TypeError("The background (Xb) has no attribute of shape: problem !")
716 if self.__Y is None: __Y_shape = (0,)
717 elif hasattr(self.__Y,"size"): __Y_shape = (self.__Y.size,)
718 elif hasattr(self.__Y,"shape"):
719 if isinstance(self.__Y.shape, tuple): __Y_shape = self.__Y.shape
720 else: __Y_shape = self.__Y.shape()
721 else: raise TypeError("The observation (Y) has no attribute of shape: problem !")
723 if self.__U is None: __U_shape = (0,)
724 elif hasattr(self.__U,"size"): __U_shape = (self.__U.size,)
725 elif hasattr(self.__U,"shape"):
726 if isinstance(self.__U.shape, tuple): __U_shape = self.__U.shape
727 else: __U_shape = self.__U.shape()
728 else: raise TypeError("The control (U) has no attribute of shape: problem !")
730 if self.__B is None: __B_shape = (0,0)
731 elif hasattr(self.__B,"shape"):
732 if isinstance(self.__B.shape, tuple): __B_shape = self.__B.shape
733 else: __B_shape = self.__B.shape()
734 else: raise TypeError("The a priori errors covariance matrix (B) has no attribute of shape: problem !")
736 if self.__R is None: __R_shape = (0,0)
737 elif hasattr(self.__R,"shape"):
738 if isinstance(self.__R.shape, tuple): __R_shape = self.__R.shape
739 else: __R_shape = self.__R.shape()
740 else: raise TypeError("The observation errors covariance matrix (R) has no attribute of shape: problem !")
742 if self.__Q is None: __Q_shape = (0,0)
743 elif hasattr(self.__Q,"shape"):
744 if isinstance(self.__Q.shape, tuple): __Q_shape = self.__Q.shape
745 else: __Q_shape = self.__Q.shape()
746 else: raise TypeError("The evolution errors covariance matrix (Q) has no attribute of shape: problem !")
748 if len(self.__HO) == 0: __HO_shape = (0,0)
749 elif isinstance(self.__HO, dict): __HO_shape = (0,0)
750 elif hasattr(self.__HO["Direct"],"shape"):
751 if isinstance(self.__HO["Direct"].shape, tuple): __HO_shape = self.__HO["Direct"].shape
752 else: __HO_shape = self.__HO["Direct"].shape()
753 else: raise TypeError("The observation operator (H) has no attribute of shape: problem !")
755 if len(self.__EM) == 0: __EM_shape = (0,0)
756 elif isinstance(self.__EM, dict): __EM_shape = (0,0)
757 elif hasattr(self.__EM["Direct"],"shape"):
758 if isinstance(self.__EM["Direct"].shape, tuple): __EM_shape = self.__EM["Direct"].shape
759 else: __EM_shape = self.__EM["Direct"].shape()
760 else: raise TypeError("The evolution model (EM) has no attribute of shape: problem !")
762 if len(self.__CM) == 0: __CM_shape = (0,0)
763 elif isinstance(self.__CM, dict): __CM_shape = (0,0)
764 elif hasattr(self.__CM["Direct"],"shape"):
765 if isinstance(self.__CM["Direct"].shape, tuple): __CM_shape = self.__CM["Direct"].shape
766 else: __CM_shape = self.__CM["Direct"].shape()
767 else: raise TypeError("The control model (CM) has no attribute of shape: problem !")
769 # Vérification des conditions
770 # ---------------------------
771 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
772 raise ValueError("Shape characteristic of background (Xb) is incorrect: \"%s\"."%(__Xb_shape,))
773 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
774 raise ValueError("Shape characteristic of observation (Y) is incorrect: \"%s\"."%(__Y_shape,))
776 if not( min(__B_shape) == max(__B_shape) ):
777 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) is incorrect: \"%s\"."%(__B_shape,))
778 if not( min(__R_shape) == max(__R_shape) ):
779 raise ValueError("Shape characteristic of observation errors covariance matrix (R) is incorrect: \"%s\"."%(__R_shape,))
780 if not( min(__Q_shape) == max(__Q_shape) ):
781 raise ValueError("Shape characteristic of evolution errors covariance matrix (Q) is incorrect: \"%s\"."%(__Q_shape,))
782 if not( min(__EM_shape) == max(__EM_shape) ):
783 raise ValueError("Shape characteristic of evolution operator (EM) is incorrect: \"%s\"."%(__EM_shape,))
785 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
786 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and state (X) \"%s\" are incompatible."%(__HO_shape,__Xb_shape))
787 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
788 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation (Y) \"%s\" are incompatible."%(__HO_shape,__Y_shape))
789 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] ):
790 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and a priori errors covariance matrix (B) \"%s\" are incompatible."%(__HO_shape,__B_shape))
791 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] ):
792 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation errors covariance matrix (R) \"%s\" are incompatible."%(__HO_shape,__R_shape))
794 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
795 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
796 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
797 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
798 for member in asPersistentVector:
799 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
800 __Xb_shape = min(__B_shape)
802 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) \"%s\" and background (Xb) \"%s\" are incompatible."%(__B_shape,__Xb_shape))
804 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
805 raise ValueError("Shape characteristic of observation errors covariance matrix (R) \"%s\" and observation (Y) \"%s\" are incompatible."%(__R_shape,__Y_shape))
807 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) ):
808 raise ValueError("Shape characteristic of evolution model (EM) \"%s\" and state (X) \"%s\" are incompatible."%(__EM_shape,__Xb_shape))
810 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) ):
811 raise ValueError("Shape characteristic of control model (CM) \"%s\" and control (U) \"%s\" are incompatible."%(__CM_shape,__U_shape))
813 if ("AlgorithmParameters" in self.__StoredInputs) \
814 and ("Bounds" in self.__StoredInputs["AlgorithmParameters"]) \
815 and (isinstance(self.__StoredInputs["AlgorithmParameters"]["Bounds"], list) or isinstance(self.__StoredInputs["AlgorithmParameters"]["Bounds"], tuple)) \
816 and (len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) != max(__Xb_shape)):
817 raise ValueError("The number \"%s\" of bound pairs for the state (X) components is different of the size \"%s\" of the state itself." \
818 %(len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]),max(__Xb_shape)))
822 # -----------------------------------------------------------
825 Permet de lancer le calcul d'assimilation.
827 Le nom de la méthode à activer est toujours "run". Les paramètres en
828 arguments de la méthode sont fixés. En sortie, on obtient les résultats
829 dans la variable de type dictionnaire "StoredVariables", qui contient en
830 particulier des objets de Persistence pour les analyses, OMA...
832 Operator.CM.clearCache()
833 self.shape_validate()
835 self.__algorithm.run(
845 Parameters = self.__Parameters,
849 # -----------------------------------------------------------
850 def get(self, key=None):
852 Renvoie les résultats disponibles après l'exécution de la méthode
853 d'assimilation, ou les diagnostics disponibles. Attention, quand un
854 diagnostic porte le même nom qu'une variable stockée, c'est la variable
855 stockée qui est renvoyée, et le diagnostic est inatteignable.
858 if key in self.__algorithm:
859 return self.__algorithm.get( key )
860 elif key in self.__StoredInputs:
861 return self.__StoredInputs[key]
862 elif key in self.__StoredDiagnostics:
863 return self.__StoredDiagnostics[key]
865 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
867 allvariables = self.__algorithm.get()
868 allvariables.update( self.__StoredDiagnostics )
869 allvariables.update( self.__StoredInputs )
872 def get_available_variables(self):
874 Renvoie les variables potentiellement utilisables pour l'étude,
875 initialement stockées comme données d'entrées ou dans les algorithmes,
876 identifiés par les chaînes de caractères. L'algorithme doit avoir été
877 préalablement choisi sinon la méthode renvoie "None".
879 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
883 if len( self.__algorithm.keys()) > 0:
884 variables.extend( self.__algorithm.get().keys() )
885 if len( self.__StoredInputs.keys() ) > 0:
886 variables.extend( self.__StoredInputs.keys() )
890 def get_available_algorithms(self):
892 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
893 par les chaînes de caractères.
896 for directory in sys.path:
897 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
898 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
899 root, ext = os.path.splitext(fname)
900 if ext == '.py' and root != '__init__':
905 def get_available_diagnostics(self):
907 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
908 par les chaînes de caractères.
911 for directory in sys.path:
912 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
913 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
914 root, ext = os.path.splitext(fname)
915 if ext == '.py' and root != '__init__':
920 # -----------------------------------------------------------
921 def get_algorithms_main_path(self):
923 Renvoie le chemin pour le répertoire principal contenant les algorithmes
924 dans un sous-répertoire "daAlgorithms"
928 def add_algorithms_path(self, asPath=None):
930 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
931 se trouve un sous-répertoire "daAlgorithms"
933 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
934 pas indispensable de le rajouter ici.
936 if not os.path.isdir(asPath):
937 raise ValueError("The given "+asPath+" argument must exist as a directory")
938 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
939 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
940 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
941 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
942 sys.path.insert(0, os.path.abspath(asPath))
943 sys.path = PlatformInfo.uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
946 def get_diagnostics_main_path(self):
948 Renvoie le chemin pour le répertoire principal contenant les diagnostics
949 dans un sous-répertoire "daDiagnostics"
953 def add_diagnostics_path(self, asPath=None):
955 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
956 se trouve un sous-répertoire "daDiagnostics"
958 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
959 pas indispensable de le rajouter ici.
961 if not os.path.isdir(asPath):
962 raise ValueError("The given "+asPath+" argument must exist as a directory")
963 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
964 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
965 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
966 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
967 sys.path.insert(0, os.path.abspath(asPath))
968 sys.path = PlatformInfo.uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
971 # -----------------------------------------------------------
972 def setDataObserver(self,
975 HookParameters = None,
979 Permet d'associer un observer à une ou des variables nommées gérées en
980 interne, activable selon des règles définies dans le Scheduler. A chaque
981 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
982 les arguments (variable persistante VariableName, paramètres HookParameters).
985 if isinstance(self.__algorithm, dict):
986 raise ValueError("No observer can be build before choosing an algorithm.")
988 # Vérification du nom de variable et typage
989 # -----------------------------------------
990 if isinstance(VariableName, str):
991 VariableNames = [VariableName,]
992 elif isinstance(VariableName, 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 n not in self.__algorithm:
1001 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
1003 self.__algorithm.StoredVariables[ n ].setDataObserver(
1004 Scheduler = Scheduler,
1005 HookFunction = HookFunction,
1006 HookParameters = HookParameters,
1009 def removeDataObserver(self,
1010 VariableName = None,
1011 HookFunction = None,
1014 Permet de retirer un observer à une ou des variable nommée.
1017 if isinstance(self.__algorithm, dict):
1018 raise ValueError("No observer can be removed before choosing an algorithm.")
1020 # Vérification du nom de variable et typage
1021 # -----------------------------------------
1022 if isinstance(VariableName, str):
1023 VariableNames = [VariableName,]
1024 elif isinstance(VariableName, list):
1025 VariableNames = map( str, VariableName )
1027 raise ValueError("The observer requires a name or a list of names of variables.")
1029 # Association interne de l'observer à la variable
1030 # -----------------------------------------------
1031 for n in VariableNames:
1032 if n not in self.__algorithm:
1033 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
1035 self.__algorithm.StoredVariables[ n ].removeDataObserver(
1036 HookFunction = HookFunction,
1039 # -----------------------------------------------------------
1040 def setDebug(self, level=10):
1042 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
1043 appel pour changer le niveau de verbosité, avec :
1044 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
1046 log = logging.getLogger()
1047 log.setLevel( level )
1049 def unsetDebug(self):
1051 Remet le logger au niveau par défaut
1053 log = logging.getLogger()
1054 log.setLevel( logging.WARNING )
1057 # return set(self.__dict__.keys() + dir(self.__class__))
1058 return ['get', '__doc__', '__init__', '__module__']
1060 def prepare_to_pickle(self):
1062 Retire les variables non pickelisables
1065 del self.__CM # non pickelisable
1066 del self.__EM # non pickelisable
1067 del self.__HO # non pickelisable
1068 del self.__Parameters
1075 del self.__algorithmFile # non pickelisable
1076 del self.__algorithmName
1077 del self.__diagnosticFile # non pickelisable
1078 self.__class__.__doc__ = ""
1080 # ==============================================================================
1081 if __name__ == "__main__":
1082 print '\n AUTODIAGNOSTIC \n'