1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2012 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 Définit les outils généraux élémentaires.
26 Ce module est destiné à etre appelée par AssimilationStudy pour constituer
27 les objets élémentaires de l'étude.
29 __author__ = "Jean-Philippe ARGAUD"
33 import Logging ; Logging.Logging() # A importer en premier
35 from BasicObjects import Operator
37 # ==============================================================================
38 class AssimilationStudy:
40 Cette classe sert d'interface pour l'utilisation de l'assimilation de
41 données. Elle contient les méthodes ou accesseurs nécessaires à la
42 construction d'un calcul d'assimilation.
44 def __init__(self, name=""):
46 Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
47 élémentaire. Ces variables sont ensuite disponibles pour implémenter un
48 algorithme élémentaire particulier.
50 Background............: vecteur Xb
51 Observation...........: vecteur Y (potentiellement temporel)
53 State.................: vecteur d'état dont une partie est le vecteur de
54 contrôle. Cette information n'est utile que si l'on veut faire des
55 calculs sur l'état complet, mais elle n'est pas indispensable pour
57 Control...............: vecteur X contenant toutes les variables de
58 contrôle, i.e. les paramètres ou l'état dont on veut estimer la
59 valeur pour obtenir les observations
60 ObservationOperator...: opérateur d'observation H
62 Les observations présentent une erreur dont la matrice de covariance est
63 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
66 self.__name = str(name)
75 self.__X = Persistence.OneVector()
76 self.__Parameters = {}
77 self.__StoredDiagnostics = {}
78 self.__StoredInputs = {}
80 self.__B_scalar = None
81 self.__R_scalar = None
82 self.__Q_scalar = None
83 self.__Parameters["B_scalar"] = None
84 self.__Parameters["R_scalar"] = None
85 self.__Parameters["Q_scalar"] = None
87 # Variables temporaires
89 self.__algorithmFile = None
90 self.__algorithmName = None
91 self.__diagnosticFile = None
93 # Récupère le chemin du répertoire parent et l'ajoute au path
94 # (Cela complète l'action de la classe PathManagement dans PlatformInfo,
95 # qui est activée dans Persistence)
96 self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
97 sys.path.insert(0, self.__parent)
98 sys.path = list(set(sys.path)) # Conserve en unique exemplaire chaque chemin
100 # ---------------------------------------------------------
101 def setBackground(self,
103 asPersistentVector = None,
108 Permet de définir l'estimation a priori :
109 - asVector : entrée des données, comme un vecteur compatible avec le
110 constructeur de numpy.matrix
111 - asPersistentVector : entrée des données, comme un vecteur de type
112 persistent contruit avec la classe ad-hoc "Persistence"
113 - Scheduler est le contrôle temporel des données
114 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
115 être rendue disponible au même titre que les variables de calcul
117 if asVector is not None:
118 if type( asVector ) is type( numpy.matrix([]) ):
119 self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
121 self.__Xb = numpy.matrix( asVector, numpy.float ).T
122 elif asPersistentVector is not None:
123 if type(asPersistentVector) in [type([]),type(()),type(numpy.array([])),type(numpy.matrix([]))]:
124 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
125 for member in asPersistentVector:
126 self.__Xb.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
128 self.__Xb = asPersistentVector
130 raise ValueError("Error: improperly defined background")
132 self.__StoredInputs["Background"] = self.__Xb
135 def setBackgroundError(self,
137 asEyeByScalar = None,
138 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 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
152 être rendue disponible au même titre que les variables de calcul
154 if asEyeByScalar is not None:
155 self.__B_scalar = float(asEyeByScalar)
157 elif asEyeByVector is not None:
158 self.__B_scalar = None
159 self.__B = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
160 elif asCovariance is not None:
161 self.__B_scalar = None
162 self.__B = numpy.matrix( asCovariance, numpy.float )
164 raise ValueError("Background error covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix")
166 self.__Parameters["B_scalar"] = self.__B_scalar
168 self.__StoredInputs["BackgroundError"] = self.__B
171 # -----------------------------------------------------------
172 def setObservation(self,
174 asPersistentVector = None,
179 Permet de définir les observations :
180 - asVector : entrée des données, comme un vecteur compatible avec le
181 constructeur de numpy.matrix
182 - asPersistentVector : entrée des données, comme un vecteur de type
183 persistent contruit avec la classe ad-hoc "Persistence"
184 - Scheduler est le contrôle temporel des données disponibles
185 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
186 être rendue disponible au même titre que les variables de calcul
188 if asVector is not None:
189 if type( asVector ) is type( numpy.matrix([]) ):
190 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
192 self.__Y = numpy.matrix( asVector, numpy.float ).T
193 elif asPersistentVector is not None:
194 self.__Y = asPersistentVector
196 raise ValueError("Error: improperly defined observations")
198 self.__StoredInputs["Observation"] = self.__Y
201 def setObservationError(self,
203 asEyeByScalar = None,
204 asEyeByVector = None,
208 Permet de définir la covariance des erreurs d'observations :
209 - asCovariance : entrée des données, comme une matrice compatible avec
210 le constructeur de numpy.matrix
211 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
212 multiplicatif d'une matrice de corrélation identité, aucune matrice
213 n'étant donc explicitement à donner
214 - asEyeByVector : entrée des données comme un seul vecteur de variance,
215 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
216 n'étant donc explicitement à donner
217 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
218 être rendue disponible au même titre que les variables de calcul
220 if asEyeByScalar is not None:
221 self.__R_scalar = float(asEyeByScalar)
223 elif asEyeByVector is not None:
224 self.__R_scalar = None
225 self.__R = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
226 elif asCovariance is not None:
227 self.__R_scalar = None
228 self.__R = numpy.matrix( asCovariance, numpy.float )
230 raise ValueError("Observation error covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix")
232 self.__Parameters["R_scalar"] = self.__R_scalar
234 self.__StoredInputs["ObservationError"] = self.__R
237 def setObservationOperator(self,
238 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
239 "useApproximatedDerivatives":False,
240 "withCenteredDF" :False,
241 "withIncrement" :0.01,
249 Permet de définir un opérateur d'observation H. L'ordre de priorité des
250 définitions et leur sens sont les suivants :
251 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
252 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
253 "Direct" n'est pas définie, on prend la fonction "Tangent".
254 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
255 des opérateurs tangents et adjoints. On utilise par défaut des
256 différences finies non centrées ou centrées (si "withCenteredDF" est
257 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
258 du point courant ou sur le point fixe "withdX".
259 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
260 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
261 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
262 matrice fournie doit être sous une forme compatible avec le
263 constructeur de numpy.matrix.
264 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
265 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
266 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
267 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
268 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
269 être rendue disponible au même titre que les variables de calcul
271 if (type(asFunction) is type({})) and \
272 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
273 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
274 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
275 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
276 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
277 from daNumerics.ApproximatedDerivatives import FDApproximation
278 FDA = FDApproximation(
279 FunctionH = asFunction["Direct"],
280 centeredDF = asFunction["withCenteredDF"],
281 increment = asFunction["withIncrement"],
282 dX = asFunction["withdX"] )
283 self.__H["Direct"] = Operator( fromMethod = FDA.FunctionH )
284 self.__H["Tangent"] = Operator( fromMethod = FDA.TangentH )
285 self.__H["Adjoint"] = Operator( fromMethod = FDA.AdjointH )
286 elif (type(asFunction) is type({})) and \
287 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
288 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
289 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
290 self.__H["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
292 self.__H["Direct"] = Operator( fromMethod = asFunction["Direct"] )
293 self.__H["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
294 self.__H["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
295 elif asMatrix is not None:
296 matrice = numpy.matrix( asMatrix, numpy.float )
297 self.__H["Direct"] = Operator( fromMatrix = matrice )
298 self.__H["Tangent"] = Operator( fromMatrix = matrice )
299 self.__H["Adjoint"] = Operator( fromMatrix = matrice.T )
301 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
303 if appliedToX is not None:
304 self.__H["AppliedToX"] = {}
305 if type(appliedToX) is not dict:
306 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
307 for key in appliedToX.keys():
308 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
309 # Pour le cas où l'on a une vraie matrice
310 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
311 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
312 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
313 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
315 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
317 self.__H["AppliedToX"] = None
320 self.__StoredInputs["ObservationOperator"] = self.__H
323 # -----------------------------------------------------------
324 def setEvolutionModel(self,
325 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
326 "useApproximatedDerivatives":False,
327 "withCenteredDF" :False,
328 "withIncrement" :0.01,
336 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
337 définitions et leur sens sont les suivants :
338 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
339 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
340 "Direct" n'est pas définie, on prend la fonction "Tangent".
341 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
342 des opérateurs tangents et adjoints. On utilise par défaut des
343 différences finies non centrées ou centrées (si "withCenteredDF" est
344 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
345 du point courant ou sur le point fixe "withdX".
346 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
347 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
348 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
349 matrice fournie doit être sous une forme compatible avec le
350 constructeur de numpy.matrix.
351 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
352 être rendue disponible au même titre que les variables de calcul
354 if (type(asFunction) is type({})) and \
355 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
356 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
357 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
358 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
359 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
360 from daNumerics.ApproximatedDerivatives import FDApproximation
361 FDA = FDApproximation(
362 FunctionH = asFunction["Direct"],
363 centeredDF = asFunction["withCenteredDF"],
364 increment = asFunction["withIncrement"],
365 dX = asFunction["withdX"] )
366 self.__H["Direct"] = Operator( fromMethod = FDA.FunctionH )
367 self.__H["Tangent"] = Operator( fromMethod = FDA.TangentH )
368 self.__H["Adjoint"] = Operator( fromMethod = FDA.AdjointH )
369 elif (type(asFunction) is type({})) and \
370 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
371 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
372 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
373 self.__M["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
375 self.__M["Direct"] = Operator( fromMethod = asFunction["Direct"] )
376 self.__M["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
377 self.__M["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
378 elif asMatrix is not None:
379 matrice = numpy.matrix( asMatrix, numpy.float )
380 self.__M["Direct"] = Operator( fromMatrix = matrice )
381 self.__M["Tangent"] = Operator( fromMatrix = matrice )
382 self.__M["Adjoint"] = Operator( fromMatrix = matrice.T )
384 raise ValueError("Improperly defined evolution operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
387 self.__StoredInputs["EvolutionModel"] = self.__M
390 def setEvolutionError(self,
392 asEyeByScalar = None,
393 asEyeByVector = None,
397 Permet de définir la covariance des erreurs de modèle :
398 - asCovariance : entrée des données, comme une matrice compatible avec
399 le constructeur de numpy.matrix
400 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
401 multiplicatif d'une matrice de corrélation identité, aucune matrice
402 n'étant donc explicitement à donner
403 - asEyeByVector : entrée des données comme un seul vecteur de variance,
404 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
405 n'étant donc explicitement à donner
406 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
407 être rendue disponible au même titre que les variables de calcul
409 if asEyeByScalar is not None:
410 self.__Q_scalar = float(asEyeByScalar)
412 elif asEyeByVector is not None:
413 self.__Q_scalar = None
414 self.__Q = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
415 elif asCovariance is not None:
416 self.__Q_scalar = None
417 self.__Q = numpy.matrix( asCovariance, numpy.float )
419 raise ValueError("Evolution error covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix")
421 self.__Parameters["Q_scalar"] = self.__Q_scalar
423 self.__StoredInputs["EvolutionError"] = self.__Q
426 # -----------------------------------------------------------
427 def setControls (self,
432 Permet de définir la valeur initiale du vecteur X contenant toutes les
433 variables de contrôle, i.e. les paramètres ou l'état dont on veut
434 estimer la valeur pour obtenir les observations. C'est utile pour un
435 algorithme itératif/incrémental.
436 - asVector : entrée des données, comme un vecteur compatible avec le
437 constructeur de numpy.matrix.
438 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
439 être rendue disponible au même titre que les variables de calcul
441 if asVector is not None:
442 self.__X.store( asVector )
444 self.__StoredInputs["Controls"] = self.__X
447 # -----------------------------------------------------------
448 def setAlgorithm(self, choice = None ):
450 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
451 d'assimilation. L'argument est un champ caractère se rapportant au nom
452 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
453 d'assimilation sur les arguments fixes.
456 raise ValueError("Error: algorithm choice has to be given")
457 if self.__algorithmName is not None:
458 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
459 daDirectory = "daAlgorithms"
461 # Recherche explicitement le fichier complet
462 # ------------------------------------------
464 for directory in sys.path:
465 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
466 module_path = os.path.abspath(os.path.join(directory, daDirectory))
467 if module_path is None:
468 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
470 # Importe le fichier complet comme un module
471 # ------------------------------------------
473 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
474 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
475 self.__algorithmName = str(choice)
476 sys.path = sys_path_tmp ; del sys_path_tmp
477 except ImportError, e:
478 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
480 # Instancie un objet du type élémentaire du fichier
481 # -------------------------------------------------
482 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
483 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
486 def setAlgorithmParameters(self, asDico=None):
488 Permet de définir les paramètres de l'algorithme, sous la forme d'un
491 if asDico is not None:
492 self.__Parameters.update( dict( asDico ) )
494 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
497 def getAlgorithmParameters(self, noDetails=True):
499 Renvoie la liste des paramètres requis selon l'algorithme
501 return self.__algorithm.getRequiredParameters(noDetails)
503 # -----------------------------------------------------------
504 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
506 Permet de sélectionner un diagnostic a effectuer.
509 raise ValueError("Error: diagnostic choice has to be given")
510 daDirectory = "daDiagnostics"
512 # Recherche explicitement le fichier complet
513 # ------------------------------------------
515 for directory in sys.path:
516 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
517 module_path = os.path.abspath(os.path.join(directory, daDirectory))
518 if module_path is None:
519 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
521 # Importe le fichier complet comme un module
522 # ------------------------------------------
524 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
525 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
526 sys.path = sys_path_tmp ; del sys_path_tmp
527 except ImportError, e:
528 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
530 # Instancie un objet du type élémentaire du fichier
531 # -------------------------------------------------
532 if self.__StoredInputs.has_key(name):
533 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
534 elif self.__StoredDiagnostics.has_key(name):
535 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
537 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
541 parameters = parameters )
544 # -----------------------------------------------------------
545 def shape_validate(self):
547 Validation de la correspondance correcte des tailles des variables et
548 des matrices s'il y en a.
550 if self.__Xb is None: __Xb_shape = (0,)
551 elif hasattr(self.__Xb,"shape"):
552 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
553 else: __Xb_shape = self.__Xb.shape()
554 else: raise TypeError("Xb has no attribute of shape: problem !")
556 if self.__Y is None: __Y_shape = (0,)
557 elif hasattr(self.__Y,"shape"):
558 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
559 else: __Y_shape = self.__Y.shape()
560 else: raise TypeError("Y has no attribute of shape: problem !")
562 if self.__B is None and self.__B_scalar is None:
564 elif self.__B is None and self.__B_scalar is not None:
565 __B_shape = (max(__Xb_shape),max(__Xb_shape))
566 elif hasattr(self.__B,"shape"):
567 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
568 else: __B_shape = self.__B.shape()
569 else: raise TypeError("B has no attribute of shape: problem !")
571 if self.__R is None and self.__R_scalar is None:
573 elif self.__R is None and self.__R_scalar is not None:
574 __R_shape = (max(__Y_shape),max(__Y_shape))
575 elif hasattr(self.__R,"shape"):
576 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
577 else: __R_shape = self.__R.shape()
578 else: raise TypeError("R has no attribute of shape: problem !")
580 if self.__Q is None: __Q_shape = (0,0)
581 elif hasattr(self.__Q,"shape"):
582 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
583 else: __Q_shape = self.__Q.shape()
584 else: raise TypeError("Q has no attribute of shape: problem !")
586 if len(self.__H) == 0: __H_shape = (0,0)
587 elif type(self.__H) is type({}): __H_shape = (0,0)
588 elif hasattr(self.__H["Direct"],"shape"):
589 if type(self.__H["Direct"].shape) is tuple: __H_shape = self.__H["Direct"].shape
590 else: __H_shape = self.__H["Direct"].shape()
591 else: raise TypeError("H has no attribute of shape: problem !")
593 if len(self.__M) == 0: __M_shape = (0,0)
594 elif type(self.__M) is type({}): __M_shape = (0,0)
595 elif hasattr(self.__M["Direct"],"shape"):
596 if type(self.__M["Direct"].shape) is tuple: __M_shape = self.__M["Direct"].shape
597 else: __M_shape = self.__M["Direct"].shape()
598 else: raise TypeError("M has no attribute of shape: problem !")
600 # Vérification des conditions
601 # ---------------------------
602 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
603 raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
604 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
605 raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
607 if not( min(__B_shape) == max(__B_shape) ):
608 raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
609 if not( min(__R_shape) == max(__R_shape) ):
610 raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
611 if not( min(__Q_shape) == max(__Q_shape) ):
612 raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
613 if not( min(__M_shape) == max(__M_shape) ):
614 raise ValueError("Shape characteristic of M is incorrect: \"%s\""%(__M_shape,))
616 if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[1] == max(__Xb_shape) ):
617 raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__H_shape,__Xb_shape))
618 if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[0] == max(__Y_shape) ):
619 raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__H_shape,__Y_shape))
620 if len(self.__H) > 0 and not(type(self.__H) is type({})) and len(self.__B) > 0 and not( __H_shape[1] == __B_shape[0] ):
621 raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__H_shape,__B_shape))
622 if len(self.__H) > 0 and not(type(self.__H) is type({})) and len(self.__R) > 0 and not( __H_shape[0] == __R_shape[1] ):
623 raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__H_shape,__R_shape))
625 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
626 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
627 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
628 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
629 for member in asPersistentVector:
630 self.__Xb.store( numpy.matrix( numpy.asarray(member).flatten(), numpy.float ).T )
631 __Xb_shape = min(__B_shape)
633 raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
635 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
636 raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
638 if self.__M is not None and len(self.__M) > 0 and not(type(self.__M) is type({})) and not( __M_shape[1] == max(__Xb_shape) ):
639 raise ValueError("Shape characteristic of M \"%s\" and X \"%s\" are incompatible"%(__M_shape,__Xb_shape))
643 # -----------------------------------------------------------
646 Permet de lancer le calcul d'assimilation.
648 Le nom de la méthode à activer est toujours "run". Les paramètres en
649 arguments de la méthode sont fixés. En sortie, on obtient les résultats
650 dans la variable de type dictionnaire "StoredVariables", qui contient en
651 particulier des objets de Persistence pour les analyses, OMA...
653 self.shape_validate()
655 self.__algorithm.run(
663 Parameters = self.__Parameters,
667 # -----------------------------------------------------------
668 def get(self, key=None):
670 Renvoie les résultats disponibles après l'exécution de la méthode
671 d'assimilation, ou les diagnostics disponibles. Attention, quand un
672 diagnostic porte le même nom qu'une variable stockée, c'est la variable
673 stockée qui est renvoyée, et le diagnostic est inatteignable.
676 if self.__algorithm.has_key(key):
677 return self.__algorithm.get( key )
678 elif self.__StoredInputs.has_key(key):
679 return self.__StoredInputs[key]
680 elif self.__StoredDiagnostics.has_key(key):
681 return self.__StoredDiagnostics[key]
683 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
685 allvariables = self.__algorithm.get()
686 allvariables.update( self.__StoredDiagnostics )
687 allvariables.update( self.__StoredInputs )
690 def get_available_variables(self):
692 Renvoie les variables potentiellement utilisables pour l'étude,
693 initialement stockées comme données d'entrées ou dans les algorithmes,
694 identifiés par les chaînes de caractères. L'algorithme doit avoir été
695 préalablement choisi sinon la méthode renvoie "None".
697 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
701 if len( self.__algorithm.keys()) > 0:
702 variables.extend( self.__algorithm.get().keys() )
703 if len( self.__StoredInputs.keys() ) > 0:
704 variables.extend( self.__StoredInputs.keys() )
708 def get_available_algorithms(self):
710 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
711 par les chaînes de caractères.
714 for directory in sys.path:
715 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
716 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
717 root, ext = os.path.splitext(fname)
718 if ext == '.py' and root != '__init__':
723 def get_available_diagnostics(self):
725 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
726 par les chaînes de caractères.
729 for directory in sys.path:
730 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
731 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
732 root, ext = os.path.splitext(fname)
733 if ext == '.py' and root != '__init__':
738 # -----------------------------------------------------------
739 def get_algorithms_main_path(self):
741 Renvoie le chemin pour le répertoire principal contenant les algorithmes
742 dans un sous-répertoire "daAlgorithms"
746 def add_algorithms_path(self, asPath=None):
748 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
749 se trouve un sous-répertoire "daAlgorithms"
751 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
752 pas indispensable de le rajouter ici.
754 if not os.path.isdir(asPath):
755 raise ValueError("The given "+asPath+" argument must exist as a directory")
756 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
757 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
758 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
759 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
760 sys.path.insert(0, os.path.abspath(asPath))
761 sys.path = list(set(sys.path)) # Conserve en unique exemplaire chaque chemin
764 def get_diagnostics_main_path(self):
766 Renvoie le chemin pour le répertoire principal contenant les diagnostics
767 dans un sous-répertoire "daDiagnostics"
771 def add_diagnostics_path(self, asPath=None):
773 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
774 se trouve un sous-répertoire "daDiagnostics"
776 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
777 pas indispensable de le rajouter ici.
779 if not os.path.isdir(asPath):
780 raise ValueError("The given "+asPath+" argument must exist as a directory")
781 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
782 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
783 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
784 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
785 sys.path.insert(0, os.path.abspath(asPath))
786 sys.path = list(set(sys.path)) # Conserve en unique exemplaire chaque chemin
789 # -----------------------------------------------------------
790 def setDataObserver(self,
793 HookParameters = None,
797 Permet d'associer un observer à une ou des variables nommées gérées en
798 interne, activable selon des règles définies dans le Scheduler. A chaque
799 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
800 les arguments (variable persistante VariableName, paramètres HookParameters).
803 if type( self.__algorithm ) is dict:
804 raise ValueError("No observer can be build before choosing an algorithm.")
806 # Vérification du nom de variable et typage
807 # -----------------------------------------
808 if type( VariableName ) is str:
809 VariableNames = [VariableName,]
810 elif type( VariableName ) is list:
811 VariableNames = map( str, VariableName )
813 raise ValueError("The observer requires a name or a list of names of variables.")
815 # Association interne de l'observer à la variable
816 # -----------------------------------------------
817 for n in VariableNames:
818 if not self.__algorithm.has_key( n ):
819 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
821 self.__algorithm.StoredVariables[ n ].setDataObserver(
822 Scheduler = Scheduler,
823 HookFunction = HookFunction,
824 HookParameters = HookParameters,
827 def removeDataObserver(self,
832 Permet de retirer un observer à une ou des variable nommée.
835 if type( self.__algorithm ) is dict:
836 raise ValueError("No observer can be removed before choosing an algorithm.")
838 # Vérification du nom de variable et typage
839 # -----------------------------------------
840 if type( VariableName ) is str:
841 VariableNames = [VariableName,]
842 elif type( VariableName ) is list:
843 VariableNames = map( str, VariableName )
845 raise ValueError("The observer requires a name or a list of names of variables.")
847 # Association interne de l'observer à la variable
848 # -----------------------------------------------
849 for n in VariableNames:
850 if not self.__algorithm.has_key( n ):
851 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
853 self.__algorithm.StoredVariables[ n ].removeDataObserver(
854 HookFunction = HookFunction,
857 # -----------------------------------------------------------
858 def setDebug(self, level=10):
860 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
861 appel pour changer le niveau de verbosité, avec :
862 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
865 log = logging.getLogger()
866 log.setLevel( level )
868 def unsetDebug(self):
870 Remet le logger au niveau par défaut
873 log = logging.getLogger()
874 log.setLevel( logging.WARNING )
876 def prepare_to_pickle(self):
877 self.__algorithmFile = None
878 self.__diagnosticFile = None
882 # ==============================================================================
883 if __name__ == "__main__":
884 print '\n AUTODIAGNOSTIC \n'
886 ADD = AssimilationStudy("Ma premiere etude BLUE")
888 ADD.setBackground (asVector = [0, 1, 2])
889 ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
890 ADD.setObservation (asVector = [0.5, 1.5, 2.5])
891 ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
892 ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
894 ADD.setAlgorithm(choice="Blue")
898 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
899 print "Ebauche :", [0, 1, 2]
900 print "Observation :", [0.5, 1.5, 2.5]
901 print "Demi-somme :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
902 print " qui doit être identique à :"
903 print "Analyse résultante :", ADD.get("Analysis").valueserie(0)
904 print "Innovation :", ADD.get("Innovation").valueserie(0)
907 print "Algorithmes disponibles.......................:", ADD.get_available_algorithms()
908 # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path()
909 print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics()
910 # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path()
911 print "Variables disponibles.........................:", ADD.get_available_variables()
914 print "Paramètres requis par algorithme :"
915 for algo in ADD.get_available_algorithms():
916 tmpADD = AssimilationStudy("Un algorithme")
917 tmpADD.setAlgorithm(choice=algo)
918 print " %25s : %s"%(algo,tmpADD.getAlgorithmParameters())
922 ADD.setDiagnostic("RMS", "Ma RMS")
924 liste = ADD.get().keys()
926 print "Variables et diagnostics nommés disponibles...:", liste
929 print "Exemple de mise en place d'un observeur :"
930 def obs(var=None,info=None):
931 print " ---> Mise en oeuvre de l'observer"
932 print " var =",var.valueserie(-1)
934 ADD.setDataObserver( 'Analysis', HookFunction=obs, Scheduler = [2, 4], HookParameters = "Second observer")
935 # Attention, il faut décaler le stockage de 1 pour suivre le pas interne
936 # car le pas 0 correspond à l'analyse ci-dessus.
939 print "Action sur la variable observée, étape :",i
940 ADD.get('Analysis').store( [i, i, i] )
943 print "Mise en debug et hors debug"
944 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
948 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
950 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()