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 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,
142 Permet de définir la covariance des erreurs d'ébauche :
143 - asCovariance : entrée des données, comme une matrice compatible avec
144 le constructeur de numpy.matrix
145 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
146 multiplicatif d'une matrice de corrélation identité, aucune matrice
147 n'étant donc explicitement à donner
148 - asEyeByVector : entrée des données comme un seul vecteur de variance,
149 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
150 n'étant donc explicitement à donner
151 - asCovObject : entrée des données comme un objet ayant des méthodes
152 particulieres de type matriciel
153 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
154 être rendue disponible au même titre que les variables de calcul
156 self.__B = Covariance(
157 name = "BackgroundError",
158 asCovariance = asCovariance,
159 asEyeByScalar = asEyeByScalar,
160 asEyeByVector = asEyeByVector,
161 asCovObject = asCovObject,
164 self.__StoredInputs["BackgroundError"] = self.__B
167 # -----------------------------------------------------------
168 def setObservation(self,
170 asPersistentVector = None,
175 Permet de définir les observations :
176 - asVector : entrée des données, comme un vecteur compatible avec le
177 constructeur de numpy.matrix
178 - asPersistentVector : entrée des données, comme un vecteur de type
179 persistent contruit avec la classe ad-hoc "Persistence"
180 - Scheduler est le contrôle temporel des données disponibles
181 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
182 être rendue disponible au même titre que les variables de calcul
184 if asVector is not None:
185 self.__Y = numpy.matrix( numpy.ravel(numpy.matrix(asVector)), numpy.float ).T
186 elif asPersistentVector is not None:
187 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
188 self.__Y = Persistence.OneVector("Observation", basetype=numpy.matrix)
189 for member in asPersistentVector:
190 self.__Y.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
192 self.__Y = asPersistentVector
194 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
196 self.__StoredInputs["Observation"] = self.__Y
199 def setObservationError(self,
201 asEyeByScalar = None,
202 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 - asCovObject : entrée des données comme un objet ayant des méthodes
217 particulieres de type matriciel
218 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
219 être rendue disponible au même titre que les variables de calcul
221 self.__R = Covariance(
222 name = "ObservationError",
223 asCovariance = asCovariance,
224 asEyeByScalar = asEyeByScalar,
225 asEyeByVector = asEyeByVector,
226 asCovObject = asCovObject,
229 self.__StoredInputs["ObservationError"] = self.__R
232 def setObservationOperator(self,
240 Permet de définir un opérateur d'observation H. L'ordre de priorité des
241 définitions et leur sens sont les suivants :
242 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
243 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
244 "Direct" n'est pas définie, on prend la fonction "Tangent".
245 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
246 des opérateurs tangents et adjoints. On utilise par défaut des
247 différences finies non centrées ou centrées (si "withCenteredDF" est
248 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
249 du point courant ou sur le point fixe "withdX".
250 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
251 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
252 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
253 matrice fournie doit être sous une forme compatible avec le
254 constructeur de numpy.matrix.
255 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
256 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
257 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
258 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
259 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
260 être rendue disponible au même titre que les variables de calcul
261 L'argument "asFunction" peut prendre la forme complète suivante, avec
262 les valeurs par défaut standards :
263 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
264 "useApproximatedDerivatives":False,
265 "withCenteredDF" :False,
266 "withIncrement" :0.01,
268 "withAvoidingRedundancy" :True,
269 "withToleranceInRedundancy" :1.e-18,
270 "withLenghtOfRedundancy" :-1,
271 "withmpEnabled" :False,
272 "withmpWorkers" :None,
275 if (type(asFunction) is type({})) and \
276 ("useApproximatedDerivatives" in asFunction) and bool(asFunction["useApproximatedDerivatives"]) and \
277 ("Direct" in asFunction) and (asFunction["Direct"] is not None):
278 if "withCenteredDF" not in asFunction: asFunction["withCenteredDF"] = False
279 if "withIncrement" not in asFunction: asFunction["withIncrement"] = 0.01
280 if "withdX" not in asFunction: asFunction["withdX"] = None
281 if "withAvoidingRedundancy" not in asFunction: asFunction["withAvoidingRedundancy"] = True
282 if "withToleranceInRedundancy" not in asFunction: asFunction["withToleranceInRedundancy"] = 1.e-18
283 if "withLenghtOfRedundancy" not in asFunction: asFunction["withLenghtOfRedundancy"] = -1
284 if "withmpEnabled" not in asFunction: asFunction["withmpEnabled"] = False
285 if "withmpWorkers" not in asFunction: asFunction["withmpWorkers"] = None
286 from daNumerics.ApproximatedDerivatives import FDApproximation
287 FDA = FDApproximation(
288 Function = asFunction["Direct"],
289 centeredDF = asFunction["withCenteredDF"],
290 increment = asFunction["withIncrement"],
291 dX = asFunction["withdX"],
292 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
293 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
294 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
295 mpEnabled = asFunction["withmpEnabled"],
296 mpWorkers = asFunction["withmpWorkers"],
298 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
299 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
300 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
301 elif (type(asFunction) is type({})) and \
302 ("Tangent" in asFunction) and ("Adjoint" in asFunction) and \
303 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
304 if ("Direct" not in asFunction) or (asFunction["Direct"] is None):
305 self.__HO["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
307 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
308 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
309 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
310 elif asMatrix is not None:
311 matrice = numpy.matrix( asMatrix, numpy.float )
312 self.__HO["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
313 self.__HO["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
314 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
317 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
319 if appliedToX is not None:
320 self.__HO["AppliedToX"] = {}
321 if type(appliedToX) is not dict:
322 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
323 for key in appliedToX.keys():
324 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
325 # Pour le cas où l'on a une vraie matrice
326 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
327 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
328 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
329 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
331 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
333 self.__HO["AppliedToX"] = None
336 self.__StoredInputs["ObservationOperator"] = self.__HO
339 # -----------------------------------------------------------
340 def setEvolutionModel(self,
348 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
349 définitions et leur sens sont les suivants :
350 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
351 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
352 "Direct" n'est pas définie, on prend la fonction "Tangent".
353 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
354 des opérateurs tangents et adjoints. On utilise par défaut des
355 différences finies non centrées ou centrées (si "withCenteredDF" est
356 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
357 du point courant ou sur le point fixe "withdX".
358 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
359 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
360 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
361 matrice fournie doit être sous une forme compatible avec le
362 constructeur de numpy.matrix.
363 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
364 être rendue disponible au même titre que les variables de calcul
365 L'argument "asFunction" peut prendre la forme complète suivante, avec
366 les valeurs par défaut standards :
367 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
368 "useApproximatedDerivatives":False,
369 "withCenteredDF" :False,
370 "withIncrement" :0.01,
372 "withAvoidingRedundancy" :True,
373 "withToleranceInRedundancy" :1.e-18,
374 "withLenghtOfRedundancy" :-1,
375 "withmpEnabled" :False,
376 "withmpWorkers" :None,
379 if (type(asFunction) is type({})) and \
380 ("useApproximatedDerivatives" in asFunction) and bool(asFunction["useApproximatedDerivatives"]) and \
381 ("Direct" in asFunction) and (asFunction["Direct"] is not None):
382 if "withCenteredDF" not in asFunction: asFunction["withCenteredDF"] = False
383 if "withIncrement" not in asFunction: asFunction["withIncrement"] = 0.01
384 if "withdX" not in asFunction: asFunction["withdX"] = None
385 if "withAvoidingRedundancy" not in asFunction: asFunction["withAvoidingRedundancy"] = True
386 if "withToleranceInRedundancy" not in asFunction: asFunction["withToleranceInRedundancy"] = 1.e-18
387 if "withLenghtOfRedundancy" not in asFunction: asFunction["withLenghtOfRedundancy"] = -1
388 if "withmpEnabled" not in asFunction: asFunction["withmpEnabled"] = False
389 if "withmpWorkers" not in asFunction: asFunction["withmpWorkers"] = None
390 from daNumerics.ApproximatedDerivatives import FDApproximation
391 FDA = FDApproximation(
392 Function = asFunction["Direct"],
393 centeredDF = asFunction["withCenteredDF"],
394 increment = asFunction["withIncrement"],
395 dX = asFunction["withdX"],
396 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
397 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
398 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
399 mpEnabled = asFunction["withmpEnabled"],
400 mpWorkers = asFunction["withmpWorkers"],
402 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
403 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
404 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
405 elif (type(asFunction) is type({})) and \
406 ("Tangent" in asFunction) and ("Adjoint" in asFunction) and \
407 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
408 if ("Direct" not in asFunction) or (asFunction["Direct"] is None):
409 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
411 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
412 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
413 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
414 elif asMatrix is not None:
415 matrice = numpy.matrix( asMatrix, numpy.float )
416 self.__EM["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
417 self.__EM["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
418 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
421 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
424 self.__StoredInputs["EvolutionModel"] = self.__EM
427 def setEvolutionError(self,
429 asEyeByScalar = None,
430 asEyeByVector = None,
435 Permet de définir la covariance des erreurs de modèle :
436 - asCovariance : entrée des données, comme une matrice compatible avec
437 le constructeur de numpy.matrix
438 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
439 multiplicatif d'une matrice de corrélation identité, aucune matrice
440 n'étant donc explicitement à donner
441 - asEyeByVector : entrée des données comme un seul vecteur de variance,
442 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
443 n'étant donc explicitement à donner
444 - asCovObject : entrée des données comme un objet ayant des méthodes
445 particulieres de type matriciel
446 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
447 être rendue disponible au même titre que les variables de calcul
449 self.__Q = Covariance(
450 name = "EvolutionError",
451 asCovariance = asCovariance,
452 asEyeByScalar = asEyeByScalar,
453 asEyeByVector = asEyeByVector,
454 asCovObject = asCovObject,
457 self.__StoredInputs["EvolutionError"] = self.__Q
460 # -----------------------------------------------------------
461 def setControlModel(self,
462 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
463 "useApproximatedDerivatives":False,
464 "withCenteredDF" :False,
465 "withIncrement" :0.01,
474 Permet de définir un opérateur de controle C. L'ordre de priorité des
475 définitions et leur sens sont les suivants :
476 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
477 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
478 "Direct" n'est pas définie, on prend la fonction "Tangent".
479 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
480 des opérateurs tangents et adjoints. On utilise par défaut des
481 différences finies non centrées ou centrées (si "withCenteredDF" est
482 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
483 du point courant ou sur le point fixe "withdX".
484 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
485 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
486 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
487 matrice fournie doit être sous une forme compatible avec le
488 constructeur de numpy.matrix.
489 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
490 être rendue disponible au même titre que les variables de calcul
491 L'argument "asFunction" peut prendre la forme complète suivante, avec
492 les valeurs par défaut standards :
493 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
494 "useApproximatedDerivatives":False,
495 "withCenteredDF" :False,
496 "withIncrement" :0.01,
498 "withAvoidingRedundancy" :True,
499 "withToleranceInRedundancy" :1.e-18,
500 "withLenghtOfRedundancy" :-1,
501 "withmpEnabled" :False,
502 "withmpWorkers" :None,
505 if (type(asFunction) is type({})) and \
506 ("useApproximatedDerivatives" in asFunction) and bool(asFunction["useApproximatedDerivatives"]) and \
507 ("Direct" in asFunction) and (asFunction["Direct"] is not None):
508 if "withCenteredDF" not in asFunction: asFunction["withCenteredDF"] = False
509 if "withIncrement" not in asFunction: asFunction["withIncrement"] = 0.01
510 if "withdX" not in asFunction: asFunction["withdX"] = None
511 if "withAvoidingRedundancy" not in asFunction: asFunction["withAvoidingRedundancy"] = True
512 if "withToleranceInRedundancy" not in asFunction: asFunction["withToleranceInRedundancy"] = 1.e-18
513 if "withLenghtOfRedundancy" not in asFunction: asFunction["withLenghtOfRedundancy"] = -1
514 if "withmpEnabled" not in asFunction: asFunction["withmpEnabled"] = False
515 if "withmpWorkers" not in asFunction: asFunction["withmpWorkers"] = None
516 from daNumerics.ApproximatedDerivatives import FDApproximation
517 FDA = FDApproximation(
518 Function = asFunction["Direct"],
519 centeredDF = asFunction["withCenteredDF"],
520 increment = asFunction["withIncrement"],
521 dX = asFunction["withdX"],
522 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
523 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
524 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
525 mpEnabled = asFunction["withmpEnabled"],
526 mpWorkers = asFunction["withmpWorkers"],
528 self.__CM["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
529 self.__CM["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
530 self.__CM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
531 elif (type(asFunction) is type({})) and \
532 ("Tangent" in asFunction) and ("Adjoint" in asFunction) and \
533 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
534 if ("Direct" not in asFunction) or (asFunction["Direct"] is None):
535 self.__CM["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
537 self.__CM["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
538 self.__CM["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
539 self.__CM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
540 elif asMatrix is not None:
541 matrice = numpy.matrix( asMatrix, numpy.float )
542 self.__CM["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
543 self.__CM["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
544 self.__CM["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
547 raise ValueError("Improperly defined input control model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
550 self.__StoredInputs["ControlModel"] = self.__CM
553 def setControlInput(self,
555 asPersistentVector = None,
560 Permet de définir le controle en entree :
561 - asVector : entrée des données, comme un vecteur compatible avec le
562 constructeur de numpy.matrix
563 - asPersistentVector : entrée des données, comme un vecteur de type
564 persistent contruit avec la classe ad-hoc "Persistence"
565 - Scheduler est le contrôle temporel des données disponibles
566 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
567 être rendue disponible au même titre que les variables de calcul
569 if asVector is not None:
570 self.__U = numpy.matrix( numpy.ravel(numpy.matrix(asVector)), numpy.float ).T
571 elif asPersistentVector is not None:
572 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
573 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.matrix)
574 for member in asPersistentVector:
575 self.__U.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
577 self.__U = asPersistentVector
579 raise ValueError("Error: improperly defined control input, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
581 self.__StoredInputs["ControlInput"] = self.__U
584 # -----------------------------------------------------------
585 def setControls (self,
590 Permet de définir la valeur initiale du vecteur X contenant toutes les
591 variables de contrôle, i.e. les paramètres ou l'état dont on veut
592 estimer la valeur pour obtenir les observations. C'est utile pour un
593 algorithme itératif/incrémental.
594 - asVector : entrée des données, comme un vecteur compatible avec le
595 constructeur de numpy.matrix.
596 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
597 être rendue disponible au même titre que les variables de calcul
599 if asVector is not None:
600 self.__X.store( asVector )
602 self.__StoredInputs["Controls"] = self.__X
605 # -----------------------------------------------------------
606 def setAlgorithm(self, choice = None ):
608 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
609 d'assimilation. L'argument est un champ caractère se rapportant au nom
610 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
611 d'assimilation sur les arguments fixes.
614 raise ValueError("Error: algorithm choice has to be given")
615 if self.__algorithmName is not None:
616 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
617 daDirectory = "daAlgorithms"
619 # Recherche explicitement le fichier complet
620 # ------------------------------------------
622 for directory in sys.path:
623 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
624 module_path = os.path.abspath(os.path.join(directory, daDirectory))
625 if module_path is None:
626 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
628 # Importe le fichier complet comme un module
629 # ------------------------------------------
631 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
632 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
633 self.__algorithmName = str(choice)
634 sys.path = sys_path_tmp ; del sys_path_tmp
635 except ImportError, e:
636 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
638 # Instancie un objet du type élémentaire du fichier
639 # -------------------------------------------------
640 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
641 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
644 def setAlgorithmParameters(self, asDico=None):
646 Permet de définir les paramètres de l'algorithme, sous la forme d'un
649 if asDico is not None:
650 self.__Parameters.update( dict( asDico ) )
652 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
655 def getAlgorithmParameters(self, noDetails=True):
657 Renvoie la liste des paramètres requis selon l'algorithme
659 return self.__algorithm.getRequiredParameters(noDetails)
661 # -----------------------------------------------------------
662 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
664 Permet de sélectionner un diagnostic a effectuer.
667 raise ValueError("Error: diagnostic choice has to be given")
668 daDirectory = "daDiagnostics"
670 # Recherche explicitement le fichier complet
671 # ------------------------------------------
673 for directory in sys.path:
674 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
675 module_path = os.path.abspath(os.path.join(directory, daDirectory))
676 if module_path is None:
677 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
679 # Importe le fichier complet comme un module
680 # ------------------------------------------
682 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
683 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
684 sys.path = sys_path_tmp ; del sys_path_tmp
685 except ImportError, e:
686 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
688 # Instancie un objet du type élémentaire du fichier
689 # -------------------------------------------------
690 if name in self.__StoredInputs:
691 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
692 elif name in self.__StoredDiagnostics:
693 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
695 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
699 parameters = parameters )
702 # -----------------------------------------------------------
703 def shape_validate(self):
705 Validation de la correspondance correcte des tailles des variables et
706 des matrices s'il y en a.
708 if self.__Xb is None: __Xb_shape = (0,)
709 elif hasattr(self.__Xb,"size"): __Xb_shape = (self.__Xb.size,)
710 elif hasattr(self.__Xb,"shape"):
711 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
712 else: __Xb_shape = self.__Xb.shape()
713 else: raise TypeError("The background (Xb) has no attribute of shape: problem !")
715 if self.__Y is None: __Y_shape = (0,)
716 elif hasattr(self.__Y,"size"): __Y_shape = (self.__Y.size,)
717 elif hasattr(self.__Y,"shape"):
718 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
719 else: __Y_shape = self.__Y.shape()
720 else: raise TypeError("The observation (Y) has no attribute of shape: problem !")
722 if self.__U is None: __U_shape = (0,)
723 elif hasattr(self.__U,"size"): __U_shape = (self.__U.size,)
724 elif hasattr(self.__U,"shape"):
725 if type(self.__U.shape) is tuple: __U_shape = self.__U.shape
726 else: __U_shape = self.__U.shape()
727 else: raise TypeError("The control (U) has no attribute of shape: problem !")
729 if self.__B is None: __B_shape = (0,0)
730 elif hasattr(self.__B,"shape"):
731 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
732 else: __B_shape = self.__B.shape()
733 else: raise TypeError("The a priori errors covariance matrix (B) has no attribute of shape: problem !")
735 if self.__R is None: __R_shape = (0,0)
736 elif hasattr(self.__R,"shape"):
737 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
738 else: __R_shape = self.__R.shape()
739 else: raise TypeError("The observation errors covariance matrix (R) has no attribute of shape: problem !")
741 if self.__Q is None: __Q_shape = (0,0)
742 elif hasattr(self.__Q,"shape"):
743 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
744 else: __Q_shape = self.__Q.shape()
745 else: raise TypeError("The evolution errors covariance matrix (Q) has no attribute of shape: problem !")
747 if len(self.__HO) == 0: __HO_shape = (0,0)
748 elif type(self.__HO) is type({}): __HO_shape = (0,0)
749 elif hasattr(self.__HO["Direct"],"shape"):
750 if type(self.__HO["Direct"].shape) is tuple: __HO_shape = self.__HO["Direct"].shape
751 else: __HO_shape = self.__HO["Direct"].shape()
752 else: raise TypeError("The observation operator (H) has no attribute of shape: problem !")
754 if len(self.__EM) == 0: __EM_shape = (0,0)
755 elif type(self.__EM) is type({}): __EM_shape = (0,0)
756 elif hasattr(self.__EM["Direct"],"shape"):
757 if type(self.__EM["Direct"].shape) is tuple: __EM_shape = self.__EM["Direct"].shape
758 else: __EM_shape = self.__EM["Direct"].shape()
759 else: raise TypeError("The evolution model (EM) has no attribute of shape: problem !")
761 if len(self.__CM) == 0: __CM_shape = (0,0)
762 elif type(self.__CM) is type({}): __CM_shape = (0,0)
763 elif hasattr(self.__CM["Direct"],"shape"):
764 if type(self.__CM["Direct"].shape) is tuple: __CM_shape = self.__CM["Direct"].shape
765 else: __CM_shape = self.__CM["Direct"].shape()
766 else: raise TypeError("The control model (CM) has no attribute of shape: problem !")
768 # Vérification des conditions
769 # ---------------------------
770 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
771 raise ValueError("Shape characteristic of background (Xb) is incorrect: \"%s\"."%(__Xb_shape,))
772 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
773 raise ValueError("Shape characteristic of observation (Y) is incorrect: \"%s\"."%(__Y_shape,))
775 if not( min(__B_shape) == max(__B_shape) ):
776 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) is incorrect: \"%s\"."%(__B_shape,))
777 if not( min(__R_shape) == max(__R_shape) ):
778 raise ValueError("Shape characteristic of observation errors covariance matrix (R) is incorrect: \"%s\"."%(__R_shape,))
779 if not( min(__Q_shape) == max(__Q_shape) ):
780 raise ValueError("Shape characteristic of evolution errors covariance matrix (Q) is incorrect: \"%s\"."%(__Q_shape,))
781 if not( min(__EM_shape) == max(__EM_shape) ):
782 raise ValueError("Shape characteristic of evolution operator (EM) is incorrect: \"%s\"."%(__EM_shape,))
784 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
785 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and state (X) \"%s\" are incompatible."%(__HO_shape,__Xb_shape))
786 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
787 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation (Y) \"%s\" are incompatible."%(__HO_shape,__Y_shape))
788 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] ):
789 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and a priori errors covariance matrix (B) \"%s\" are incompatible."%(__HO_shape,__B_shape))
790 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] ):
791 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation errors covariance matrix (R) \"%s\" are incompatible."%(__HO_shape,__R_shape))
793 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
794 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
795 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
796 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
797 for member in asPersistentVector:
798 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
799 __Xb_shape = min(__B_shape)
801 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) \"%s\" and background (Xb) \"%s\" are incompatible."%(__B_shape,__Xb_shape))
803 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
804 raise ValueError("Shape characteristic of observation errors covariance matrix (R) \"%s\" and observation (Y) \"%s\" are incompatible."%(__R_shape,__Y_shape))
806 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) ):
807 raise ValueError("Shape characteristic of evolution model (EM) \"%s\" and state (X) \"%s\" are incompatible."%(__EM_shape,__Xb_shape))
809 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) ):
810 raise ValueError("Shape characteristic of control model (CM) \"%s\" and control (U) \"%s\" are incompatible."%(__CM_shape,__U_shape))
812 if ("AlgorithmParameters" in self.__StoredInputs) \
813 and ("Bounds" in self.__StoredInputs["AlgorithmParameters"]) \
814 and (type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type([]) or type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type(())) \
815 and (len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) != max(__Xb_shape)):
816 raise ValueError("The number \"%s\" of bound pairs for the state (X) components is different of the size \"%s\" of the state itself." \
817 %(len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]),max(__Xb_shape)))
821 # -----------------------------------------------------------
824 Permet de lancer le calcul d'assimilation.
826 Le nom de la méthode à activer est toujours "run". Les paramètres en
827 arguments de la méthode sont fixés. En sortie, on obtient les résultats
828 dans la variable de type dictionnaire "StoredVariables", qui contient en
829 particulier des objets de Persistence pour les analyses, OMA...
831 Operator.CM.clearCache()
832 self.shape_validate()
834 self.__algorithm.run(
844 Parameters = self.__Parameters,
848 # -----------------------------------------------------------
849 def get(self, key=None):
851 Renvoie les résultats disponibles après l'exécution de la méthode
852 d'assimilation, ou les diagnostics disponibles. Attention, quand un
853 diagnostic porte le même nom qu'une variable stockée, c'est la variable
854 stockée qui est renvoyée, et le diagnostic est inatteignable.
857 if key in self.__algorithm:
858 return self.__algorithm.get( key )
859 elif key in self.__StoredInputs:
860 return self.__StoredInputs[key]
861 elif key in self.__StoredDiagnostics:
862 return self.__StoredDiagnostics[key]
864 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
866 allvariables = self.__algorithm.get()
867 allvariables.update( self.__StoredDiagnostics )
868 allvariables.update( self.__StoredInputs )
871 def get_available_variables(self):
873 Renvoie les variables potentiellement utilisables pour l'étude,
874 initialement stockées comme données d'entrées ou dans les algorithmes,
875 identifiés par les chaînes de caractères. L'algorithme doit avoir été
876 préalablement choisi sinon la méthode renvoie "None".
878 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
882 if len( self.__algorithm.keys()) > 0:
883 variables.extend( self.__algorithm.get().keys() )
884 if len( self.__StoredInputs.keys() ) > 0:
885 variables.extend( self.__StoredInputs.keys() )
889 def get_available_algorithms(self):
891 Renvoie la liste des algorithmes 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,"daAlgorithms")):
897 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
898 root, ext = os.path.splitext(fname)
899 if ext == '.py' and root != '__init__':
904 def get_available_diagnostics(self):
906 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
907 par les chaînes de caractères.
910 for directory in sys.path:
911 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
912 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
913 root, ext = os.path.splitext(fname)
914 if ext == '.py' and root != '__init__':
919 # -----------------------------------------------------------
920 def get_algorithms_main_path(self):
922 Renvoie le chemin pour le répertoire principal contenant les algorithmes
923 dans un sous-répertoire "daAlgorithms"
927 def add_algorithms_path(self, asPath=None):
929 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
930 se trouve un sous-répertoire "daAlgorithms"
932 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
933 pas indispensable de le rajouter ici.
935 if not os.path.isdir(asPath):
936 raise ValueError("The given "+asPath+" argument must exist as a directory")
937 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
938 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
939 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
940 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
941 sys.path.insert(0, os.path.abspath(asPath))
942 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
945 def get_diagnostics_main_path(self):
947 Renvoie le chemin pour le répertoire principal contenant les diagnostics
948 dans un sous-répertoire "daDiagnostics"
952 def add_diagnostics_path(self, asPath=None):
954 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
955 se trouve un sous-répertoire "daDiagnostics"
957 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
958 pas indispensable de le rajouter ici.
960 if not os.path.isdir(asPath):
961 raise ValueError("The given "+asPath+" argument must exist as a directory")
962 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
963 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
964 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
965 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
966 sys.path.insert(0, os.path.abspath(asPath))
967 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
970 # -----------------------------------------------------------
971 def setDataObserver(self,
974 HookParameters = None,
978 Permet d'associer un observer à une ou des variables nommées gérées en
979 interne, activable selon des règles définies dans le Scheduler. A chaque
980 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
981 les arguments (variable persistante VariableName, paramètres HookParameters).
984 if type( self.__algorithm ) is dict:
985 raise ValueError("No observer can be build before choosing an algorithm.")
987 # Vérification du nom de variable et typage
988 # -----------------------------------------
989 if type( VariableName ) is str:
990 VariableNames = [VariableName,]
991 elif type( VariableName ) is list:
992 VariableNames = map( str, VariableName )
994 raise ValueError("The observer requires a name or a list of names of variables.")
996 # Association interne de l'observer à la variable
997 # -----------------------------------------------
998 for n in VariableNames:
999 if n not in self.__algorithm:
1000 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
1002 self.__algorithm.StoredVariables[ n ].setDataObserver(
1003 Scheduler = Scheduler,
1004 HookFunction = HookFunction,
1005 HookParameters = HookParameters,
1008 def removeDataObserver(self,
1009 VariableName = None,
1010 HookFunction = None,
1013 Permet de retirer un observer à une ou des variable nommée.
1016 if type( self.__algorithm ) is dict:
1017 raise ValueError("No observer can be removed before choosing an algorithm.")
1019 # Vérification du nom de variable et typage
1020 # -----------------------------------------
1021 if type( VariableName ) is str:
1022 VariableNames = [VariableName,]
1023 elif type( VariableName ) is list:
1024 VariableNames = map( str, VariableName )
1026 raise ValueError("The observer requires a name or a list of names of variables.")
1028 # Association interne de l'observer à la variable
1029 # -----------------------------------------------
1030 for n in VariableNames:
1031 if n not in self.__algorithm:
1032 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
1034 self.__algorithm.StoredVariables[ n ].removeDataObserver(
1035 HookFunction = HookFunction,
1038 # -----------------------------------------------------------
1039 def setDebug(self, level=10):
1041 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
1042 appel pour changer le niveau de verbosité, avec :
1043 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
1045 log = logging.getLogger()
1046 log.setLevel( level )
1048 def unsetDebug(self):
1050 Remet le logger au niveau par défaut
1052 log = logging.getLogger()
1053 log.setLevel( logging.WARNING )
1056 # return set(self.__dict__.keys() + dir(self.__class__))
1057 return ['get', '__doc__', '__init__', '__module__']
1059 def prepare_to_pickle(self):
1061 Retire les variables non pickelisables
1064 del self.__CM # non pickelisable
1065 del self.__EM # non pickelisable
1066 del self.__HO # non pickelisable
1067 del self.__Parameters
1074 del self.__algorithmFile # non pickelisable
1075 del self.__algorithmName
1076 del self.__diagnosticFile # non pickelisable
1077 self.__class__.__doc__ = ""
1079 # ==============================================================================
1080 if __name__ == "__main__":
1081 print '\n AUTODIAGNOSTIC \n'