Salome HOME
Updating version and copyright date information
[modules/adao.git] / doc / fr / ref_operator_requirements.rst
1 ..
2    Copyright (C) 2008-2020 EDF R&D
3
4    This file is part of SALOME ADAO module.
5
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.
10
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.
15
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
19
20    See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21
22    Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
23
24 .. _section_ref_operator_requirements:
25
26 Exigences pour les fonctions décrivant un opérateur
27 ---------------------------------------------------
28
29 .. index:: single: setObservationOperator
30 .. index:: single: setEvolutionModel
31 .. index:: single: setControlModel
32
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.
40
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:
46
47 .. math:: \mathbf{y} = O( \mathbf{x} )
48
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.
55
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).
60
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
63 fonctionnelle**.
64
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.
72
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:
75
76 - :ref:`section_ref_operator_one`
77 - :ref:`section_ref_operator_funcs`
78 - :ref:`section_ref_operator_switch`
79
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:
85
86   .. eficas_operator_function:
87   .. image:: images/eficas_operator_function.png
88     :align: center
89     :width: 100%
90   .. centered::
91     **Choisir graphiquement une représentation fonctionnelle de l'opérateur**
92
93 .. _section_ref_operator_one:
94
95 Première forme fonctionnelle : un seul opérateur direct
96 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
97
98 .. index:: single: OneFunction
99 .. index:: single: ScriptWithOneFunction
100 .. index:: single: DirectOperator
101 .. index:: single: DifferentialIncrement
102 .. index:: single: CenteredFiniteDifference
103
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.
106
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::
114
115     def DirectOperator( X ):
116         """ Opérateur direct de simulation non-linéaire """
117         ...
118         ...
119         ...
120         return Y=O(X)
121
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
136 calculs redondants.
137
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.
143
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
147 tous ces cas.
148
149 .. _section_ref_operator_funcs:
150
151 Seconde forme fonctionnelle : trois opérateurs direct, tangent et adjoint
152 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
153
154 .. index:: single: ThreeFunctions
155 .. index:: single: ScriptWithFunctions
156 .. index:: single: DirectOperator
157 .. index:: single: TangentOperator
158 .. index:: single: AdjointOperator
159
160 .. warning::
161
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
165   fonctionnelle.
166
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::
176
177     def DirectOperator( X ):
178         """ Opérateur direct de simulation non-linéaire """
179         ...
180         ...
181         ...
182         return quelque chose comme Y
183
184     def TangentOperator( paire = (X, dX) ):
185         """ Opérateur linéaire tangent, autour de X, appliqué à dX """
186         X, dX = paire
187         ...
188         ...
189         ...
190         return quelque chose comme Y
191
192     def AdjointOperator( paire = (X, Y) ):
193         """ Opérateur adjoint, autour de X, appliqué à Y """
194         X, Y = paire
195         ...
196         ...
197         ...
198         return quelque chose comme X
199
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.
203
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.
208
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.
213
214 .. _section_ref_operator_switch:
215
216 Troisième forme fonctionnelle : trois opérateurs avec un branchement
217 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
218
219 .. index:: single: ScriptWithSwitch
220 .. index:: single: DirectOperator
221 .. index:: single: TangentOperator
222 .. index:: single: AdjointOperator
223
224 .. warning::
225
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.
231
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.
244
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::
251
252     import Physical_simulation_functions
253     import numpy, logging, codecs, pickle
254     def loads( data ):
255         return pickle.loads(codecs.decode(data.encode(), "base64"))
256     #
257     method = ""
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)
264     #
265     logging.info("Loading operator functions")
266     Function = Physical_simulation_functions.DirectOperator
267     Tangent  = Physical_simulation_functions.TangentOperator
268     Adjoint  = Physical_simulation_functions.AdjointOperator
269     #
270     logging.info("Executing the possible computations")
271     data = []
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))
286     #
287     logging.info("Formatting the output")
288     it = numpy.ravel(data)
289     outputValues = [[[[]]]]
290     for val in it:
291       outputValues[0][0][0].append(val)
292     #
293     result = {}
294     result["outputValues"]        = outputValues
295     result["specificOutputInfos"] = []
296     result["returnCode"]          = 0
297     result["errorMessage"]        = ""
298
299 Toutes les modifications envisageables peuvent être faites à partir de cette
300 hypothèse de squelette.
301
302 .. _section_ref_operator_control:
303
304 Cas spécial d'un opérateur d'évolution avec contrôle
305 ++++++++++++++++++++++++++++++++++++++++++++++++++++
306
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
310 suit:
311
312 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
313
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::
318
319     def DirectOperator( paire = (X, U) ):
320         """ Opérateur direct de simulation non-linéaire """
321         X, U = paire
322         ...
323         ...
324         ...
325         return quelque chose comme X(n+1) (évolution) ou Y(n+1) (observation)
326
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.
332
333 .. _section_ref_operator_dimensionless:
334
335 Remarques complémentaires sur l'adimensionnement des opérateurs
336 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
337
338 .. index:: single: Adimensionnement
339 .. index:: single: Sans dimension
340
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.
348
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:
357
358 .. math:: \mathbf{x} = \mathbf{\alpha}\mathbf{x}^b
359
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:
364
365 .. math:: \mathbf{x} =\mathbf{x}^b + \mathbf{\alpha}
366
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.
369
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.
373
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
376 l'optimisation.
377
378 Gestion explicite de fonctions "multiples"
379 ++++++++++++++++++++++++++++++++++++++++++
380
381 .. warning::
382
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
389   performances.
390
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::
397
398     def MultiFunctionO( Inputs ):
399         """ Multiple ! """
400         Outputs = []
401         for X in Inputs:
402             Y = DirectOperator( X )
403             Outputs.append( Y )
404         return Outputs
405
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).
409
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".