1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2014 EDF R&D
5 # This library is free software; you can redistribute it and/or
6 # modify it under the terms of the GNU Lesser General Public
7 # License as published by the Free Software Foundation; either
8 # version 2.1 of the License.
10 # This library is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 # Lesser General Public License for more details.
15 # You should have received a copy of the GNU Lesser General Public
16 # License along with this library; if not, write to the Free Software
17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 Classe principale pour la préparation, la réalisation et la restitution de
25 calculs d'assimilation de données.
27 Ce module est destiné à être appelé par AssimilationStudy pour constituer
28 les objets élémentaires de l'étude.
30 __author__ = "Jean-Philippe ARGAUD"
34 import ExtendedLogging ; ExtendedLogging.ExtendedLogging() # A importer en premier
38 logging.debug("Succeed initial import of scipy.optimize with Scipy %s", scipy.version.version)
40 logging.debug("Fail initial import of scipy.optimize")
42 from BasicObjects import Operator, Covariance
43 from PlatformInfo import uniq
45 # ==============================================================================
46 class AssimilationStudy:
48 Cette classe sert d'interface pour l'utilisation de l'assimilation de
49 données. Elle contient les méthodes ou accesseurs nécessaires à la
50 construction d'un calcul d'assimilation.
52 def __init__(self, name=""):
54 Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
55 élémentaire. Ces variables sont ensuite disponibles pour implémenter un
56 algorithme élémentaire particulier.
58 - Background : vecteur Xb
59 - Observation : vecteur Y (potentiellement temporel) d'observations
60 - State : vecteur d'état dont une partie est le vecteur de contrôle.
61 Cette information n'est utile que si l'on veut faire des calculs sur
62 l'état complet, mais elle n'est pas indispensable pour l'assimilation.
63 - Control : vecteur X contenant toutes les variables de contrôle, i.e.
64 les paramètres ou l'état dont on veut estimer la valeur pour obtenir
66 - ObservationOperator...: opérateur d'observation H
68 Les observations présentent une erreur dont la matrice de covariance est
69 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
72 self.__name = str(name)
83 self.__X = Persistence.OneVector()
84 self.__Parameters = {}
85 self.__StoredDiagnostics = {}
86 self.__StoredInputs = {}
88 # Variables temporaires
90 self.__algorithmFile = None
91 self.__algorithmName = None
92 self.__diagnosticFile = None
94 # Récupère le chemin du répertoire parent et l'ajoute au path
95 # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
96 # qui est activée dans Persistence)
97 self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
98 sys.path.insert(0, self.__parent)
99 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
101 # ---------------------------------------------------------
102 def setBackground(self,
104 asPersistentVector = None,
109 Permet de définir l'estimation a priori :
110 - asVector : entrée des données, comme un vecteur compatible avec le
111 constructeur de numpy.matrix
112 - asPersistentVector : entrée des données, comme un vecteur de type
113 persistent contruit avec la classe ad-hoc "Persistence"
114 - Scheduler est le contrôle temporel des données
115 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
116 être rendue disponible au même titre que les variables de calcul
118 if asVector is not None:
119 if type( asVector ) is type( numpy.matrix([]) ):
120 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
122 self.__Xb = numpy.matrix( asVector, numpy.float ).T
123 elif asPersistentVector is not None:
124 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
125 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
126 for member in asPersistentVector:
127 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
129 self.__Xb = asPersistentVector
131 raise ValueError("Error: improperly defined background, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
133 self.__StoredInputs["Background"] = self.__Xb
136 def setBackgroundError(self,
138 asEyeByScalar = None,
139 asEyeByVector = None,
143 Permet de définir la covariance des erreurs d'ébauche :
144 - asCovariance : entrée des données, comme une matrice compatible avec
145 le constructeur de numpy.matrix
146 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
147 multiplicatif d'une matrice de corrélation identité, aucune matrice
148 n'étant donc explicitement à donner
149 - asEyeByVector : entrée des données comme un seul vecteur de variance,
150 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
151 n'étant donc explicitement à donner
152 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
153 être rendue disponible au même titre que les variables de calcul
155 self.__B = Covariance(
156 name = "BackgroundError",
157 asCovariance = asCovariance,
158 asEyeByScalar = asEyeByScalar,
159 asEyeByVector = asEyeByVector,
162 self.__StoredInputs["BackgroundError"] = self.__B
165 # -----------------------------------------------------------
166 def setObservation(self,
168 asPersistentVector = None,
173 Permet de définir les observations :
174 - asVector : entrée des données, comme un vecteur compatible avec le
175 constructeur de numpy.matrix
176 - asPersistentVector : entrée des données, comme un vecteur de type
177 persistent contruit avec la classe ad-hoc "Persistence"
178 - Scheduler est le contrôle temporel des données disponibles
179 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
180 être rendue disponible au même titre que les variables de calcul
182 if asVector is not None:
183 if type( asVector ) is type( numpy.matrix([]) ):
184 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
186 self.__Y = numpy.matrix( asVector, numpy.float ).T
187 elif asPersistentVector is not None:
188 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
189 self.__Y = Persistence.OneVector("Observation", basetype=numpy.matrix)
190 for member in asPersistentVector:
191 self.__Y.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
193 self.__Y = asPersistentVector
195 raise ValueError("Error: improperly defined observations, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
197 self.__StoredInputs["Observation"] = self.__Y
200 def setObservationError(self,
202 asEyeByScalar = None,
203 asEyeByVector = None,
207 Permet de définir la covariance des erreurs d'observations :
208 - asCovariance : entrée des données, comme une matrice compatible avec
209 le constructeur de numpy.matrix
210 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
211 multiplicatif d'une matrice de corrélation identité, aucune matrice
212 n'étant donc explicitement à donner
213 - asEyeByVector : entrée des données comme un seul vecteur de variance,
214 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
215 n'étant donc explicitement à donner
216 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
217 être rendue disponible au même titre que les variables de calcul
219 self.__R = Covariance(
220 name = "ObservationError",
221 asCovariance = asCovariance,
222 asEyeByScalar = asEyeByScalar,
223 asEyeByVector = asEyeByVector,
226 self.__StoredInputs["ObservationError"] = self.__R
229 def setObservationOperator(self,
237 Permet de définir un opérateur d'observation H. L'ordre de priorité des
238 définitions et leur sens sont les suivants :
239 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
240 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
241 "Direct" n'est pas définie, on prend la fonction "Tangent".
242 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
243 des opérateurs tangents et adjoints. On utilise par défaut des
244 différences finies non centrées ou centrées (si "withCenteredDF" est
245 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
246 du point courant ou sur le point fixe "withdX".
247 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
248 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
249 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
250 matrice fournie doit être sous une forme compatible avec le
251 constructeur de numpy.matrix.
252 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
253 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
254 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
255 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
256 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
257 être rendue disponible au même titre que les variables de calcul
258 L'argument "asFunction" peut prendre la forme complète suivante, avec
259 les valeurs par défaut standards :
260 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
261 "useApproximatedDerivatives":False,
262 "withCenteredDF" :False,
263 "withIncrement" :0.01,
265 "withAvoidingRedundancy" :True,
266 "withToleranceInRedundancy" :1.e-18,
267 "withLenghtOfRedundancy" :-1,
268 "withmpEnabled" :False,
269 "withmpWorkers" :None,
272 if (type(asFunction) is type({})) and \
273 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
274 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
275 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
276 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
277 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
278 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
279 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
280 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
281 if not asFunction.has_key("withmpEnabled"): asFunction["withmpEnabled"] = False
282 if not asFunction.has_key("withmpWorkers"): asFunction["withmpWorkers"] = None
283 from daNumerics.ApproximatedDerivatives import FDApproximation
284 FDA = FDApproximation(
285 Function = asFunction["Direct"],
286 centeredDF = asFunction["withCenteredDF"],
287 increment = asFunction["withIncrement"],
288 dX = asFunction["withdX"],
289 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
290 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
291 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
292 mpEnabled = asFunction["withmpEnabled"],
293 mpWorkers = asFunction["withmpWorkers"],
295 self.__HO["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
296 self.__HO["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
297 self.__HO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
298 elif (type(asFunction) is type({})) and \
299 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
300 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
301 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
302 self.__HO["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
304 self.__HO["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
305 self.__HO["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
306 self.__HO["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
307 elif asMatrix is not None:
308 matrice = numpy.matrix( asMatrix, numpy.float )
309 self.__HO["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
310 self.__HO["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
311 self.__HO["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
314 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
316 if appliedToX is not None:
317 self.__HO["AppliedToX"] = {}
318 if type(appliedToX) is not dict:
319 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
320 for key in appliedToX.keys():
321 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
322 # Pour le cas où l'on a une vraie matrice
323 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
324 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
325 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
326 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
328 self.__HO["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
330 self.__HO["AppliedToX"] = None
333 self.__StoredInputs["ObservationOperator"] = self.__HO
336 # -----------------------------------------------------------
337 def setEvolutionModel(self,
345 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
346 définitions et leur sens sont les suivants :
347 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
348 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
349 "Direct" n'est pas définie, on prend la fonction "Tangent".
350 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
351 des opérateurs tangents et adjoints. On utilise par défaut des
352 différences finies non centrées ou centrées (si "withCenteredDF" est
353 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
354 du point courant ou sur le point fixe "withdX".
355 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
356 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
357 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
358 matrice fournie doit être sous une forme compatible avec le
359 constructeur de numpy.matrix.
360 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
361 être rendue disponible au même titre que les variables de calcul
362 L'argument "asFunction" peut prendre la forme complète suivante, avec
363 les valeurs par défaut standards :
364 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
365 "useApproximatedDerivatives":False,
366 "withCenteredDF" :False,
367 "withIncrement" :0.01,
369 "withAvoidingRedundancy" :True,
370 "withToleranceInRedundancy" :1.e-18,
371 "withLenghtOfRedundancy" :-1,
372 "withmpEnabled" :False,
373 "withmpWorkers" :None,
376 if (type(asFunction) is type({})) and \
377 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
378 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
379 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
380 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
381 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
382 if not asFunction.has_key("withAvoidingRedundancy"): asFunction["withAvoidingRedundancy"] = True
383 if not asFunction.has_key("withToleranceInRedundancy"): asFunction["withToleranceInRedundancy"] = 1.e-18
384 if not asFunction.has_key("withLenghtOfRedundancy"): asFunction["withLenghtOfRedundancy"] = -1
385 if not asFunction.has_key("withmpEnabled"): asFunction["withmpEnabled"] = False
386 if not asFunction.has_key("withmpWorkers"): asFunction["withmpWorkers"] = None
387 from daNumerics.ApproximatedDerivatives import FDApproximation
388 FDA = FDApproximation(
389 Function = asFunction["Direct"],
390 centeredDF = asFunction["withCenteredDF"],
391 increment = asFunction["withIncrement"],
392 dX = asFunction["withdX"],
393 avoidingRedundancy = asFunction["withAvoidingRedundancy"],
394 toleranceInRedundancy = asFunction["withToleranceInRedundancy"],
395 lenghtOfRedundancy = asFunction["withLenghtOfRedundancy"],
396 mpEnabled = asFunction["withmpEnabled"],
397 mpWorkers = asFunction["withmpWorkers"],
399 self.__EM["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
400 self.__EM["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
401 self.__EM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
402 elif (type(asFunction) is type({})) and \
403 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
404 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
405 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
406 self.__EM["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
408 self.__EM["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
409 self.__EM["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
410 self.__EM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
411 elif asMatrix is not None:
412 matrice = numpy.matrix( asMatrix, numpy.float )
413 self.__EM["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
414 self.__EM["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
415 self.__EM["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
418 raise ValueError("Improperly defined evolution model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
421 self.__StoredInputs["EvolutionModel"] = self.__EM
424 def setEvolutionError(self,
426 asEyeByScalar = None,
427 asEyeByVector = None,
431 Permet de définir la covariance des erreurs de modèle :
432 - asCovariance : entrée des données, comme une matrice compatible avec
433 le constructeur de numpy.matrix
434 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
435 multiplicatif d'une matrice de corrélation identité, aucune matrice
436 n'étant donc explicitement à donner
437 - asEyeByVector : entrée des données comme un seul vecteur de variance,
438 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
439 n'étant donc explicitement à donner
440 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
441 être rendue disponible au même titre que les variables de calcul
443 self.__Q = Covariance(
444 name = "EvolutionError",
445 asCovariance = asCovariance,
446 asEyeByScalar = asEyeByScalar,
447 asEyeByVector = asEyeByVector,
450 self.__StoredInputs["EvolutionError"] = self.__Q
453 # -----------------------------------------------------------
454 def setControlModel(self,
455 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
456 "useApproximatedDerivatives":False,
457 "withCenteredDF" :False,
458 "withIncrement" :0.01,
467 Permet de définir un opérateur de controle C. L'ordre de priorité des
468 définitions et leur sens sont les suivants :
469 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
470 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
471 "Direct" n'est pas définie, on prend la fonction "Tangent".
472 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
473 des opérateurs tangents et adjoints. On utilise par défaut des
474 différences finies non centrées ou centrées (si "withCenteredDF" est
475 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
476 du point courant ou sur le point fixe "withdX".
477 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
478 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
479 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
480 matrice fournie doit être sous une forme compatible avec le
481 constructeur de numpy.matrix.
482 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
483 être rendue disponible au même titre que les variables de calcul
485 if (type(asFunction) is type({})) and \
486 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
487 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
488 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
489 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
490 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
491 from daNumerics.ApproximatedDerivatives import FDApproximation
492 FDA = FDApproximation(
493 Function = asFunction["Direct"],
494 centeredDF = asFunction["withCenteredDF"],
495 increment = asFunction["withIncrement"],
496 dX = asFunction["withdX"] )
497 self.__CM["Direct"] = Operator( fromMethod = FDA.DirectOperator, avoidingRedundancy = avoidRC )
498 self.__CM["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC )
499 self.__CM["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC )
500 elif (type(asFunction) is type({})) and \
501 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
502 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
503 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
504 self.__CM["Direct"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
506 self.__CM["Direct"] = Operator( fromMethod = asFunction["Direct"], avoidingRedundancy = avoidRC )
507 self.__CM["Tangent"] = Operator( fromMethod = asFunction["Tangent"], avoidingRedundancy = avoidRC )
508 self.__CM["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"], avoidingRedundancy = avoidRC )
509 elif asMatrix is not None:
510 matrice = numpy.matrix( asMatrix, numpy.float )
511 self.__CM["Direct"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
512 self.__CM["Tangent"] = Operator( fromMatrix = matrice, avoidingRedundancy = avoidRC )
513 self.__CM["Adjoint"] = Operator( fromMatrix = matrice.T, avoidingRedundancy = avoidRC )
516 raise ValueError("Improperly defined input control model, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
519 self.__StoredInputs["ControlModel"] = self.__CM
522 def setControlInput(self,
524 asPersistentVector = None,
529 Permet de définir le controle en entree :
530 - asVector : entrée des données, comme un vecteur compatible avec le
531 constructeur de numpy.matrix
532 - asPersistentVector : entrée des données, comme un vecteur de type
533 persistent contruit avec la classe ad-hoc "Persistence"
534 - Scheduler est le contrôle temporel des données disponibles
535 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
536 être rendue disponible au même titre que les variables de calcul
538 if asVector is not None:
539 if isinstance(asVector,numpy.matrix):
540 self.__U = numpy.matrix( asVector.A1, numpy.float ).T
542 self.__U = numpy.matrix( asVector, numpy.float ).T
543 elif asPersistentVector is not None:
544 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
545 self.__U = Persistence.OneVector("ControlInput", basetype=numpy.matrix)
546 for member in asPersistentVector:
547 self.__U.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
549 self.__U = asPersistentVector
551 raise ValueError("Error: improperly defined control input, it requires at minima either a vector, a list/tuple of vectors or a persistent object")
553 self.__StoredInputs["ControlInput"] = self.__U
556 # -----------------------------------------------------------
557 def setControls (self,
562 Permet de définir la valeur initiale du vecteur X contenant toutes les
563 variables de contrôle, i.e. les paramètres ou l'état dont on veut
564 estimer la valeur pour obtenir les observations. C'est utile pour un
565 algorithme itératif/incrémental.
566 - asVector : entrée des données, comme un vecteur compatible avec le
567 constructeur de numpy.matrix.
568 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
569 être rendue disponible au même titre que les variables de calcul
571 if asVector is not None:
572 self.__X.store( asVector )
574 self.__StoredInputs["Controls"] = self.__X
577 # -----------------------------------------------------------
578 def setAlgorithm(self, choice = None ):
580 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
581 d'assimilation. L'argument est un champ caractère se rapportant au nom
582 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
583 d'assimilation sur les arguments fixes.
586 raise ValueError("Error: algorithm choice has to be given")
587 if self.__algorithmName is not None:
588 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
589 daDirectory = "daAlgorithms"
591 # Recherche explicitement le fichier complet
592 # ------------------------------------------
594 for directory in sys.path:
595 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
596 module_path = os.path.abspath(os.path.join(directory, daDirectory))
597 if module_path is None:
598 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
600 # Importe le fichier complet comme un module
601 # ------------------------------------------
603 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
604 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
605 self.__algorithmName = str(choice)
606 sys.path = sys_path_tmp ; del sys_path_tmp
607 except ImportError, e:
608 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
610 # Instancie un objet du type élémentaire du fichier
611 # -------------------------------------------------
612 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
613 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
616 def setAlgorithmParameters(self, asDico=None):
618 Permet de définir les paramètres de l'algorithme, sous la forme d'un
621 if asDico is not None:
622 self.__Parameters.update( dict( asDico ) )
624 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
627 def getAlgorithmParameters(self, noDetails=True):
629 Renvoie la liste des paramètres requis selon l'algorithme
631 return self.__algorithm.getRequiredParameters(noDetails)
633 # -----------------------------------------------------------
634 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
636 Permet de sélectionner un diagnostic a effectuer.
639 raise ValueError("Error: diagnostic choice has to be given")
640 daDirectory = "daDiagnostics"
642 # Recherche explicitement le fichier complet
643 # ------------------------------------------
645 for directory in sys.path:
646 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
647 module_path = os.path.abspath(os.path.join(directory, daDirectory))
648 if module_path is None:
649 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
651 # Importe le fichier complet comme un module
652 # ------------------------------------------
654 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
655 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
656 sys.path = sys_path_tmp ; del sys_path_tmp
657 except ImportError, e:
658 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
660 # Instancie un objet du type élémentaire du fichier
661 # -------------------------------------------------
662 if self.__StoredInputs.has_key(name):
663 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
664 elif self.__StoredDiagnostics.has_key(name):
665 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
667 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
671 parameters = parameters )
674 # -----------------------------------------------------------
675 def shape_validate(self):
677 Validation de la correspondance correcte des tailles des variables et
678 des matrices s'il y en a.
680 if self.__Xb is None: __Xb_shape = (0,)
681 elif hasattr(self.__Xb,"size"): __Xb_shape = (self.__Xb.size,)
682 elif hasattr(self.__Xb,"shape"):
683 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
684 else: __Xb_shape = self.__Xb.shape()
685 else: raise TypeError("The background (Xb) has no attribute of shape: problem !")
687 if self.__Y is None: __Y_shape = (0,)
688 elif hasattr(self.__Y,"size"): __Y_shape = (self.__Y.size,)
689 elif hasattr(self.__Y,"shape"):
690 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
691 else: __Y_shape = self.__Y.shape()
692 else: raise TypeError("The observation (Y) has no attribute of shape: problem !")
694 if self.__U is None: __U_shape = (0,)
695 elif hasattr(self.__U,"size"): __U_shape = (self.__U.size,)
696 elif hasattr(self.__U,"shape"):
697 if type(self.__U.shape) is tuple: __U_shape = self.__U.shape
698 else: __U_shape = self.__U.shape()
699 else: raise TypeError("The control (U) has no attribute of shape: problem !")
701 if self.__B is None: __B_shape = (0,0)
702 elif hasattr(self.__B,"shape"):
703 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
704 else: __B_shape = self.__B.shape()
705 else: raise TypeError("The a priori errors covariance matrix (B) has no attribute of shape: problem !")
707 if self.__R is None: __R_shape = (0,0)
708 elif hasattr(self.__R,"shape"):
709 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
710 else: __R_shape = self.__R.shape()
711 else: raise TypeError("The observation errors covariance matrix (R) has no attribute of shape: problem !")
713 if self.__Q is None: __Q_shape = (0,0)
714 elif hasattr(self.__Q,"shape"):
715 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
716 else: __Q_shape = self.__Q.shape()
717 else: raise TypeError("The evolution errors covariance matrix (Q) has no attribute of shape: problem !")
719 if len(self.__HO) == 0: __HO_shape = (0,0)
720 elif type(self.__HO) is type({}): __HO_shape = (0,0)
721 elif hasattr(self.__HO["Direct"],"shape"):
722 if type(self.__HO["Direct"].shape) is tuple: __HO_shape = self.__HO["Direct"].shape
723 else: __HO_shape = self.__HO["Direct"].shape()
724 else: raise TypeError("The observation operator (H) has no attribute of shape: problem !")
726 if len(self.__EM) == 0: __EM_shape = (0,0)
727 elif type(self.__EM) is type({}): __EM_shape = (0,0)
728 elif hasattr(self.__EM["Direct"],"shape"):
729 if type(self.__EM["Direct"].shape) is tuple: __EM_shape = self.__EM["Direct"].shape
730 else: __EM_shape = self.__EM["Direct"].shape()
731 else: raise TypeError("The evolution model (EM) has no attribute of shape: problem !")
733 if len(self.__CM) == 0: __CM_shape = (0,0)
734 elif type(self.__CM) is type({}): __CM_shape = (0,0)
735 elif hasattr(self.__CM["Direct"],"shape"):
736 if type(self.__CM["Direct"].shape) is tuple: __CM_shape = self.__CM["Direct"].shape
737 else: __CM_shape = self.__CM["Direct"].shape()
738 else: raise TypeError("The control model (CM) has no attribute of shape: problem !")
740 # Vérification des conditions
741 # ---------------------------
742 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
743 raise ValueError("Shape characteristic of background (Xb) is incorrect: \"%s\"."%(__Xb_shape,))
744 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
745 raise ValueError("Shape characteristic of observation (Y) is incorrect: \"%s\"."%(__Y_shape,))
747 if not( min(__B_shape) == max(__B_shape) ):
748 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) is incorrect: \"%s\"."%(__B_shape,))
749 if not( min(__R_shape) == max(__R_shape) ):
750 raise ValueError("Shape characteristic of observation errors covariance matrix (R) is incorrect: \"%s\"."%(__R_shape,))
751 if not( min(__Q_shape) == max(__Q_shape) ):
752 raise ValueError("Shape characteristic of evolution errors covariance matrix (Q) is incorrect: \"%s\"."%(__Q_shape,))
753 if not( min(__EM_shape) == max(__EM_shape) ):
754 raise ValueError("Shape characteristic of evolution operator (EM) is incorrect: \"%s\"."%(__EM_shape,))
756 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[1] == max(__Xb_shape) ):
757 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and state (X) \"%s\" are incompatible."%(__HO_shape,__Xb_shape))
758 if len(self.__HO) > 0 and not(type(self.__HO) is type({})) and not( __HO_shape[0] == max(__Y_shape) ):
759 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation (Y) \"%s\" are incompatible."%(__HO_shape,__Y_shape))
760 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] ):
761 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and a priori errors covariance matrix (B) \"%s\" are incompatible."%(__HO_shape,__B_shape))
762 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] ):
763 raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation errors covariance matrix (R) \"%s\" are incompatible."%(__HO_shape,__R_shape))
765 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
766 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
767 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
768 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
769 for member in asPersistentVector:
770 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
771 __Xb_shape = min(__B_shape)
773 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) \"%s\" and background (Xb) \"%s\" are incompatible."%(__B_shape,__Xb_shape))
775 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
776 raise ValueError("Shape characteristic of observation errors covariance matrix (R) \"%s\" and observation (Y) \"%s\" are incompatible."%(__R_shape,__Y_shape))
778 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) ):
779 raise ValueError("Shape characteristic of evolution model (EM) \"%s\" and state (X) \"%s\" are incompatible."%(__EM_shape,__Xb_shape))
781 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) ):
782 raise ValueError("Shape characteristic of control model (CM) \"%s\" and control (U) \"%s\" are incompatible."%(__CM_shape,__U_shape))
784 if self.__StoredInputs.has_key("AlgorithmParameters") \
785 and self.__StoredInputs["AlgorithmParameters"].has_key("Bounds") \
786 and (type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type([]) or type(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) is type(())) \
787 and (len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]) != max(__Xb_shape)):
788 raise ValueError("The number \"%s\" of bound pairs for the state (X) components is different of the size \"%s\" of the state itself." \
789 %(len(self.__StoredInputs["AlgorithmParameters"]["Bounds"]),max(__Xb_shape)))
793 # -----------------------------------------------------------
796 Permet de lancer le calcul d'assimilation.
798 Le nom de la méthode à activer est toujours "run". Les paramètres en
799 arguments de la méthode sont fixés. En sortie, on obtient les résultats
800 dans la variable de type dictionnaire "StoredVariables", qui contient en
801 particulier des objets de Persistence pour les analyses, OMA...
803 Operator.CM.clearCache()
804 self.shape_validate()
806 self.__algorithm.run(
816 Parameters = self.__Parameters,
820 # -----------------------------------------------------------
821 def get(self, key=None):
823 Renvoie les résultats disponibles après l'exécution de la méthode
824 d'assimilation, ou les diagnostics disponibles. Attention, quand un
825 diagnostic porte le même nom qu'une variable stockée, c'est la variable
826 stockée qui est renvoyée, et le diagnostic est inatteignable.
829 if self.__algorithm.has_key(key):
830 return self.__algorithm.get( key )
831 elif self.__StoredInputs.has_key(key):
832 return self.__StoredInputs[key]
833 elif self.__StoredDiagnostics.has_key(key):
834 return self.__StoredDiagnostics[key]
836 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
838 allvariables = self.__algorithm.get()
839 allvariables.update( self.__StoredDiagnostics )
840 allvariables.update( self.__StoredInputs )
843 def get_available_variables(self):
845 Renvoie les variables potentiellement utilisables pour l'étude,
846 initialement stockées comme données d'entrées ou dans les algorithmes,
847 identifiés par les chaînes de caractères. L'algorithme doit avoir été
848 préalablement choisi sinon la méthode renvoie "None".
850 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
854 if len( self.__algorithm.keys()) > 0:
855 variables.extend( self.__algorithm.get().keys() )
856 if len( self.__StoredInputs.keys() ) > 0:
857 variables.extend( self.__StoredInputs.keys() )
861 def get_available_algorithms(self):
863 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
864 par les chaînes de caractères.
867 for directory in sys.path:
868 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
869 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
870 root, ext = os.path.splitext(fname)
871 if ext == '.py' and root != '__init__':
876 def get_available_diagnostics(self):
878 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
879 par les chaînes de caractères.
882 for directory in sys.path:
883 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
884 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
885 root, ext = os.path.splitext(fname)
886 if ext == '.py' and root != '__init__':
891 # -----------------------------------------------------------
892 def get_algorithms_main_path(self):
894 Renvoie le chemin pour le répertoire principal contenant les algorithmes
895 dans un sous-répertoire "daAlgorithms"
899 def add_algorithms_path(self, asPath=None):
901 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
902 se trouve un sous-répertoire "daAlgorithms"
904 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
905 pas indispensable de le rajouter ici.
907 if not os.path.isdir(asPath):
908 raise ValueError("The given "+asPath+" argument must exist as a directory")
909 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
910 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
911 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
912 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
913 sys.path.insert(0, os.path.abspath(asPath))
914 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
917 def get_diagnostics_main_path(self):
919 Renvoie le chemin pour le répertoire principal contenant les diagnostics
920 dans un sous-répertoire "daDiagnostics"
924 def add_diagnostics_path(self, asPath=None):
926 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
927 se trouve un sous-répertoire "daDiagnostics"
929 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
930 pas indispensable de le rajouter ici.
932 if not os.path.isdir(asPath):
933 raise ValueError("The given "+asPath+" argument must exist as a directory")
934 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
935 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
936 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
937 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
938 sys.path.insert(0, os.path.abspath(asPath))
939 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
942 # -----------------------------------------------------------
943 def setDataObserver(self,
946 HookParameters = None,
950 Permet d'associer un observer à une ou des variables nommées gérées en
951 interne, activable selon des règles définies dans le Scheduler. A chaque
952 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
953 les arguments (variable persistante VariableName, paramètres HookParameters).
956 if type( self.__algorithm ) is dict:
957 raise ValueError("No observer can be build before choosing an algorithm.")
959 # Vérification du nom de variable et typage
960 # -----------------------------------------
961 if type( VariableName ) is str:
962 VariableNames = [VariableName,]
963 elif type( VariableName ) is list:
964 VariableNames = map( str, VariableName )
966 raise ValueError("The observer requires a name or a list of names of variables.")
968 # Association interne de l'observer à la variable
969 # -----------------------------------------------
970 for n in VariableNames:
971 if not self.__algorithm.has_key( n ):
972 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
974 self.__algorithm.StoredVariables[ n ].setDataObserver(
975 Scheduler = Scheduler,
976 HookFunction = HookFunction,
977 HookParameters = HookParameters,
980 def removeDataObserver(self,
985 Permet de retirer un observer à une ou des variable nommée.
988 if type( self.__algorithm ) is dict:
989 raise ValueError("No observer can be removed before choosing an algorithm.")
991 # Vérification du nom de variable et typage
992 # -----------------------------------------
993 if type( VariableName ) is str:
994 VariableNames = [VariableName,]
995 elif type( VariableName ) is list:
996 VariableNames = map( str, VariableName )
998 raise ValueError("The observer requires a name or a list of names of variables.")
1000 # Association interne de l'observer à la variable
1001 # -----------------------------------------------
1002 for n in VariableNames:
1003 if not self.__algorithm.has_key( n ):
1004 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
1006 self.__algorithm.StoredVariables[ n ].removeDataObserver(
1007 HookFunction = HookFunction,
1010 # -----------------------------------------------------------
1011 def setDebug(self, level=10):
1013 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
1014 appel pour changer le niveau de verbosité, avec :
1015 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
1017 log = logging.getLogger()
1018 log.setLevel( level )
1020 def unsetDebug(self):
1022 Remet le logger au niveau par défaut
1024 log = logging.getLogger()
1025 log.setLevel( logging.WARNING )
1027 def prepare_to_pickle(self):
1029 Retire les variables non pickelisables
1031 self.__algorithmFile = None
1032 self.__diagnosticFile = None
1037 # ==============================================================================
1038 if __name__ == "__main__":
1039 print '\n AUTODIAGNOSTIC \n'