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
36 from PlatformInfo import uniq
38 # ==============================================================================
39 class AssimilationStudy:
41 Cette classe sert d'interface pour l'utilisation de l'assimilation de
42 données. Elle contient les méthodes ou accesseurs nécessaires à la
43 construction d'un calcul d'assimilation.
45 def __init__(self, name=""):
47 Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
48 élémentaire. Ces variables sont ensuite disponibles pour implémenter un
49 algorithme élémentaire particulier.
51 Background............: vecteur Xb
52 Observation...........: vecteur Y (potentiellement temporel)
54 State.................: vecteur d'état dont une partie est le vecteur de
55 contrôle. Cette information n'est utile que si l'on veut faire des
56 calculs sur l'état complet, mais elle n'est pas indispensable pour
58 Control...............: vecteur X contenant toutes les variables de
59 contrôle, i.e. les paramètres ou l'état dont on veut estimer la
60 valeur pour obtenir les observations
61 ObservationOperator...: opérateur d'observation H
63 Les observations présentent une erreur dont la matrice de covariance est
64 R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
67 self.__name = str(name)
76 self.__X = Persistence.OneVector()
77 self.__Parameters = {}
78 self.__StoredDiagnostics = {}
79 self.__StoredInputs = {}
81 self.__B_scalar = None
82 self.__R_scalar = None
83 self.__Q_scalar = None
84 self.__Parameters["B_scalar"] = None
85 self.__Parameters["R_scalar"] = None
86 self.__Parameters["Q_scalar"] = None
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")
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 if asEyeByScalar is not None:
156 self.__B_scalar = float(asEyeByScalar)
158 elif asEyeByVector is not None:
159 self.__B_scalar = None
160 self.__B = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
161 elif asCovariance is not None:
162 self.__B_scalar = None
163 self.__B = numpy.matrix( asCovariance, numpy.float )
165 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")
167 self.__Parameters["B_scalar"] = self.__B_scalar
169 self.__StoredInputs["BackgroundError"] = self.__B
172 # -----------------------------------------------------------
173 def setObservation(self,
175 asPersistentVector = None,
180 Permet de définir les observations :
181 - asVector : entrée des données, comme un vecteur compatible avec le
182 constructeur de numpy.matrix
183 - asPersistentVector : entrée des données, comme un vecteur de type
184 persistent contruit avec la classe ad-hoc "Persistence"
185 - Scheduler est le contrôle temporel des données disponibles
186 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
187 être rendue disponible au même titre que les variables de calcul
189 if asVector is not None:
190 if type( asVector ) is type( numpy.matrix([]) ):
191 self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
193 self.__Y = numpy.matrix( asVector, numpy.float ).T
194 elif asPersistentVector is not None:
195 if type( asPersistentVector ) is list or type( asPersistentVector ) is tuple:
196 from Persistence import OneVector
197 self.__Y = OneVector("Observation", basetype=numpy.array)
198 for y in asPersistentVector:
201 self.__Y = asPersistentVector
203 raise ValueError("Error: improperly defined observations")
205 self.__StoredInputs["Observation"] = self.__Y
208 def setObservationError(self,
210 asEyeByScalar = None,
211 asEyeByVector = None,
215 Permet de définir la covariance des erreurs d'observations :
216 - asCovariance : entrée des données, comme une matrice compatible avec
217 le constructeur de numpy.matrix
218 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
219 multiplicatif d'une matrice de corrélation identité, aucune matrice
220 n'étant donc explicitement à donner
221 - asEyeByVector : entrée des données comme un seul vecteur de variance,
222 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
223 n'étant donc explicitement à donner
224 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
225 être rendue disponible au même titre que les variables de calcul
227 if asEyeByScalar is not None:
228 self.__R_scalar = float(asEyeByScalar)
230 elif asEyeByVector is not None:
231 self.__R_scalar = None
232 self.__R = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
233 elif asCovariance is not None:
234 self.__R_scalar = None
235 self.__R = numpy.matrix( asCovariance, numpy.float )
237 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")
239 self.__Parameters["R_scalar"] = self.__R_scalar
241 self.__StoredInputs["ObservationError"] = self.__R
244 def setObservationOperator(self,
245 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
246 "useApproximatedDerivatives":False,
247 "withCenteredDF" :False,
248 "withIncrement" :0.01,
256 Permet de définir un opérateur d'observation H. L'ordre de priorité des
257 définitions et leur sens sont les suivants :
258 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
259 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
260 "Direct" n'est pas définie, on prend la fonction "Tangent".
261 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
262 des opérateurs tangents et adjoints. On utilise par défaut des
263 différences finies non centrées ou centrées (si "withCenteredDF" est
264 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
265 du point courant ou sur le point fixe "withdX".
266 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
267 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
268 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
269 matrice fournie doit être sous une forme compatible avec le
270 constructeur de numpy.matrix.
271 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
272 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
273 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
274 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
275 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
276 être rendue disponible au même titre que les variables de calcul
278 if (type(asFunction) is type({})) and \
279 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
280 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
281 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
282 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
283 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
284 from daNumerics.ApproximatedDerivatives import FDApproximation
285 FDA = FDApproximation(
286 FunctionH = asFunction["Direct"],
287 centeredDF = asFunction["withCenteredDF"],
288 increment = asFunction["withIncrement"],
289 dX = asFunction["withdX"] )
290 self.__H["Direct"] = Operator( fromMethod = FDA.FunctionH )
291 self.__H["Tangent"] = Operator( fromMethod = FDA.TangentH )
292 self.__H["Adjoint"] = Operator( fromMethod = FDA.AdjointH )
293 elif (type(asFunction) is type({})) and \
294 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
295 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
296 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
297 self.__H["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
299 self.__H["Direct"] = Operator( fromMethod = asFunction["Direct"] )
300 self.__H["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
301 self.__H["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
302 elif asMatrix is not None:
303 matrice = numpy.matrix( asMatrix, numpy.float )
304 self.__H["Direct"] = Operator( fromMatrix = matrice )
305 self.__H["Tangent"] = Operator( fromMatrix = matrice )
306 self.__H["Adjoint"] = Operator( fromMatrix = matrice.T )
308 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
310 if appliedToX is not None:
311 self.__H["AppliedToX"] = {}
312 if type(appliedToX) is not dict:
313 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
314 for key in appliedToX.keys():
315 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
316 # Pour le cas où l'on a une vraie matrice
317 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
318 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
319 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
320 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
322 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
324 self.__H["AppliedToX"] = None
327 self.__StoredInputs["ObservationOperator"] = self.__H
330 # -----------------------------------------------------------
331 def setEvolutionModel(self,
332 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
333 "useApproximatedDerivatives":False,
334 "withCenteredDF" :False,
335 "withIncrement" :0.01,
343 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
344 définitions et leur sens sont les suivants :
345 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
346 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
347 "Direct" n'est pas définie, on prend la fonction "Tangent".
348 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
349 des opérateurs tangents et adjoints. On utilise par défaut des
350 différences finies non centrées ou centrées (si "withCenteredDF" est
351 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
352 du point courant ou sur le point fixe "withdX".
353 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
354 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
355 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
356 matrice fournie doit être sous une forme compatible avec le
357 constructeur de numpy.matrix.
358 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
359 être rendue disponible au même titre que les variables de calcul
361 if (type(asFunction) is type({})) and \
362 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
363 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
364 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
365 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
366 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
367 from daNumerics.ApproximatedDerivatives import FDApproximation
368 FDA = FDApproximation(
369 FunctionH = asFunction["Direct"],
370 centeredDF = asFunction["withCenteredDF"],
371 increment = asFunction["withIncrement"],
372 dX = asFunction["withdX"] )
373 self.__M["Direct"] = Operator( fromMethod = FDA.FunctionH )
374 self.__M["Tangent"] = Operator( fromMethod = FDA.TangentH )
375 self.__M["Adjoint"] = Operator( fromMethod = FDA.AdjointH )
376 elif (type(asFunction) is type({})) and \
377 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
378 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
379 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
380 self.__M["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
382 self.__M["Direct"] = Operator( fromMethod = asFunction["Direct"] )
383 self.__M["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
384 self.__M["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
385 elif asMatrix is not None:
386 matrice = numpy.matrix( asMatrix, numpy.float )
387 self.__M["Direct"] = Operator( fromMatrix = matrice )
388 self.__M["Tangent"] = Operator( fromMatrix = matrice )
389 self.__M["Adjoint"] = Operator( fromMatrix = matrice.T )
391 raise ValueError("Improperly defined evolution operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
394 self.__StoredInputs["EvolutionModel"] = self.__M
397 def setEvolutionError(self,
399 asEyeByScalar = None,
400 asEyeByVector = None,
404 Permet de définir la covariance des erreurs de modèle :
405 - asCovariance : entrée des données, comme une matrice compatible avec
406 le constructeur de numpy.matrix
407 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
408 multiplicatif d'une matrice de corrélation identité, aucune matrice
409 n'étant donc explicitement à donner
410 - asEyeByVector : entrée des données comme un seul vecteur de variance,
411 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
412 n'étant donc explicitement à donner
413 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
414 être rendue disponible au même titre que les variables de calcul
416 if asEyeByScalar is not None:
417 self.__Q_scalar = float(asEyeByScalar)
419 elif asEyeByVector is not None:
420 self.__Q_scalar = None
421 self.__Q = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
422 elif asCovariance is not None:
423 self.__Q_scalar = None
424 self.__Q = numpy.matrix( asCovariance, numpy.float )
426 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")
428 self.__Parameters["Q_scalar"] = self.__Q_scalar
430 self.__StoredInputs["EvolutionError"] = self.__Q
433 # -----------------------------------------------------------
434 def setControls (self,
439 Permet de définir la valeur initiale du vecteur X contenant toutes les
440 variables de contrôle, i.e. les paramètres ou l'état dont on veut
441 estimer la valeur pour obtenir les observations. C'est utile pour un
442 algorithme itératif/incrémental.
443 - asVector : entrée des données, comme un vecteur compatible avec le
444 constructeur de numpy.matrix.
445 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
446 être rendue disponible au même titre que les variables de calcul
448 if asVector is not None:
449 self.__X.store( asVector )
451 self.__StoredInputs["Controls"] = self.__X
454 # -----------------------------------------------------------
455 def setAlgorithm(self, choice = None ):
457 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
458 d'assimilation. L'argument est un champ caractère se rapportant au nom
459 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
460 d'assimilation sur les arguments fixes.
463 raise ValueError("Error: algorithm choice has to be given")
464 if self.__algorithmName is not None:
465 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
466 daDirectory = "daAlgorithms"
468 # Recherche explicitement le fichier complet
469 # ------------------------------------------
471 for directory in sys.path:
472 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
473 module_path = os.path.abspath(os.path.join(directory, daDirectory))
474 if module_path is None:
475 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
477 # Importe le fichier complet comme un module
478 # ------------------------------------------
480 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
481 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
482 self.__algorithmName = str(choice)
483 sys.path = sys_path_tmp ; del sys_path_tmp
484 except ImportError, e:
485 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
487 # Instancie un objet du type élémentaire du fichier
488 # -------------------------------------------------
489 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
490 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
493 def setAlgorithmParameters(self, asDico=None):
495 Permet de définir les paramètres de l'algorithme, sous la forme d'un
498 if asDico is not None:
499 self.__Parameters.update( dict( asDico ) )
501 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
504 def getAlgorithmParameters(self, noDetails=True):
506 Renvoie la liste des paramètres requis selon l'algorithme
508 return self.__algorithm.getRequiredParameters(noDetails)
510 # -----------------------------------------------------------
511 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
513 Permet de sélectionner un diagnostic a effectuer.
516 raise ValueError("Error: diagnostic choice has to be given")
517 daDirectory = "daDiagnostics"
519 # Recherche explicitement le fichier complet
520 # ------------------------------------------
522 for directory in sys.path:
523 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
524 module_path = os.path.abspath(os.path.join(directory, daDirectory))
525 if module_path is None:
526 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
528 # Importe le fichier complet comme un module
529 # ------------------------------------------
531 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
532 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
533 sys.path = sys_path_tmp ; del sys_path_tmp
534 except ImportError, e:
535 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
537 # Instancie un objet du type élémentaire du fichier
538 # -------------------------------------------------
539 if self.__StoredInputs.has_key(name):
540 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
541 elif self.__StoredDiagnostics.has_key(name):
542 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
544 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
548 parameters = parameters )
551 # -----------------------------------------------------------
552 def shape_validate(self):
554 Validation de la correspondance correcte des tailles des variables et
555 des matrices s'il y en a.
557 if self.__Xb is None: __Xb_shape = (0,)
558 elif hasattr(self.__Xb,"shape"):
559 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
560 else: __Xb_shape = self.__Xb.shape()
561 else: raise TypeError("Xb has no attribute of shape: problem !")
563 if self.__Y is None: __Y_shape = (0,)
564 elif hasattr(self.__Y,"shape"):
565 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
566 else: __Y_shape = self.__Y.shape()
567 else: raise TypeError("Y has no attribute of shape: problem !")
569 if self.__B is None and self.__B_scalar is None:
571 elif self.__B is None and self.__B_scalar is not None:
572 __B_shape = (max(__Xb_shape),max(__Xb_shape))
573 elif hasattr(self.__B,"shape"):
574 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
575 else: __B_shape = self.__B.shape()
576 else: raise TypeError("B has no attribute of shape: problem !")
578 if self.__R is None and self.__R_scalar is None:
580 elif self.__R is None and self.__R_scalar is not None:
581 __R_shape = (max(__Y_shape),max(__Y_shape))
582 elif hasattr(self.__R,"shape"):
583 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
584 else: __R_shape = self.__R.shape()
585 else: raise TypeError("R has no attribute of shape: problem !")
587 if self.__Q is None: __Q_shape = (0,0)
588 elif hasattr(self.__Q,"shape"):
589 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
590 else: __Q_shape = self.__Q.shape()
591 else: raise TypeError("Q has no attribute of shape: problem !")
593 if len(self.__H) == 0: __H_shape = (0,0)
594 elif type(self.__H) is type({}): __H_shape = (0,0)
595 elif hasattr(self.__H["Direct"],"shape"):
596 if type(self.__H["Direct"].shape) is tuple: __H_shape = self.__H["Direct"].shape
597 else: __H_shape = self.__H["Direct"].shape()
598 else: raise TypeError("H has no attribute of shape: problem !")
600 if len(self.__M) == 0: __M_shape = (0,0)
601 elif type(self.__M) is type({}): __M_shape = (0,0)
602 elif hasattr(self.__M["Direct"],"shape"):
603 if type(self.__M["Direct"].shape) is tuple: __M_shape = self.__M["Direct"].shape
604 else: __M_shape = self.__M["Direct"].shape()
605 else: raise TypeError("M has no attribute of shape: problem !")
607 # Vérification des conditions
608 # ---------------------------
609 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
610 raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
611 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
612 raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
614 if not( min(__B_shape) == max(__B_shape) ):
615 raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
616 if not( min(__R_shape) == max(__R_shape) ):
617 raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
618 if not( min(__Q_shape) == max(__Q_shape) ):
619 raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
620 if not( min(__M_shape) == max(__M_shape) ):
621 raise ValueError("Shape characteristic of M is incorrect: \"%s\""%(__M_shape,))
623 if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[1] == max(__Xb_shape) ):
624 raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__H_shape,__Xb_shape))
625 if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[0] == max(__Y_shape) ):
626 raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__H_shape,__Y_shape))
627 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] ):
628 raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__H_shape,__B_shape))
629 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] ):
630 raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__H_shape,__R_shape))
632 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
633 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
634 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
635 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
636 for member in asPersistentVector:
637 self.__Xb.store( numpy.matrix( numpy.asarray(member).flatten(), numpy.float ).T )
638 __Xb_shape = min(__B_shape)
640 raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
642 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
643 raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
645 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) ):
646 raise ValueError("Shape characteristic of M \"%s\" and X \"%s\" are incompatible"%(__M_shape,__Xb_shape))
650 # -----------------------------------------------------------
653 Permet de lancer le calcul d'assimilation.
655 Le nom de la méthode à activer est toujours "run". Les paramètres en
656 arguments de la méthode sont fixés. En sortie, on obtient les résultats
657 dans la variable de type dictionnaire "StoredVariables", qui contient en
658 particulier des objets de Persistence pour les analyses, OMA...
660 self.shape_validate()
662 self.__algorithm.run(
670 Parameters = self.__Parameters,
674 # -----------------------------------------------------------
675 def get(self, key=None):
677 Renvoie les résultats disponibles après l'exécution de la méthode
678 d'assimilation, ou les diagnostics disponibles. Attention, quand un
679 diagnostic porte le même nom qu'une variable stockée, c'est la variable
680 stockée qui est renvoyée, et le diagnostic est inatteignable.
683 if self.__algorithm.has_key(key):
684 return self.__algorithm.get( key )
685 elif self.__StoredInputs.has_key(key):
686 return self.__StoredInputs[key]
687 elif self.__StoredDiagnostics.has_key(key):
688 return self.__StoredDiagnostics[key]
690 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
692 allvariables = self.__algorithm.get()
693 allvariables.update( self.__StoredDiagnostics )
694 allvariables.update( self.__StoredInputs )
697 def get_available_variables(self):
699 Renvoie les variables potentiellement utilisables pour l'étude,
700 initialement stockées comme données d'entrées ou dans les algorithmes,
701 identifiés par les chaînes de caractères. L'algorithme doit avoir été
702 préalablement choisi sinon la méthode renvoie "None".
704 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
708 if len( self.__algorithm.keys()) > 0:
709 variables.extend( self.__algorithm.get().keys() )
710 if len( self.__StoredInputs.keys() ) > 0:
711 variables.extend( self.__StoredInputs.keys() )
715 def get_available_algorithms(self):
717 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
718 par les chaînes de caractères.
721 for directory in sys.path:
722 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
723 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
724 root, ext = os.path.splitext(fname)
725 if ext == '.py' and root != '__init__':
730 def get_available_diagnostics(self):
732 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
733 par les chaînes de caractères.
736 for directory in sys.path:
737 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
738 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
739 root, ext = os.path.splitext(fname)
740 if ext == '.py' and root != '__init__':
745 # -----------------------------------------------------------
746 def get_algorithms_main_path(self):
748 Renvoie le chemin pour le répertoire principal contenant les algorithmes
749 dans un sous-répertoire "daAlgorithms"
753 def add_algorithms_path(self, asPath=None):
755 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
756 se trouve un sous-répertoire "daAlgorithms"
758 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
759 pas indispensable de le rajouter ici.
761 if not os.path.isdir(asPath):
762 raise ValueError("The given "+asPath+" argument must exist as a directory")
763 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
764 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
765 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
766 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
767 sys.path.insert(0, os.path.abspath(asPath))
768 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
771 def get_diagnostics_main_path(self):
773 Renvoie le chemin pour le répertoire principal contenant les diagnostics
774 dans un sous-répertoire "daDiagnostics"
778 def add_diagnostics_path(self, asPath=None):
780 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
781 se trouve un sous-répertoire "daDiagnostics"
783 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
784 pas indispensable de le rajouter ici.
786 if not os.path.isdir(asPath):
787 raise ValueError("The given "+asPath+" argument must exist as a directory")
788 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
789 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
790 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
791 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
792 sys.path.insert(0, os.path.abspath(asPath))
793 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
796 # -----------------------------------------------------------
797 def setDataObserver(self,
800 HookParameters = None,
804 Permet d'associer un observer à une ou des variables nommées gérées en
805 interne, activable selon des règles définies dans le Scheduler. A chaque
806 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
807 les arguments (variable persistante VariableName, paramètres HookParameters).
810 if type( self.__algorithm ) is dict:
811 raise ValueError("No observer can be build before choosing an algorithm.")
813 # Vérification du nom de variable et typage
814 # -----------------------------------------
815 if type( VariableName ) is str:
816 VariableNames = [VariableName,]
817 elif type( VariableName ) is list:
818 VariableNames = map( str, VariableName )
820 raise ValueError("The observer requires a name or a list of names of variables.")
822 # Association interne de l'observer à la variable
823 # -----------------------------------------------
824 for n in VariableNames:
825 if not self.__algorithm.has_key( n ):
826 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
828 self.__algorithm.StoredVariables[ n ].setDataObserver(
829 Scheduler = Scheduler,
830 HookFunction = HookFunction,
831 HookParameters = HookParameters,
834 def removeDataObserver(self,
839 Permet de retirer un observer à une ou des variable nommée.
842 if type( self.__algorithm ) is dict:
843 raise ValueError("No observer can be removed before choosing an algorithm.")
845 # Vérification du nom de variable et typage
846 # -----------------------------------------
847 if type( VariableName ) is str:
848 VariableNames = [VariableName,]
849 elif type( VariableName ) is list:
850 VariableNames = map( str, VariableName )
852 raise ValueError("The observer requires a name or a list of names of variables.")
854 # Association interne de l'observer à la variable
855 # -----------------------------------------------
856 for n in VariableNames:
857 if not self.__algorithm.has_key( n ):
858 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
860 self.__algorithm.StoredVariables[ n ].removeDataObserver(
861 HookFunction = HookFunction,
864 # -----------------------------------------------------------
865 def setDebug(self, level=10):
867 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
868 appel pour changer le niveau de verbosité, avec :
869 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
872 log = logging.getLogger()
873 log.setLevel( level )
875 def unsetDebug(self):
877 Remet le logger au niveau par défaut
880 log = logging.getLogger()
881 log.setLevel( logging.WARNING )
883 def prepare_to_pickle(self):
884 self.__algorithmFile = None
885 self.__diagnosticFile = None
889 # ==============================================================================
890 if __name__ == "__main__":
891 print '\n AUTODIAGNOSTIC \n'
893 ADD = AssimilationStudy("Ma premiere etude BLUE")
895 ADD.setBackground (asVector = [0, 1, 2])
896 ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
897 ADD.setObservation (asVector = [0.5, 1.5, 2.5])
898 ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
899 ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
901 ADD.setAlgorithm(choice="Blue")
905 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
906 print "Ebauche :", [0, 1, 2]
907 print "Observation :", [0.5, 1.5, 2.5]
908 print "Demi-somme :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
909 print " qui doit être identique à :"
910 print "Analyse résultante :", ADD.get("Analysis").valueserie(0)
911 print "Innovation :", ADD.get("Innovation").valueserie(0)
914 print "Algorithmes disponibles.......................:", ADD.get_available_algorithms()
915 # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path()
916 print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics()
917 # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path()
918 print "Variables disponibles.........................:", ADD.get_available_variables()
921 print "Paramètres requis par algorithme :"
922 for algo in ADD.get_available_algorithms():
923 tmpADD = AssimilationStudy("Un algorithme")
924 tmpADD.setAlgorithm(choice=algo)
925 print " %25s : %s"%(algo,tmpADD.getAlgorithmParameters())
929 ADD.setDiagnostic("RMS", "Ma RMS")
931 liste = ADD.get().keys()
933 print "Variables et diagnostics nommés disponibles...:", liste
936 print "Exemple de mise en place d'un observeur :"
937 def obs(var=None,info=None):
938 print " ---> Mise en oeuvre de l'observer"
939 print " var =",var.valueserie(-1)
941 ADD.setDataObserver( 'Analysis', HookFunction=obs, Scheduler = [2, 4], HookParameters = "Second observer")
942 # Attention, il faut décaler le stockage de 1 pour suivre le pas interne
943 # car le pas 0 correspond à l'analyse ci-dessus.
946 print "Action sur la variable observée, étape :",i
947 ADD.get('Analysis').store( [i, i, i] )
950 print "Mise en debug et hors debug"
951 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
955 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
957 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()