2 Copyright (C) 2008-2022 EDF R&D
4 This file is part of SALOME ADAO module.
6 This library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
11 This library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
16 You should have received a copy of the GNU Lesser General Public
17 License along with this library; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 .. _section_ref_operator_requirements:
26 Conditions requises pour les fonctions décrivant un opérateur
27 -------------------------------------------------------------
29 .. index:: single: setObservationOperator
30 .. index:: single: setEvolutionModel
31 .. index:: single: setControlModel
33 Les opérateurs d'observation et d'évolution sont nécessaires pour mettre en
34 oeuvre les procédures d'assimilation de données ou d'optimisation. Ils
35 comprennent la simulation physique par des calculs numériques, mais aussi le
36 filtrage et de restriction pour comparer la simulation à l'observation.
37 L'opérateur d'évolution est ici considéré dans sa forme incrémentale, qui
38 représente la transition entre deux états successifs, et il est alors similaire
39 à l'opérateur d'observation.
41 Schématiquement, un opérateur :math:`O` a pour objet de restituer une solution
42 pour des paramètres d'entrée spécifiés. Une partie des paramètres d'entrée peut
43 être modifiée au cours de la procédure d'optimisation. Ainsi, la représentation
44 mathématique d'un tel processus est une fonction. Il a été brièvement décrit
45 dans la section :ref:`section_theory` et il est généralisé ici par la relation:
47 .. math:: \mathbf{y} = O( \mathbf{x} )
49 entre les pseudo-observations en sortie :math:`\mathbf{y}` et les paramètres
50 d'entrée :math:`\mathbf{x}` en utilisant l'opérateur d'observation ou
51 d'évolution :math:`O`. La même représentation fonctionnelle peut être utilisée
52 pour le modèle linéaire tangent :math:`\mathbf{O}` de :math:`O` et son adjoint
53 :math:`\mathbf{O}^*` qui sont aussi requis par certains algorithmes
54 d'assimilation de données ou d'optimisation.
56 En entrée et en sortie de ces opérateurs, les variables :math:`\mathbf{x}` et
57 :math:`\mathbf{y}`, ou leurs incréments, sont mathématiquement des vecteurs, et
58 ils peuvent donc être donnés par l'utilisateur comme des vecteurs non-orientés
59 (de type liste ou vecteur Numpy) ou orientés (de type matrice Numpy).
61 Ainsi, **pour décrire de manière complète un opérateur, l'utilisateur n'a qu'à
62 fournir une fonction qui réalise complètement et uniquement l'opération
65 Cette fonction est généralement donnée comme une fonction ou un script Python,
66 qui peuvent en particulier être exécuté comme une fonction Python indépendante
67 ou dans un noeud YACS. Cette fonction ou ce script peuvent, sans différences,
68 lancer des codes externes ou utiliser des appels et des méthodes internes
69 Python ou SALOME. Si l'algorithme nécessite les 3 aspects de l'opérateur (forme
70 directe, forme tangente et forme adjointe), l'utilisateur doit donner les 3
71 fonctions ou les approximer grâce à ADAO.
73 Il existe pour l'utilisateur 3 méthodes effectives de fournir une représentation
74 fonctionnelle de l'opérateur, qui diffèrent selon le type d'argument choisi:
76 - :ref:`section_ref_operator_one`
77 - :ref:`section_ref_operator_funcs`
78 - :ref:`section_ref_operator_switch`
80 Dans le cas de l'interface textuelle d'ADAO (TUI), seules les deux premières
81 sont nécessaires car la troisième est incluse dans la seconde. Dans le cas de
82 l'interface graphique EFICAS d'ADAO, ces méthodes sont choisies dans le champ
83 "*FROM*" de chaque opérateur ayant une valeur "*Function*" comme
84 "*INPUT_TYPE*", comme le montre la figure suivante :
86 .. eficas_operator_function:
87 .. image:: images/eficas_operator_function.png
91 **Choisir graphiquement une représentation fonctionnelle de l'opérateur**
93 En interface textuelle d'ADAO (TUI), dans le cas précis illustré ci-dessus, on
94 réalise la même démarche en écrivant :
98 case.set( 'ObservationOperator',
100 Script = 'scripts_for_JDC.py'
104 .. _section_ref_operator_one:
106 Première forme fonctionnelle : un seul opérateur direct
107 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
109 .. index:: single: OneFunction
110 .. index:: single: ScriptWithOneFunction
111 .. index:: single: DirectOperator
112 .. index:: single: DifferentialIncrement
113 .. index:: single: CenteredFiniteDifference
115 La première consiste à ne fournir qu'une seule fonction, potentiellement non
116 linéaire, et à approximer les opérateurs tangent et adjoint associés.
118 Ceci est fait dans ADAO en utilisant, dans l'interface graphique EFICAS d'ADAO,
119 le mot-clé "*ScriptWithOneFunction*" pour la description par un script. Dans
120 l'interface textuelle, c'est le mot-clé "*OneFunction*", éventuellement combiné
121 avec le mot-clé "*Script*" selon que c'est une fonction ou un script. Si c'est
122 par script externe, l'utilisateur doit fournir un fichier contenant une
123 fonction qui porte le nom obligatoire "*DirectOperator*". Par exemple, un
124 script externe peut suivre le modèle générique suivant::
126 def DirectOperator( X ):
127 """ Opérateur direct de simulation non-linéaire """
133 Dans ce cas, l'utilisateur doit aussi fournir une valeur pour l'incrément
134 différentiel (ou conserver la valeur par défaut), en utilisant dans l'interface
135 graphique (GUI) ou textuelle (TUI) le mot-clé "*DifferentialIncrement*" comme
136 paramètre, qui a une valeur par défaut de 1%. Ce coefficient est utilisé dans
137 l'approximation différences finies pour construire les opérateurs tangent et
138 adjoint. L'ordre de l'approximation différences finies peut aussi être choisi à
139 travers l'interface, en utilisant le mot-clé "*CenteredFiniteDifference*", avec
140 0 pour un schéma non centré du premier ordre (qui est la valeur par défaut), et
141 avec 1 pour un schéma centré du second ordre (et qui coûte numériquement deux
142 fois plus cher que le premier ordre). Si nécessaire et si possible, on peut
143 :ref:`subsection_ref_parallel_df`. Dans tous les cas, un mécanisme de cache
144 interne permet de limiter le nombre d'évaluations de l'opérateur au minimum
145 possible du point de vue de l'exécution séquentielle ou parallèle des
146 approximations numériques des opérateurs tangent et adjoint, pour éviter des
149 Cette première forme de définition de l'opérateur permet aisément de tester la
150 forme fonctionnelle avant son usage dans un cas ADAO, réduisant notablement la
151 complexité de l'implémentation de l'opérateur. On peut ainsi utiliser
152 l'algorithme ADAO de vérification "*FunctionTest*" (voir la section sur
153 l':ref:`section_ref_algorithm_FunctionTest`) pour ce test.
155 **Avertissement important :** le nom "*DirectOperator*" est obligatoire, et le
156 type de l'argument ``X`` peut être une liste de valeur réelles, un vecteur
157 Numpy ou une matrice Numpy. La fonction utilisateur doit accepter et traiter
160 .. _section_ref_operator_funcs:
162 Seconde forme fonctionnelle : trois opérateurs direct, tangent et adjoint
163 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
165 .. index:: single: ThreeFunctions
166 .. index:: single: ScriptWithFunctions
167 .. index:: single: DirectOperator
168 .. index:: single: TangentOperator
169 .. index:: single: AdjointOperator
173 En général, il est recommandé d'utiliser la première forme fonctionnelle
174 plutôt que la seconde. Un petit accroissement de performances n'est pas une
175 bonne raison pour utiliser l'implémentation détaillée de cette seconde forme
178 La seconde consiste à fournir directement les trois opérateurs liés :math:`O`,
179 :math:`\mathbf{O}` et :math:`\mathbf{O}^*`. C'est effectué en utilisant le
180 mot-clé "*ScriptWithFunctions*" pour la description de l'opérateur choisi dans
181 l'interface graphique EFICAS d'ADAO. Dans l'interface textuelle, c'est le
182 mot-clé "*ThreeFunctions*", éventuellement combiné avec le mot-clé "*Script*"
183 selon que c'est une fonction ou un script. L'utilisateur doit fournir dans un
184 script trois fonctions, avec les trois noms obligatoires "*DirectOperator*",
185 "*TangentOperator*" et "*AdjointOperator*". Par exemple, le script externe peut
186 suivre le squelette suivant::
188 def DirectOperator( X ):
189 """ Opérateur direct de simulation non-linéaire """
193 return "un vecteur similaire à Y"
195 def TangentOperator( paire = (X, dX) ):
196 """ Opérateur linéaire tangent, autour de X, appliqué à dX """
201 return "un vecteur similaire à Y"
203 def AdjointOperator( paire = (X, Y) ):
204 """ Opérateur adjoint, autour de X, appliqué à Y """
209 return "un vecteur similaire à X"
211 Un nouvelle fois, cette seconde définition d'opérateur permet aisément de tester
212 les formes fonctionnelles avant de les utiliser dans le cas ADAO, réduisant la
213 complexité de l'implémentation de l'opérateur.
215 Pour certains algorithmes (en particulier les filtres non ensemblistes), il
216 faut que les fonctions tangente et adjointe puisse renvoyer les matrices
217 équivalentes à l'opérateur linéaire. Dans ce cas, lorsque, respectivement, les
218 arguments ``dX`` ou ``Y`` valent ``None``, le script de l'utilisateur doit
219 renvoyer la matrice associée. Les squelettes des fonctions "*TangentOperator*"
220 et "*AdjointOperator*" deviennent alors les suivants::
222 def TangentOperator( paire = (X, dX) ):
223 """ Opérateur linéaire tangent, autour de X, appliqué à dX """
228 if dX is None or len(dX) == 0:
229 return "la matrice de l'opérateur linéaire tangent"
231 return "un vecteur similaire à Y"
233 def AdjointOperator( paire = (X, Y) ):
234 """ Opérateur adjoint, autour de X, appliqué à Y """
239 if Y is None or len(Y) == 0:
240 return "la matrice de l'opérateur linéaire adjoint"
242 return "un vecteur similaire à X"
244 **Avertissement important :** les noms "*DirectOperator*", "*TangentOperator*"
245 et "*AdjointOperator*" sont obligatoires, et le type des arguments ``X``,
246 ``Y``, ``dX`` peut être une liste de valeur réelles, un vecteur Numpy ou une
247 matrice Numpy. La fonction utilisateur doit accepter et traiter tous ces cas.
249 .. _section_ref_operator_switch:
251 Troisième forme fonctionnelle : trois opérateurs avec un branchement
252 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
254 .. index:: single: ScriptWithSwitch
255 .. index:: single: DirectOperator
256 .. index:: single: TangentOperator
257 .. index:: single: AdjointOperator
261 Il est recommandé de ne pas utiliser cette troisième forme fonctionnelle sans
262 une solide raison numérique ou physique. Un accroissement de performances
263 n'est pas une bonne raison pour utiliser la complexité de cette troisième
264 forme fonctionnelle. Seule une impossibilité à utiliser les première ou
265 seconde formes justifie l'usage de la troisième.
267 La troisième forme donne de plus grandes possibilités de contrôle de
268 l'exécution des trois fonctions représentant l'opérateur, permettant un usage
269 et un contrôle avancés sur chaque exécution du code de simulation. C'est
270 réalisable en utilisant le mot-clé "*ScriptWithSwitch*" pour la description de
271 l'opérateur à travers l'interface graphique EFICAS d'ADAO. Dans l'interface
272 textuelle, il suffit d'utiliser le mot-clé "*ThreeFunctions*" précédent pour
273 définir aussi ce cas, en indiquant les fonctions adéquates. L'utilisateur doit
274 fournir un script unique aiguillant, selon un contrôle, l'exécution des formes
275 directe, tangente et adjointe du code de simulation. L'utilisateur peut alors,
276 par exemple, utiliser des approximations pour les codes tangent et adjoint, ou
277 introduire une plus grande complexité du traitement des arguments des
278 fonctions. Mais cette démarche sera plus difficile à implémenter et à déboguer.
280 Toutefois, si vous souhaitez utiliser cette troisième forme, on recommande de
281 se baser sur le modèle suivant pour le script d'aiguillage. Il nécessite un
282 fichier script ou un code externe nommé ici
283 "*Physical_simulation_functions.py*", contenant trois fonctions nommées
284 "*DirectOperator*", "*TangentOperator*" et "*AdjointOperator*" comme
285 précédemment. Voici le squelette d'aiguillage:
288 import Physical_simulation_functions
289 import numpy, logging, codecs, pickle
291 return pickle.loads(codecs.decode(data.encode(), "base64"))
294 for param in computation["specificParameters"]:
295 if param["name"] == "method":
296 method = loads(param["value"])
297 if method not in ["Direct", "Tangent", "Adjoint"]:
298 raise ValueError("No valid computation method is given")
299 logging.info("Found method is \'%s\'"%method)
301 logging.info("Loading operator functions")
302 Function = Physical_simulation_functions.DirectOperator
303 Tangent = Physical_simulation_functions.TangentOperator
304 Adjoint = Physical_simulation_functions.AdjointOperator
306 logging.info("Executing the possible computations")
308 if method == "Direct":
309 logging.info("Direct computation")
310 Xcurrent = computation["inputValues"][0][0][0]
311 data = Function(numpy.matrix( Xcurrent ).T)
312 if method == "Tangent":
313 logging.info("Tangent computation")
314 Xcurrent = computation["inputValues"][0][0][0]
315 dXcurrent = computation["inputValues"][0][0][1]
316 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
317 if method == "Adjoint":
318 logging.info("Adjoint computation")
319 Xcurrent = computation["inputValues"][0][0][0]
320 Ycurrent = computation["inputValues"][0][0][1]
321 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
323 logging.info("Formatting the output")
324 it = numpy.ravel(data)
325 outputValues = [[[[]]]]
327 outputValues[0][0][0].append(val)
330 result["outputValues"] = outputValues
331 result["specificOutputInfos"] = []
332 result["returnCode"] = 0
333 result["errorMessage"] = ""
335 Toutes les modifications envisageables peuvent être faites à partir de cette
336 hypothèse de squelette.
338 .. _section_ref_operator_control:
340 Cas spécial d'un opérateur d'évolution avec contrôle
341 ++++++++++++++++++++++++++++++++++++++++++++++++++++
343 Dans certains cas, l'opérateur d'évolution ou d'observation doit être contrôlé
344 par un contrôle d'entrée externe, qui est donné *a priori*. Dans ce cas, la
345 forme générique du modèle incrémental :math:`O` est légèrement modifiée comme
348 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
350 où :math:`\mathbf{u}` est le contrôle sur l'incrément d'état. En effet,
351 l'opérateur direct doit être appliqué à une paire de variables :math:`(X,U)`.
352 Schématiquement, l'opérateur :math:`O` doit être construit comme une fonction
353 applicable sur une paire :math:`\mathbf{(X, U)}` comme suit :
356 def DirectOperator( paire = (X, U) ):
357 """ Opérateur direct de simulation non-linéaire """
362 return quelque chose comme X(n+1) (évolution) ou Y(n+1) (observation)
364 Les opérateurs tangent et adjoint ont la même signature que précédemment, en
365 notant que les dérivées doivent être faites seulement partiellement par rapport
366 à :math:`\mathbf{x}`. Dans un tel cas de contrôle explicite, seule la deuxième
367 forme fonctionnelle (en utilisant "*ScriptWithFunctions*") et la troisième
368 forme fonctionnelle (en utilisant "*ScriptWithSwitch*") peuvent être utilisées.
370 .. _section_ref_operator_dimensionless:
372 Remarques complémentaires sur l'adimensionnement des opérateurs
373 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
375 .. index:: single: Adimensionnement
376 .. index:: single: Sans dimension
378 Il est fréquent que les grandeurs physiques, en entrée ou en sortie des
379 opérateurs, présentent des différences notables d'ordre de grandeur ou de taux
380 de variation. Une manière d'éviter des difficultés numériques est d'utiliser,
381 ou d'établir, un adimensionnement des calculs menés dans les opérateurs
382 [WikipediaND]_. Par principe, dans la mesure où la simulation de la physique
383 devrait être la plus adimensionnée possible, il est en premier lieu recommandé
384 d'utiliser les capacités existantes d'adimensionnement du code de calcul.
386 Néanmoins, dans le cas courant où l'on ne peut en disposer, il est souvent
387 utile d'environner le calcul pour l'adimensionner en entrée ou en sortie. Une
388 manière simple de faire cela en entrée consiste à transformer les paramètres
389 :math:`\mathbf{x}` en argument d'une fonction comme "*DirectOperator*". On
390 utilise le plus souvent comme référence les valeurs par défaut
391 :math:`\mathbf{x}^b` (ébauche, ou valeur nominale). Pourvu que chaque
392 composante de :math:`\mathbf{x}^b` soit non nulle, on peut ensuite procéder par
393 correction multiplicative. Pour cela, on peut par exemple poser :
395 .. math:: \mathbf{x} = \mathbf{\alpha}\mathbf{x}^b
397 et optimiser ensuite le paramètre multiplicatif :math:`\mathbf{\alpha}`. Ce
398 paramètre a pour valeur par défaut (ou pour ébauche) un vecteur de 1. De
399 manière similaire, on peut procéder par correction additive si c'est plus
400 judicieux pour la physique sous-jacente. Ainsi, dans ce cas, on peut poser :
402 .. math:: \mathbf{x} =\mathbf{x}^b + \mathbf{\alpha}
404 et optimiser ensuite le paramètre additif :math:`\mathbf{\alpha}`. Cette fois,
405 ce paramètre a pour valeur d'ébauche un vecteur de 0.
407 Attention, l'application d'une démarche d'adimensionnement nécessite aussi la
408 modification des covariances d'erreurs associées dans la formulation globale du
409 problème d'optimisation.
411 Une telle démarche suffit rarement à éviter tous les problèmes numériques, mais
412 permet souvent d'améliorer beaucoup le conditionnement numérique de
415 Gestion explicite de fonctions "multiples"
416 ++++++++++++++++++++++++++++++++++++++++++
420 Il est fortement recommandé de ne pas utiliser cette gestion explicite de
421 fonctions "multiples" sans une très solide raison informatique pour le faire.
422 Cette gestion est déjà effectuée par défaut dans ADAO pour l'amélioration des
423 performances. Seul l'utilisateur très averti, cherchant à gérer des cas
424 particulièrement difficiles, peut s'intéresser à cette extension. En dépit de
425 sa simplicité, c'est au risque explicite de dégrader notablement les
428 Il est possible, lorsque l'on fournit des fonctions d'opérateurs, de les
429 définir comme des fonctions qui traitent non pas un seul argument, mais une
430 série d'arguments, pour restituer en sortie la série des valeurs
431 correspondantes. En pseudo-code, la fonction "multiple", ici nommée
432 ``MultiFunctionO``, représentant l'opérateur classique :math:`O` nommé
433 "*DirectOperator*", effectue :
436 def MultiFunctionO( Inputs ):
440 Y = DirectOperator( X )
444 La longueur de la sortie (c'est-à-dire le nombre de valeurs calculées) est
445 égale à la longueur de l'entrée (c'est-à-dire le nombre d'états dont on veut
446 calculer la valeur par l'opérateur).
448 Cette possibilité n'est disponible que dans l'interface textuelle d'ADAO. Pour
449 cela, lors de la définition d'une fonction d'opérateur, en même temps que l'on
450 définit de manière habituelle la fonction ou le script externe, il suffit
451 d'indiquer en plus en argument par un booléen "*InputFunctionAsMulti*" que la
452 définition est celle d'une fonction "multiple".