2 Copyright (C) 2008-2020 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 Exigences 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 .. _section_ref_operator_one:
95 Première forme fonctionnelle : un seul opérateur direct
96 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
98 .. index:: single: OneFunction
99 .. index:: single: ScriptWithOneFunction
100 .. index:: single: DirectOperator
101 .. index:: single: DifferentialIncrement
102 .. index:: single: CenteredFiniteDifference
104 La première consiste à ne fournir qu'une seule fonction, potentiellement non
105 linéaire, et à approximer les opérateurs tangent et adjoint associés.
107 Ceci est fait dans ADAO en utilisant dans l'interface graphique EFICAS le
108 mot-clé "*ScriptWithOneFunction*" pour la description par un script. Dans
109 l'interface textuelle, c'est le mot-clé "*OneFunction*", éventuellement combiné
110 avec le mot-clé "*Script*" selon que c'est une fonction ou un script. Si c'est
111 par script externe, l'utilisateur doit fournir un fichier contenant une
112 fonction qui porte le nom obligatoire "*DirectOperator*". Par exemple, un
113 script externe peut suivre le modèle générique suivant::
115 def DirectOperator( X ):
116 """ Opérateur direct de simulation non-linéaire """
122 Dans ce cas, l'utilisateur doit aussi fournir une valeur pour l'incrément
123 différentiel (ou conserver la valeur par défaut), en utilisant dans l'interface
124 graphique (GUI) ou textuelle (TUI) le mot-clé "*DifferentialIncrement*" comme
125 paramètre, qui a une valeur par défaut de 1%. Ce coefficient est utilisé dans
126 l'approximation différences finies pour construire les opérateurs tangent et
127 adjoint. L'ordre de l'approximation différences finies peut aussi être choisi à
128 travers l'interface, en utilisant le mot-clé "*CenteredFiniteDifference*", avec
129 0 pour un schéma non centré du premier ordre (qui est la valeur par défaut), et
130 avec 1 pour un schéma centré du second ordre (et qui coûte numériquement deux
131 fois plus cher que le premier ordre). Si nécessaire et si possible, on peut
132 :ref:`subsection_ref_parallel_df`. Dans tous les cas, un mécanisme de cache
133 interne permet de limiter le nombre d'évaluations de l'opérateur au minimum
134 possible du point de vue de l'exécution séquentielle ou parallèle des
135 approximations numériques des opérateurs tangent et adjoint, pour éviter des
138 Cette première forme de définition de l'opérateur permet aisément de tester la
139 forme fonctionnelle avant son usage dans un cas ADAO, réduisant notablement la
140 complexité de l'implémentation de l'opérateur. On peut ainsi utiliser
141 l'algorithme ADAO de vérification "*FunctionTest*" (voir la section sur
142 l':ref:`section_ref_algorithm_FunctionTest`) pour ce test.
144 **Avertissement important :** le nom "*DirectOperator*" est obligatoire, et le
145 type de l'argument ``X`` peut être une liste de valeur réelles, un vecteur
146 Numpy ou une matrice Numpy. La fonction utilisateur doit accepter et traiter
149 .. _section_ref_operator_funcs:
151 Seconde forme fonctionnelle : trois opérateurs direct, tangent et adjoint
152 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
154 .. index:: single: ThreeFunctions
155 .. index:: single: ScriptWithFunctions
156 .. index:: single: DirectOperator
157 .. index:: single: TangentOperator
158 .. index:: single: AdjointOperator
162 en général, il est recommandé d'utiliser la première forme fonctionnelle
163 plutôt que la seconde. Un petit accroissement de performances n'est pas une
164 bonne raison pour utiliser l'implémentation détaillée de cette seconde forme
167 La seconde consiste à fournir directement les trois opérateurs liés :math:`O`,
168 :math:`\mathbf{O}` et :math:`\mathbf{O}^*`. C'est effectué en utilisant le
169 mot-clé "*ScriptWithFunctions*" pour la description de l'opérateur choisi dans
170 l'interface graphique EFICAS d'ADAO. Dans l'interface textuelle, c'est le
171 mot-clé "*ThreeFunctions*", éventuellement combiné avec le mot-clé "*Script*"
172 selon que c'est une fonction ou un script. L'utilisateur doit fournir dans un
173 script trois fonctions, avec les trois noms obligatoires "*DirectOperator*",
174 "*TangentOperator*" et "*AdjointOperator*". Par exemple, le script externe peut
175 suivre le squelette suivant::
177 def DirectOperator( X ):
178 """ Opérateur direct de simulation non-linéaire """
182 return quelque chose comme Y
184 def TangentOperator( paire = (X, dX) ):
185 """ Opérateur linéaire tangent, autour de X, appliqué à dX """
190 return quelque chose comme Y
192 def AdjointOperator( paire = (X, Y) ):
193 """ Opérateur adjoint, autour de X, appliqué à Y """
198 return quelque chose comme X
200 Un nouvelle fois, cette seconde définition d'opérateur permet aisément de tester
201 les formes fonctionnelles avant de les utiliser dans le cas ADAO, réduisant la
202 complexité de l'implémentation de l'opérateur.
204 Pour certains algorithmes, il faut que les fonctions tangente et adjointe
205 puisse renvoyer les matrices équivalentes à l'opérateur linéaire. Dans ce cas,
206 lorsque, respectivement, les arguments ``dX`` ou ``Y`` valent ``None``, le
207 script de l'utilisateur doit renvoyer la matrice associée.
209 **Avertissement important :** les noms "*DirectOperator*", "*TangentOperator*"
210 et "*AdjointOperator*" sont obligatoires, et le type des arguments ``X``,
211 ``Y``, ``dX`` peut être une liste de valeur réelles, un vecteur Numpy ou une
212 matrice Numpy. La fonction utilisateur doit accepter et traiter tous ces cas.
214 .. _section_ref_operator_switch:
216 Troisième forme fonctionnelle : trois opérateurs avec un branchement
217 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
219 .. index:: single: ScriptWithSwitch
220 .. index:: single: DirectOperator
221 .. index:: single: TangentOperator
222 .. index:: single: AdjointOperator
226 il est recommandé de ne pas utiliser cette troisième forme fonctionnelle sans
227 une solide raison numérique ou physique. Un accroissement de performances
228 n'est pas une bonne raison pour utiliser la complexité de cette troisième
229 forme fonctionnelle. Seule une impossibilité à utiliser les première ou
230 seconde formes justifie l'usage de la troisième.
232 La troisième forme donne de plus grandes possibilités de contrôle de
233 l'exécution des trois fonctions représentant l'opérateur, permettant un usage
234 et un contrôle avancés sur chaque exécution du code de simulation. C'est
235 réalisable en utilisant le mot-clé "*ScriptWithSwitch*" pour la description de
236 l'opérateur à travers l'interface graphique EFICAS d'ADAO. Dans l'interface
237 textuelle, il suffit d'utiliser le mot-clé "*ThreeFunctions*" précédent pour
238 définir aussi ce cas, en indiquant les fonctions adéquates. L'utilisateur doit
239 fournir un script unique aiguillant, selon un contrôle, l'exécution des formes
240 directe, tangente et adjointe du code de simulation. L'utilisateur peut alors,
241 par exemple, utiliser des approximations pour les codes tangent et adjoint, ou
242 introduire une plus grande complexité du traitement des arguments des
243 fonctions. Mais cette démarche sera plus difficile à implémenter et à déboguer.
245 Toutefois, si vous souhaitez utiliser cette troisième forme, on recommande de
246 se baser sur le modèle suivant pour le script d'aiguillage. Il nécessite un
247 fichier script ou un code externe nommé ici
248 "*Physical_simulation_functions.py*", contenant trois fonctions nommées
249 "*DirectOperator*", "*TangentOperator*" et "*AdjointOperator*" comme
250 précédemment. Voici le squelette d'aiguillage::
252 import Physical_simulation_functions
253 import numpy, logging, codecs, pickle
255 return pickle.loads(codecs.decode(data.encode(), "base64"))
258 for param in computation["specificParameters"]:
259 if param["name"] == "method":
260 method = loads(param["value"])
261 if method not in ["Direct", "Tangent", "Adjoint"]:
262 raise ValueError("No valid computation method is given")
263 logging.info("Found method is \'%s\'"%method)
265 logging.info("Loading operator functions")
266 Function = Physical_simulation_functions.DirectOperator
267 Tangent = Physical_simulation_functions.TangentOperator
268 Adjoint = Physical_simulation_functions.AdjointOperator
270 logging.info("Executing the possible computations")
272 if method == "Direct":
273 logging.info("Direct computation")
274 Xcurrent = computation["inputValues"][0][0][0]
275 data = Function(numpy.matrix( Xcurrent ).T)
276 if method == "Tangent":
277 logging.info("Tangent computation")
278 Xcurrent = computation["inputValues"][0][0][0]
279 dXcurrent = computation["inputValues"][0][0][1]
280 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
281 if method == "Adjoint":
282 logging.info("Adjoint computation")
283 Xcurrent = computation["inputValues"][0][0][0]
284 Ycurrent = computation["inputValues"][0][0][1]
285 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
287 logging.info("Formatting the output")
288 it = numpy.ravel(data)
289 outputValues = [[[[]]]]
291 outputValues[0][0][0].append(val)
294 result["outputValues"] = outputValues
295 result["specificOutputInfos"] = []
296 result["returnCode"] = 0
297 result["errorMessage"] = ""
299 Toutes les modifications envisageables peuvent être faites à partir de cette
300 hypothèse de squelette.
302 .. _section_ref_operator_control:
304 Cas spécial d'un opérateur d'évolution avec contrôle
305 ++++++++++++++++++++++++++++++++++++++++++++++++++++
307 Dans certains cas, l'opérateur d'évolution ou d'observation doit être contrôlé
308 par un contrôle d'entrée externe, qui est donné *a priori*. Dans ce cas, la
309 forme générique du modèle incrémental :math:`O` est légèrement modifiée comme
312 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
314 où :math:`\mathbf{u}` est le contrôle sur l'incrément d'état. En effet,
315 l'opérateur direct doit être appliqué à une paire de variables :math:`(X,U)`.
316 Schématiquement, l'opérateur :math:`O` doit être construit comme une fonction
317 applicable sur une paire:math:`\mathbf{(X, U)}` suit::
319 def DirectOperator( paire = (X, U) ):
320 """ Opérateur direct de simulation non-linéaire """
325 return quelque chose comme X(n+1) (évolution) ou Y(n+1) (observation)
327 Les opérateurs tangent et adjoint ont la même signature que précédemment, en
328 notant que les dérivées doivent être faites seulement partiellement par rapport
329 à :math:`\mathbf{x}`. Dans un tel cas de contrôle explicite, seule la deuxième
330 forme fonctionnelle (en utilisant "*ScriptWithFunctions*") et la troisième
331 forme fonctionnelle (en utilisant "*ScriptWithSwitch*") peuvent être utilisées.
333 .. _section_ref_operator_dimensionless:
335 Remarques complémentaires sur l'adimensionnement des opérateurs
336 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
338 .. index:: single: Adimensionnement
339 .. index:: single: Sans dimension
341 Il est fréquent que les grandeurs physiques, en entrée ou en sortie des
342 opérateurs, présentent des différences notables d'ordre de grandeur ou de taux
343 de variation. Une manière d'éviter des difficultés numériques est d'utiliser,
344 ou d'établir, un adimensionnement des calculs menés dans les opérateurs
345 [WikipediaND]_. Par principe, dans la mesure où la simulation de la physique
346 devrait être la plus adimensionnée possible, il est en premier lieu recommandé
347 d'utiliser les capacités existantes d'adimensionnement du code de calcul.
349 Néanmoins, dans le cas courant où l'on ne peut en disposer, il est souvent
350 utile d'environner le calcul pour l'adimensionner en entrée ou en sortie. Une
351 manière simple de faire cela en entrée consiste à transformer les paramètres
352 :math:`\mathbf{x}` en argument d'une fonction comme "*DirectOperator*". On
353 utilise le plus souvent comme référence les valeurs par défaut
354 :math:`\mathbf{x}^b` (ébauche, ou valeur nominale). Pourvu que chaque
355 composante de :math:`\mathbf{x}^b` soit non nulle, on peut ensuite procéder par
356 correction multiplicative. Pour cela, on peut par exemple poser:
358 .. math:: \mathbf{x} = \mathbf{\alpha}\mathbf{x}^b
360 et optimiser ensuite le paramètre multiplicatif :math:`\mathbf{\alpha}`. Ce
361 paramètre a pour valeur par défaut (ou pour ébauche) un vecteur de 1. De
362 manière similaire, on peut procéder par correction additive si c'est plus
363 judicieux pour la physique sous-jacente. Ainsi, dans ce cas, on peut poser:
365 .. math:: \mathbf{x} =\mathbf{x}^b + \mathbf{\alpha}
367 et optimiser ensuite le paramètre additif :math:`\mathbf{\alpha}`. Cette fois,
368 ce paramètre a pour valeur d'ébauche un vecteur de 0.
370 Attention, l'application d'une démarche d'adimensionnement nécessite aussi la
371 modification des covariances d'erreurs associées dans la formulation globale du
372 problème d'optimisation.
374 Une telle démarche suffit rarement à éviter tous les problèmes numériques, mais
375 permet souvent d'améliorer beaucoup le conditionnement numérique de
378 Gestion explicite de fonctions "multiples"
379 ++++++++++++++++++++++++++++++++++++++++++
383 il est fortement recommandé de ne pas utiliser cette gestion explicite de
384 fonctions "multiples" sans une très solide raison informatique pour le faire.
385 Cette gestion est déjà effectuée par défaut dans ADAO pour l'amélioration des
386 performances. Seul l'utilisateur très averti, cherchant à gérer des cas
387 particulièrement difficiles, peut s'intéresser à cette extension. En dépit de
388 sa simplicité, c'est au risque explicite de dégrader notablement les
391 Il est possible, lorsque l'on fournit des fonctions d'opérateurs, de les
392 définir comme des fonctions qui traitent non pas un seul argument, mais une
393 série d'arguments, pour restituer en sortie la série des valeurs
394 correspondantes. En pseudo-code, la fonction "multiple", ici nommée
395 ``MultiFunctionO``, représentant l'opérateur classique :math:`O` nommé
396 "*DirectOperator*", effectue::
398 def MultiFunctionO( Inputs ):
402 Y = DirectOperator( X )
406 La longueur de la sortie (c'est-à-dire le nombre de valeurs calculées) est
407 égale à la longueur de l'entrée (c'est-à-dire le nombre d'états dont on veut
408 calculer la valeur par l'opérateur).
410 Cette possibilité n'est disponible que dans l'interface textuelle d'ADAO. Pour
411 cela, lors de la définition d'une fonction d'opérateur, en même temps que l'on
412 définit de manière habituelle la fonction ou le script externe, il suffit
413 d'indiquer en plus en argument par un booléen "*InputFunctionAsMulti*" que la
414 définition est celle d'une fonction "multiple".