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 Function = asFunction["Direct"],
287 centeredDF = asFunction["withCenteredDF"],
288 increment = asFunction["withIncrement"],
289 dX = asFunction["withdX"] )
290 self.__H["Direct"] = Operator( fromMethod = FDA.DirectOperator )
291 self.__H["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
292 self.__H["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
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 )
309 raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
311 if appliedToX is not None:
312 self.__H["AppliedToX"] = {}
313 if type(appliedToX) is not dict:
314 raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
315 for key in appliedToX.keys():
316 if type( appliedToX[key] ) is type( numpy.matrix([]) ):
317 # Pour le cas où l'on a une vraie matrice
318 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
319 elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
320 # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
321 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
323 self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
325 self.__H["AppliedToX"] = None
328 self.__StoredInputs["ObservationOperator"] = self.__H
331 # -----------------------------------------------------------
332 def setEvolutionModel(self,
333 asFunction = {"Direct":None, "Tangent":None, "Adjoint":None,
334 "useApproximatedDerivatives":False,
335 "withCenteredDF" :False,
336 "withIncrement" :0.01,
344 Permet de définir un opérateur d'évolution M. L'ordre de priorité des
345 définitions et leur sens sont les suivants :
346 - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
347 alors on définit l'opérateur à l'aide de fonctions. Si la fonction
348 "Direct" n'est pas définie, on prend la fonction "Tangent".
349 Si "useApproximatedDerivatives" est vrai, on utilise une approximation
350 des opérateurs tangents et adjoints. On utilise par défaut des
351 différences finies non centrées ou centrées (si "withCenteredDF" est
352 vrai) avec un incrément multiplicatif "withIncrement" de 1% autour
353 du point courant ou sur le point fixe "withdX".
354 - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
355 None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
356 la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
357 matrice fournie doit être sous une forme compatible avec le
358 constructeur de numpy.matrix.
359 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
360 être rendue disponible au même titre que les variables de calcul
362 if (type(asFunction) is type({})) and \
363 asFunction.has_key("useApproximatedDerivatives") and bool(asFunction["useApproximatedDerivatives"]) and \
364 asFunction.has_key("Direct") and (asFunction["Direct"] is not None):
365 if not asFunction.has_key("withCenteredDF"): asFunction["withCenteredDF"] = False
366 if not asFunction.has_key("withIncrement"): asFunction["withIncrement"] = 0.01
367 if not asFunction.has_key("withdX"): asFunction["withdX"] = None
368 from daNumerics.ApproximatedDerivatives import FDApproximation
369 FDA = FDApproximation(
370 Function = asFunction["Direct"],
371 centeredDF = asFunction["withCenteredDF"],
372 increment = asFunction["withIncrement"],
373 dX = asFunction["withdX"] )
374 self.__M["Direct"] = Operator( fromMethod = FDA.DirectOperator )
375 self.__M["Tangent"] = Operator( fromMethod = FDA.TangentOperator )
376 self.__M["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator )
377 elif (type(asFunction) is type({})) and \
378 asFunction.has_key("Tangent") and asFunction.has_key("Adjoint") and \
379 (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
380 if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
381 self.__M["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
383 self.__M["Direct"] = Operator( fromMethod = asFunction["Direct"] )
384 self.__M["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
385 self.__M["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
386 elif asMatrix is not None:
387 matrice = numpy.matrix( asMatrix, numpy.float )
388 self.__M["Direct"] = Operator( fromMatrix = matrice )
389 self.__M["Tangent"] = Operator( fromMatrix = matrice )
390 self.__M["Adjoint"] = Operator( fromMatrix = matrice.T )
393 raise ValueError("Improperly defined evolution operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
396 self.__StoredInputs["EvolutionModel"] = self.__M
399 def setEvolutionError(self,
401 asEyeByScalar = None,
402 asEyeByVector = None,
406 Permet de définir la covariance des erreurs de modèle :
407 - asCovariance : entrée des données, comme une matrice compatible avec
408 le constructeur de numpy.matrix
409 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
410 multiplicatif d'une matrice de corrélation identité, aucune matrice
411 n'étant donc explicitement à donner
412 - asEyeByVector : entrée des données comme un seul vecteur de variance,
413 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
414 n'étant donc explicitement à donner
415 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
416 être rendue disponible au même titre que les variables de calcul
418 if asEyeByScalar is not None:
419 self.__Q_scalar = float(asEyeByScalar)
421 elif asEyeByVector is not None:
422 self.__Q_scalar = None
423 self.__Q = numpy.matrix( numpy.diagflat( asEyeByVector ), numpy.float )
424 elif asCovariance is not None:
425 self.__Q_scalar = None
426 self.__Q = numpy.matrix( asCovariance, numpy.float )
428 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")
430 self.__Parameters["Q_scalar"] = self.__Q_scalar
432 self.__StoredInputs["EvolutionError"] = self.__Q
435 # -----------------------------------------------------------
436 def setControls (self,
441 Permet de définir la valeur initiale du vecteur X contenant toutes les
442 variables de contrôle, i.e. les paramètres ou l'état dont on veut
443 estimer la valeur pour obtenir les observations. C'est utile pour un
444 algorithme itératif/incrémental.
445 - asVector : entrée des données, comme un vecteur compatible avec le
446 constructeur de numpy.matrix.
447 - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
448 être rendue disponible au même titre que les variables de calcul
450 if asVector is not None:
451 self.__X.store( asVector )
453 self.__StoredInputs["Controls"] = self.__X
456 # -----------------------------------------------------------
457 def setAlgorithm(self, choice = None ):
459 Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
460 d'assimilation. L'argument est un champ caractère se rapportant au nom
461 d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
462 d'assimilation sur les arguments fixes.
465 raise ValueError("Error: algorithm choice has to be given")
466 if self.__algorithmName is not None:
467 raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
468 daDirectory = "daAlgorithms"
470 # Recherche explicitement le fichier complet
471 # ------------------------------------------
473 for directory in sys.path:
474 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
475 module_path = os.path.abspath(os.path.join(directory, daDirectory))
476 if module_path is None:
477 raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
479 # Importe le fichier complet comme un module
480 # ------------------------------------------
482 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
483 self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
484 self.__algorithmName = str(choice)
485 sys.path = sys_path_tmp ; del sys_path_tmp
486 except ImportError, e:
487 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
489 # Instancie un objet du type élémentaire du fichier
490 # -------------------------------------------------
491 self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
492 self.__StoredInputs["AlgorithmName"] = self.__algorithmName
495 def setAlgorithmParameters(self, asDico=None):
497 Permet de définir les paramètres de l'algorithme, sous la forme d'un
500 if asDico is not None:
501 self.__Parameters.update( dict( asDico ) )
503 self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
506 def getAlgorithmParameters(self, noDetails=True):
508 Renvoie la liste des paramètres requis selon l'algorithme
510 return self.__algorithm.getRequiredParameters(noDetails)
512 # -----------------------------------------------------------
513 def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
515 Permet de sélectionner un diagnostic a effectuer.
518 raise ValueError("Error: diagnostic choice has to be given")
519 daDirectory = "daDiagnostics"
521 # Recherche explicitement le fichier complet
522 # ------------------------------------------
524 for directory in sys.path:
525 if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
526 module_path = os.path.abspath(os.path.join(directory, daDirectory))
527 if module_path is None:
528 raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
530 # Importe le fichier complet comme un module
531 # ------------------------------------------
533 sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
534 self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
535 sys.path = sys_path_tmp ; del sys_path_tmp
536 except ImportError, e:
537 raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
539 # Instancie un objet du type élémentaire du fichier
540 # -------------------------------------------------
541 if self.__StoredInputs.has_key(name):
542 raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
543 elif self.__StoredDiagnostics.has_key(name):
544 raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
546 self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
550 parameters = parameters )
553 # -----------------------------------------------------------
554 def shape_validate(self):
556 Validation de la correspondance correcte des tailles des variables et
557 des matrices s'il y en a.
559 if self.__Xb is None: __Xb_shape = (0,)
560 elif hasattr(self.__Xb,"shape"):
561 if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
562 else: __Xb_shape = self.__Xb.shape()
563 else: raise TypeError("Xb has no attribute of shape: problem !")
565 if self.__Y is None: __Y_shape = (0,)
566 elif hasattr(self.__Y,"shape"):
567 if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
568 else: __Y_shape = self.__Y.shape()
569 else: raise TypeError("Y has no attribute of shape: problem !")
571 if self.__B is None and self.__B_scalar is None:
573 elif self.__B is None and self.__B_scalar is not None:
574 __B_shape = (max(__Xb_shape),max(__Xb_shape))
575 elif hasattr(self.__B,"shape"):
576 if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
577 else: __B_shape = self.__B.shape()
578 else: raise TypeError("B has no attribute of shape: problem !")
580 if self.__R is None and self.__R_scalar is None:
582 elif self.__R is None and self.__R_scalar is not None:
583 __R_shape = (max(__Y_shape),max(__Y_shape))
584 elif hasattr(self.__R,"shape"):
585 if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
586 else: __R_shape = self.__R.shape()
587 else: raise TypeError("R has no attribute of shape: problem !")
589 if self.__Q is None: __Q_shape = (0,0)
590 elif hasattr(self.__Q,"shape"):
591 if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
592 else: __Q_shape = self.__Q.shape()
593 else: raise TypeError("Q has no attribute of shape: problem !")
595 if len(self.__H) == 0: __H_shape = (0,0)
596 elif type(self.__H) is type({}): __H_shape = (0,0)
597 elif hasattr(self.__H["Direct"],"shape"):
598 if type(self.__H["Direct"].shape) is tuple: __H_shape = self.__H["Direct"].shape
599 else: __H_shape = self.__H["Direct"].shape()
600 else: raise TypeError("H has no attribute of shape: problem !")
602 if len(self.__M) == 0: __M_shape = (0,0)
603 elif type(self.__M) is type({}): __M_shape = (0,0)
604 elif hasattr(self.__M["Direct"],"shape"):
605 if type(self.__M["Direct"].shape) is tuple: __M_shape = self.__M["Direct"].shape
606 else: __M_shape = self.__M["Direct"].shape()
607 else: raise TypeError("M has no attribute of shape: problem !")
609 # Vérification des conditions
610 # ---------------------------
611 if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
612 raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
613 if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
614 raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
616 if not( min(__B_shape) == max(__B_shape) ):
617 raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
618 if not( min(__R_shape) == max(__R_shape) ):
619 raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
620 if not( min(__Q_shape) == max(__Q_shape) ):
621 raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
622 if not( min(__M_shape) == max(__M_shape) ):
623 raise ValueError("Shape characteristic of M is incorrect: \"%s\""%(__M_shape,))
625 if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[1] == max(__Xb_shape) ):
626 raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__H_shape,__Xb_shape))
627 if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[0] == max(__Y_shape) ):
628 raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__H_shape,__Y_shape))
629 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] ):
630 raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__H_shape,__B_shape))
631 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] ):
632 raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__H_shape,__R_shape))
634 if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
635 if self.__StoredInputs["AlgorithmName"] in ["EnsembleBlue",]:
636 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
637 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
638 for member in asPersistentVector:
639 self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
640 __Xb_shape = min(__B_shape)
642 raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
644 if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
645 raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
647 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) ):
648 raise ValueError("Shape characteristic of M \"%s\" and X \"%s\" are incompatible"%(__M_shape,__Xb_shape))
652 # -----------------------------------------------------------
655 Permet de lancer le calcul d'assimilation.
657 Le nom de la méthode à activer est toujours "run". Les paramètres en
658 arguments de la méthode sont fixés. En sortie, on obtient les résultats
659 dans la variable de type dictionnaire "StoredVariables", qui contient en
660 particulier des objets de Persistence pour les analyses, OMA...
662 self.shape_validate()
664 self.__algorithm.run(
672 Parameters = self.__Parameters,
676 # -----------------------------------------------------------
677 def get(self, key=None):
679 Renvoie les résultats disponibles après l'exécution de la méthode
680 d'assimilation, ou les diagnostics disponibles. Attention, quand un
681 diagnostic porte le même nom qu'une variable stockée, c'est la variable
682 stockée qui est renvoyée, et le diagnostic est inatteignable.
685 if self.__algorithm.has_key(key):
686 return self.__algorithm.get( key )
687 elif self.__StoredInputs.has_key(key):
688 return self.__StoredInputs[key]
689 elif self.__StoredDiagnostics.has_key(key):
690 return self.__StoredDiagnostics[key]
692 raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
694 allvariables = self.__algorithm.get()
695 allvariables.update( self.__StoredDiagnostics )
696 allvariables.update( self.__StoredInputs )
699 def get_available_variables(self):
701 Renvoie les variables potentiellement utilisables pour l'étude,
702 initialement stockées comme données d'entrées ou dans les algorithmes,
703 identifiés par les chaînes de caractères. L'algorithme doit avoir été
704 préalablement choisi sinon la méthode renvoie "None".
706 if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
710 if len( self.__algorithm.keys()) > 0:
711 variables.extend( self.__algorithm.get().keys() )
712 if len( self.__StoredInputs.keys() ) > 0:
713 variables.extend( self.__StoredInputs.keys() )
717 def get_available_algorithms(self):
719 Renvoie la liste des algorithmes potentiellement utilisables, identifiés
720 par les chaînes de caractères.
723 for directory in sys.path:
724 if os.path.isdir(os.path.join(directory,"daAlgorithms")):
725 for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
726 root, ext = os.path.splitext(fname)
727 if ext == '.py' and root != '__init__':
732 def get_available_diagnostics(self):
734 Renvoie la liste des diagnostics potentiellement utilisables, identifiés
735 par les chaînes de caractères.
738 for directory in sys.path:
739 if os.path.isdir(os.path.join(directory,"daDiagnostics")):
740 for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
741 root, ext = os.path.splitext(fname)
742 if ext == '.py' and root != '__init__':
747 # -----------------------------------------------------------
748 def get_algorithms_main_path(self):
750 Renvoie le chemin pour le répertoire principal contenant les algorithmes
751 dans un sous-répertoire "daAlgorithms"
755 def add_algorithms_path(self, asPath=None):
757 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
758 se trouve un sous-répertoire "daAlgorithms"
760 Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
761 pas indispensable de le rajouter ici.
763 if not os.path.isdir(asPath):
764 raise ValueError("The given "+asPath+" argument must exist as a directory")
765 if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
766 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
767 if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
768 raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
769 sys.path.insert(0, os.path.abspath(asPath))
770 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
773 def get_diagnostics_main_path(self):
775 Renvoie le chemin pour le répertoire principal contenant les diagnostics
776 dans un sous-répertoire "daDiagnostics"
780 def add_diagnostics_path(self, asPath=None):
782 Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
783 se trouve un sous-répertoire "daDiagnostics"
785 Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
786 pas indispensable de le rajouter ici.
788 if not os.path.isdir(asPath):
789 raise ValueError("The given "+asPath+" argument must exist as a directory")
790 if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
791 raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
792 if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
793 raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
794 sys.path.insert(0, os.path.abspath(asPath))
795 sys.path = uniq( sys.path ) # Conserve en unique exemplaire chaque chemin
798 # -----------------------------------------------------------
799 def setDataObserver(self,
802 HookParameters = None,
806 Permet d'associer un observer à une ou des variables nommées gérées en
807 interne, activable selon des règles définies dans le Scheduler. A chaque
808 pas demandé dans le Scheduler, il effectue la fonction HookFunction avec
809 les arguments (variable persistante VariableName, paramètres HookParameters).
812 if type( self.__algorithm ) is dict:
813 raise ValueError("No observer can be build before choosing an algorithm.")
815 # Vérification du nom de variable et typage
816 # -----------------------------------------
817 if type( VariableName ) is str:
818 VariableNames = [VariableName,]
819 elif type( VariableName ) is list:
820 VariableNames = map( str, VariableName )
822 raise ValueError("The observer requires a name or a list of names of variables.")
824 # Association interne de l'observer à la variable
825 # -----------------------------------------------
826 for n in VariableNames:
827 if not self.__algorithm.has_key( n ):
828 raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
830 self.__algorithm.StoredVariables[ n ].setDataObserver(
831 Scheduler = Scheduler,
832 HookFunction = HookFunction,
833 HookParameters = HookParameters,
836 def removeDataObserver(self,
841 Permet de retirer un observer à une ou des variable nommée.
844 if type( self.__algorithm ) is dict:
845 raise ValueError("No observer can be removed before choosing an algorithm.")
847 # Vérification du nom de variable et typage
848 # -----------------------------------------
849 if type( VariableName ) is str:
850 VariableNames = [VariableName,]
851 elif type( VariableName ) is list:
852 VariableNames = map( str, VariableName )
854 raise ValueError("The observer requires a name or a list of names of variables.")
856 # Association interne de l'observer à la variable
857 # -----------------------------------------------
858 for n in VariableNames:
859 if not self.__algorithm.has_key( n ):
860 raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%n)
862 self.__algorithm.StoredVariables[ n ].removeDataObserver(
863 HookFunction = HookFunction,
866 # -----------------------------------------------------------
867 def setDebug(self, level=10):
869 Utiliser par exemple "import logging ; level = logging.DEBUG" avant cet
870 appel pour changer le niveau de verbosité, avec :
871 NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
874 log = logging.getLogger()
875 log.setLevel( level )
877 def unsetDebug(self):
879 Remet le logger au niveau par défaut
882 log = logging.getLogger()
883 log.setLevel( logging.WARNING )
885 def prepare_to_pickle(self):
886 self.__algorithmFile = None
887 self.__diagnosticFile = None
891 # ==============================================================================
892 if __name__ == "__main__":
893 print '\n AUTODIAGNOSTIC \n'
895 ADD = AssimilationStudy("Ma premiere etude BLUE")
897 ADD.setBackground (asVector = [0, 1, 2])
898 ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
899 ADD.setObservation (asVector = [0.5, 1.5, 2.5])
900 ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
901 ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
903 ADD.setAlgorithm(choice="Blue")
907 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
908 print "Ebauche :", [0, 1, 2]
909 print "Observation :", [0.5, 1.5, 2.5]
910 print "Demi-somme :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
911 print " qui doit être identique à :"
912 print "Analyse résultante :", ADD.get("Analysis").valueserie(0)
915 print "Algorithmes disponibles.......................:", ADD.get_available_algorithms()
916 # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path()
917 print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics()
918 # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path()
919 print "Variables disponibles.........................:", ADD.get_available_variables()
922 print "Paramètres requis par algorithme :"
923 for algo in ADD.get_available_algorithms():
924 tmpADD = AssimilationStudy("Un algorithme")
925 tmpADD.setAlgorithm(choice=algo)
926 print " %25s : %s"%(algo,tmpADD.getAlgorithmParameters())
930 ADD.setDiagnostic("RMS", "Ma RMS")
932 liste = ADD.get().keys()
934 print "Variables et diagnostics nommés disponibles...:", liste
937 print "Exemple de mise en place d'un observeur :"
938 def obs(var=None,info=None):
939 print " ---> Mise en oeuvre de l'observer"
940 print " var =",var.valueserie(-1)
942 ADD.setDataObserver( 'Analysis', HookFunction=obs, Scheduler = [2, 4], HookParameters = "Second observer")
943 # Attention, il faut décaler le stockage de 1 pour suivre le pas interne
944 # car le pas 0 correspond à l'analyse ci-dessus.
947 print "Action sur la variable observée, étape :",i
948 ADD.get('Analysis').store( [i, i, i] )
951 print "Mise en debug et hors debug"
952 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
956 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
958 print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()