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 self.__Y = Persistence.OneVector("Observation", basetype=numpy.array)
197 for y in asPersistentVector:
200 self.__Y = asPersistentVector
202 raise ValueError("Error: improperly defined observations")
204 self.__StoredInputs["Observation"] = self.__Y
207 def setObservationError(self,
209 asEyeByScalar = None,
210 asEyeByVector = None,
214 Permet de définir la covariance des erreurs d'observations :
215 - asCovariance : entrée des données, comme une matrice compatible avec
216 le constructeur de numpy.matrix
217 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
218 multiplicatif d'une matrice de corrélation identité, aucune matrice
219 n'étant donc explicitement à donner
220 - asEyeByVector : entrée des données comme un seul vecteur de variance,
221 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
222 n'étant donc explicitement à donner
223 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
224 être rendue disponible au même titre que les variables de calcul
226 if asEyeByScalar is not None:
227 self.__R_scalar = float(asEyeByScalar)
229 elif asEyeByVector is not None:
230 self.__R_scalar = None
231 self.__R = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
232 elif asCovariance is not None:
233 self.__R_scalar = None
234 self.__R = numpy.matrix( asCovariance, numpy.float )
236 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")
238 self.__Parameters["R_scalar"] = self.__R_scalar
240 self.__StoredInputs["ObservationError"] = self.__R
243 def setObservationOperator(self,
244 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
245 "useApproximatedDerivatives":False,
246 "withCenteredDF" :False,
247 "withIncrement" :0.01,
255 Permet de définir un opérateur d'observation H. L'ordre de priorité des
256 définitions et leur sens sont les suivants :
257 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
258 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
259 "Direct" n'est pas définie, on prend la fonction "Tangent".
260 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
261 des opérateurs tangents et adjoints. On utilise par défaut des
262 différences finies non centrées ou centrées (si "withCenteredDF" est
263 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
264 du point courant ou sur le point fixe "withdX".
265 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
266 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
267 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
268 matrice fournie doit être sous une forme compatible avec le
269 constructeur de numpy.matrix.
270 - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
271 X divers, l'opérateur par sa valeur appliquée à cet X particulier,
272 sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
273 L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
274 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
275 être rendue disponible au même titre que les variables de calcul
277 if (type(asFunction) is type({})) and \
278 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
279 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
280 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
281 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
282 if not asFunction.has_key("withdX"): asFunction["withdX"] = 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 self.__H["Direct"] = Operator( fromMethod = FDA.DirectOperator )
290 self.__H["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
291 self.__H["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
292 elif (type(asFunction) is type({})) and \
293 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
294 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
295 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
296 self.__H["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
298 self.__H["Direct"] = Operator( fromMethod = asFunction["Direct"] )
299 self.__H["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
300 self.__H["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
301 elif asMatrix is not None:
302 matrice = numpy.matrix( asMatrix, numpy.float )
303 self.__H["Direct"] = Operator( fromMatrix = matrice )
304 self.__H["Tangent"] = Operator( fromMatrix = matrice )
305 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 Function = asFunction["Direct"],
370 centeredDF = asFunction["withCenteredDF"],
371 increment = asFunction["withIncrement"],
372 dX = asFunction["withdX"] )
373 self.__M["Direct"] = Operator( fromMethod = FDA.DirectOperator )
374 self.__M["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
375 self.__M["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
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 )
392 raise ValueError("Improperly defined evolution operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
395 self.__StoredInputs["EvolutionModel"] = self.__M
398 def setEvolutionError(self,
400 asEyeByScalar = None,
401 asEyeByVector = None,
405 Permet de définir la covariance des erreurs de modèle :
406 - asCovariance : entrée des données, comme une matrice compatible avec
407 le constructeur de numpy.matrix
408 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
409 multiplicatif d'une matrice de corrélation identité, aucune matrice
410 n'étant donc explicitement à donner
411 - asEyeByVector : entrée des données comme un seul vecteur de variance,
412 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
413 n'étant donc explicitement à donner
414 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
415 être rendue disponible au même titre que les variables de calcul
417 if asEyeByScalar is not None:
418 self.__Q_scalar = float(asEyeByScalar)
420 elif asEyeByVector is not None:
421 self.__Q_scalar = None
422 self.__Q = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
423 elif asCovariance is not None:
424 self.__Q_scalar = None
425 self.__Q = numpy.matrix( asCovariance, numpy.float )
427 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")
429 self.__Parameters["Q_scalar"] = self.__Q_scalar
431 self.__StoredInputs["EvolutionError"] = self.__Q
434 # -----------------------------------------------------------
435 def setControls (self,
440 Permet de définir la valeur initiale du vecteur X contenant toutes les
441 variables de contrôle, i.e. les paramètres ou l'état dont on veut
442 estimer la valeur pour obtenir les observations. C'est utile pour un
443 algorithme itératif/incrémental.
444 - asVector : entrée des données, comme un vecteur compatible avec le
445 constructeur de numpy.matrix.
446 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
447 être rendue disponible au même titre que les variables de calcul
449 if asVector is not None:
450 self.__X.store( asVector )
452 self.__StoredInputs["Controls"] = self.__X
455 # -----------------------------------------------------------
456 def setAlgorithm(self, choice = None ):
458 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
459 d'assimilation. L'argument est un champ caractère se rapportant au nom
460 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
461 d'assimilation sur les arguments fixes.
464 raise ValueError("Error: algorithm choice has to be given")
465 if self.__algorithmName is not None:
466 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
467 daDirectory = "daAlgorithms"
469 # Recherche explicitement le fichier complet
470 # ------------------------------------------
472 for directory in sys.path:
473 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
474 module_path = os.path.abspath(os.path.join(directory, daDirectory))
475 if module_path is None:
476 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
478 # Importe le fichier complet comme un module
479 # ------------------------------------------
481 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
482 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
483 self.__algorithmName = str(choice)
484 sys.path = sys_path_tmp ; del sys_path_tmp
485 except ImportError, e:
486 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
488 # Instancie un objet du type élémentaire du fichier
489 # -------------------------------------------------
490 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
491 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
494 def setAlgorithmParameters(self, asDico=None):
496 Permet de définir les paramètres de l'algorithme, sous la forme d'un
499 if asDico is not None:
500 self.__Parameters.update( dict( asDico ) )
502 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
505 def getAlgorithmParameters(self, noDetails=True):
507 Renvoie la liste des paramètres requis selon l'algorithme
509 return self.__algorithm.getRequiredParameters(noDetails)
511 # -----------------------------------------------------------
512 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
514 Permet de sélectionner un diagnostic a effectuer.
517 raise ValueError("Error: diagnostic choice has to be given")
518 daDirectory = "daDiagnostics"
520 # Recherche explicitement le fichier complet
521 # ------------------------------------------
523 for directory in sys.path:
524 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
525 module_path = os.path.abspath(os.path.join(directory, daDirectory))
526 if module_path is None:
527 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
529 # Importe le fichier complet comme un module
530 # ------------------------------------------
532 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
533 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
534 sys.path = sys_path_tmp ; del sys_path_tmp
535 except ImportError, e:
536 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
538 # Instancie un objet du type élémentaire du fichier
539 # -------------------------------------------------
540 if self.__StoredInputs.has_key(name):
541 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
542 elif self.__StoredDiagnostics.has_key(name):
543 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
545 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
549 parameters = parameters )
552 # -----------------------------------------------------------
553 def shape_validate(self):
555 Validation de la correspondance correcte des tailles des variables et
556 des matrices s'il y en a.
558 if self.__Xb is None: __Xb_shape = (0,)
559 elif hasattr(self.__Xb,"shape"):
560 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
561 else: __Xb_shape = self.__Xb.shape()
562 else: raise TypeError("Xb has no attribute of shape: problem !")
564 if self.__Y is None: __Y_shape = (0,)
565 elif hasattr(self.__Y,"shape"):
566 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
567 else: __Y_shape = self.__Y.shape()
568 else: raise TypeError("Y has no attribute of shape: problem !")
570 if self.__B is None and self.__B_scalar is None:
572 elif self.__B is None and self.__B_scalar is not None:
573 __B_shape = (max(__Xb_shape),max(__Xb_shape))
574 elif hasattr(self.__B,"shape"):
575 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
576 else: __B_shape = self.__B.shape()
577 else: raise TypeError("B has no attribute of shape: problem !")
579 if self.__R is None and self.__R_scalar is None:
581 elif self.__R is None and self.__R_scalar is not None:
582 __R_shape = (max(__Y_shape),max(__Y_shape))
583 elif hasattr(self.__R,"shape"):
584 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
585 else: __R_shape = self.__R.shape()
586 else: raise TypeError("R has no attribute of shape: problem !")
588 if self.__Q is None: __Q_shape = (0,0)
589 elif hasattr(self.__Q,"shape"):
590 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
591 else: __Q_shape = self.__Q.shape()
592 else: raise TypeError("Q has no attribute of shape: problem !")
594 if len(self.__H) == 0: __H_shape = (0,0)
595 elif type(self.__H) is type({}): __H_shape = (0,0)
596 elif hasattr(self.__H["Direct"],"shape"):
597 if type(self.__H["Direct"].shape) is tuple: __H_shape = self.__H["Direct"].shape
598 else: __H_shape = self.__H["Direct"].shape()
599 else: raise TypeError("H has no attribute of shape: problem !")
601 if len(self.__M) == 0: __M_shape = (0,0)
602 elif type(self.__M) is type({}): __M_shape = (0,0)
603 elif hasattr(self.__M["Direct"],"shape"):
604 if type(self.__M["Direct"].shape) is tuple: __M_shape = self.__M["Direct"].shape
605 else: __M_shape = self.__M["Direct"].shape()
606 else: raise TypeError("M has no attribute of shape: problem !")
608 # Vérification des conditions
609 # ---------------------------
610 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
611 raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
612 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
613 raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
615 if not( min(__B_shape) == max(__B_shape) ):
616 raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
617 if not( min(__R_shape) == max(__R_shape) ):
618 raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
619 if not( min(__Q_shape) == max(__Q_shape) ):
620 raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
621 if not( min(__M_shape) == max(__M_shape) ):
622 raise ValueError("Shape characteristic of M is incorrect: \"%s\""%(__M_shape,))
624 if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[1] == max(__Xb_shape) ):
625 raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__H_shape,__Xb_shape))
626 if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[0] == max(__Y_shape) ):
627 raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__H_shape,__Y_shape))
628 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] ):
629 raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__H_shape,__B_shape))
630 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] ):
631 raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__H_shape,__R_shape))
633 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
634 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
635 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
636 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
637 for member in asPersistentVector:
638 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
639 __Xb_shape = min(__B_shape)
641 raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
643 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
644 raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
646 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) ):
647 raise ValueError("Shape characteristic of M \"%s\" and X \"%s\" are incompatible"%(__M_shape,__Xb_shape))
651 # -----------------------------------------------------------
654 Permet de lancer le calcul d'assimilation.
656 Le nom de la méthode à activer est toujours "run". Les paramètres en
657 arguments de la méthode sont fixés. En sortie, on obtient les résultats
658 dans la variable de type dictionnaire "StoredVariables", qui contient en
659 particulier des objets de Persistence pour les analyses, OMA...
661 self.shape_validate()
663 self.__algorithm.run(
671 Parameters = self.__Parameters,
675 # -----------------------------------------------------------
676 def get(self, key=None):
678 Renvoie les résultats disponibles après l'exécution de la méthode
679 d'assimilation, ou les diagnostics disponibles. Attention, quand un
680 diagnostic porte le même nom qu'une variable stockée, c'est la variable
681 stockée qui est renvoyée, et le diagnostic est inatteignable.
684 if self.__algorithm.has_key(key):
685 return self.__algorithm.get( key )
686 elif self.__StoredInputs.has_key(key):
687 return self.__StoredInputs[key]
688 elif self.__StoredDiagnostics.has_key(key):
689 return self.__StoredDiagnostics[key]
691 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
693 allvariables = self.__algorithm.get()
694 allvariables.update( self.__StoredDiagnostics )
695 allvariables.update( self.__StoredInputs )
698 def get_available_variables(self):
700 Renvoie les variables potentiellement utilisables pour l'étude,
701 initialement stockées comme données d'entrées ou dans les algorithmes,
702 identifiés par les chaînes de caractères. L'algorithme doit avoir été
703 préalablement choisi sinon la méthode renvoie "None".
705 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
709 if len( self.__algorithm.keys()) > 0:
710 variables.extend( self.__algorithm.get().keys() )
711 if len( self.__StoredInputs.keys() ) > 0:
712 variables.extend( self.__StoredInputs.keys() )
716 def get_available_algorithms(self):
718 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
719 par les chaînes de caractères.
722 for directory in sys.path:
723 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
724 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
725 root, ext = os.path.splitext(fname)
726 if ext == '.py' and root != '__init__':
731 def get_available_diagnostics(self):
733 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
734 par les chaînes de caractères.
737 for directory in sys.path:
738 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
739 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
740 root, ext = os.path.splitext(fname)
741 if ext == '.py' and root != '__init__':
746 # -----------------------------------------------------------
747 def get_algorithms_main_path(self):
749 Renvoie le chemin pour le répertoire principal contenant les algorithmes
750 dans un sous-répertoire "daAlgorithms"
754 def add_algorithms_path(self, asPath=None):
756 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
757 se trouve un sous-répertoire "daAlgorithms"
759 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
760 pas indispensable de le rajouter ici.
762 if not os.path.isdir(asPath):
763 raise ValueError("The given "+asPath+" argument must exist as a directory")
764 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
765 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
766 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
767 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
768 sys.path.insert(0, os.path.abspath(asPath))
769 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
772 def get_diagnostics_main_path(self):
774 Renvoie le chemin pour le répertoire principal contenant les diagnostics
775 dans un sous-répertoire "daDiagnostics"
779 def add_diagnostics_path(self, asPath=None):
781 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
782 se trouve un sous-répertoire "daDiagnostics"
784 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
785 pas indispensable de le rajouter ici.
787 if not os.path.isdir(asPath):
788 raise ValueError("The given "+asPath+" argument must exist as a directory")
789 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
790 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
791 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
792 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
793 sys.path.insert(0, os.path.abspath(asPath))
794 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
797 # -----------------------------------------------------------
798 def setDataObserver(self,
801 HookParameters = None,
805 Permet d'associer un observer à une ou des variables nommées gérées en
806 interne, activable selon des règles définies dans le Scheduler. A chaque
807 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
808 les arguments (variable persistante VariableName, paramètres HookParameters).
811 if type( self.__algorithm ) is dict:
812 raise ValueError("No observer can be build before choosing an algorithm.")
814 # Vérification du nom de variable et typage
815 # -----------------------------------------
816 if type( VariableName ) is str:
817 VariableNames = [VariableName,]
818 elif type( VariableName ) is list:
819 VariableNames = map( str, VariableName )
821 raise ValueError("The observer requires a name or a list of names of variables.")
823 # Association interne de l'observer à la variable
824 # -----------------------------------------------
825 for n in VariableNames:
826 if not self.__algorithm.has_key( n ):
827 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
829 self.__algorithm.StoredVariables[ n ].setDataObserver(
830 Scheduler = Scheduler,
831 HookFunction = HookFunction,
832 HookParameters = HookParameters,
835 def removeDataObserver(self,
840 Permet de retirer un observer à une ou des variable nommée.
843 if type( self.__algorithm ) is dict:
844 raise ValueError("No observer can be removed before choosing an algorithm.")
846 # Vérification du nom de variable et typage
847 # -----------------------------------------
848 if type( VariableName ) is str:
849 VariableNames = [VariableName,]
850 elif type( VariableName ) is list:
851 VariableNames = map( str, VariableName )
853 raise ValueError("The observer requires a name or a list of names of variables.")
855 # Association interne de l'observer à la variable
856 # -----------------------------------------------
857 for n in VariableNames:
858 if not self.__algorithm.has_key( n ):
859 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
861 self.__algorithm.StoredVariables[ n ].removeDataObserver(
862 HookFunction = HookFunction,
865 # -----------------------------------------------------------
866 def setDebug(self, level=10):
868 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
869 appel pour changer le niveau de verbosité, avec :
870 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
873 log = logging.getLogger()
874 log.setLevel( level )
876 def unsetDebug(self):
878 Remet le logger au niveau par défaut
881 log = logging.getLogger()
882 log.setLevel( logging.WARNING )
884 def prepare_to_pickle(self):
885 self.__algorithmFile = None
886 self.__diagnosticFile = None
890 # ==============================================================================
891 if __name__ == "__main__":
892 print '\n AUTODIAGNOSTIC \n'
894 ADD = AssimilationStudy("Ma premiere etude BLUE")
896 ADD.setBackground (asVector = [0, 1, 2])
897 ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
898 ADD.setObservation (asVector = [0.5, 1.5, 2.5])
899 ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
900 ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
902 ADD.setAlgorithm(choice="Blue")
906 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
907 print "Ebauche :", [0, 1, 2]
908 print "Observation :", [0.5, 1.5, 2.5]
909 print "Demi-somme :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
910 print " qui doit être identique à :"
911 print "Analyse résultante :", ADD.get("Analysis").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()