2 Copyright (C) 2008-2021 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 d'ADAO,
108 le 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 "un vecteur similaire à Y"
184 def TangentOperator( paire = (X, dX) ):
185 """ Opérateur linéaire tangent, autour de X, appliqué à dX """
190 return "un vecteur similaire à Y"
192 def AdjointOperator( paire = (X, Y) ):
193 """ Opérateur adjoint, autour de X, appliqué à Y """
198 return "un vecteur similaire à 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 (en particulier les filtres non ensemblistes), il
205 faut que les fonctions tangente et adjointe puisse renvoyer les matrices
206 équivalentes à l'opérateur linéaire. Dans ce cas, lorsque, respectivement, les
207 arguments ``dX`` ou ``Y`` valent ``None``, le script de l'utilisateur doit
208 renvoyer la matrice associée. Les squelettes des fonctions "*TangentOperator*"
209 et "*AdjointOperator*" deviennent alors les suivants::
211 def TangentOperator( paire = (X, dX) ):
212 """ Opérateur linéaire tangent, autour de X, appliqué à dX """
217 if dX is None or len(dX) == 0:
218 return "la matrice de l'opérateur linéaire tangent"
220 return "un vecteur similaire à Y"
222 def AdjointOperator( paire = (X, Y) ):
223 """ Opérateur adjoint, autour de X, appliqué à Y """
228 if Y is None or len(Y) == 0:
229 return "la matrice de l'opérateur linéaire adjoint"
231 return "un vecteur similaire à X"
233 **Avertissement important :** les noms "*DirectOperator*", "*TangentOperator*"
234 et "*AdjointOperator*" sont obligatoires, et le type des arguments ``X``,
235 ``Y``, ``dX`` peut être une liste de valeur réelles, un vecteur Numpy ou une
236 matrice Numpy. La fonction utilisateur doit accepter et traiter tous ces cas.
238 .. _section_ref_operator_switch:
240 Troisième forme fonctionnelle : trois opérateurs avec un branchement
241 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
243 .. index:: single: ScriptWithSwitch
244 .. index:: single: DirectOperator
245 .. index:: single: TangentOperator
246 .. index:: single: AdjointOperator
250 il est recommandé de ne pas utiliser cette troisième forme fonctionnelle sans
251 une solide raison numérique ou physique. Un accroissement de performances
252 n'est pas une bonne raison pour utiliser la complexité de cette troisième
253 forme fonctionnelle. Seule une impossibilité à utiliser les première ou
254 seconde formes justifie l'usage de la troisième.
256 La troisième forme donne de plus grandes possibilités de contrôle de
257 l'exécution des trois fonctions représentant l'opérateur, permettant un usage
258 et un contrôle avancés sur chaque exécution du code de simulation. C'est
259 réalisable en utilisant le mot-clé "*ScriptWithSwitch*" pour la description de
260 l'opérateur à travers l'interface graphique EFICAS d'ADAO. Dans l'interface
261 textuelle, il suffit d'utiliser le mot-clé "*ThreeFunctions*" précédent pour
262 définir aussi ce cas, en indiquant les fonctions adéquates. L'utilisateur doit
263 fournir un script unique aiguillant, selon un contrôle, l'exécution des formes
264 directe, tangente et adjointe du code de simulation. L'utilisateur peut alors,
265 par exemple, utiliser des approximations pour les codes tangent et adjoint, ou
266 introduire une plus grande complexité du traitement des arguments des
267 fonctions. Mais cette démarche sera plus difficile à implémenter et à déboguer.
269 Toutefois, si vous souhaitez utiliser cette troisième forme, on recommande de
270 se baser sur le modèle suivant pour le script d'aiguillage. Il nécessite un
271 fichier script ou un code externe nommé ici
272 "*Physical_simulation_functions.py*", contenant trois fonctions nommées
273 "*DirectOperator*", "*TangentOperator*" et "*AdjointOperator*" comme
274 précédemment. Voici le squelette d'aiguillage::
276 import Physical_simulation_functions
277 import numpy, logging, codecs, pickle
279 return pickle.loads(codecs.decode(data.encode(), "base64"))
282 for param in computation["specificParameters"]:
283 if param["name"] == "method":
284 method = loads(param["value"])
285 if method not in ["Direct", "Tangent", "Adjoint"]:
286 raise ValueError("No valid computation method is given")
287 logging.info("Found method is \'%s\'"%method)
289 logging.info("Loading operator functions")
290 Function = Physical_simulation_functions.DirectOperator
291 Tangent = Physical_simulation_functions.TangentOperator
292 Adjoint = Physical_simulation_functions.AdjointOperator
294 logging.info("Executing the possible computations")
296 if method == "Direct":
297 logging.info("Direct computation")
298 Xcurrent = computation["inputValues"][0][0][0]
299 data = Function(numpy.matrix( Xcurrent ).T)
300 if method == "Tangent":
301 logging.info("Tangent computation")
302 Xcurrent = computation["inputValues"][0][0][0]
303 dXcurrent = computation["inputValues"][0][0][1]
304 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
305 if method == "Adjoint":
306 logging.info("Adjoint computation")
307 Xcurrent = computation["inputValues"][0][0][0]
308 Ycurrent = computation["inputValues"][0][0][1]
309 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
311 logging.info("Formatting the output")
312 it = numpy.ravel(data)
313 outputValues = [[[[]]]]
315 outputValues[0][0][0].append(val)
318 result["outputValues"] = outputValues
319 result["specificOutputInfos"] = []
320 result["returnCode"] = 0
321 result["errorMessage"] = ""
323 Toutes les modifications envisageables peuvent être faites à partir de cette
324 hypothèse de squelette.
326 .. _section_ref_operator_control:
328 Cas spécial d'un opérateur d'évolution avec contrôle
329 ++++++++++++++++++++++++++++++++++++++++++++++++++++
331 Dans certains cas, l'opérateur d'évolution ou d'observation doit être contrôlé
332 par un contrôle d'entrée externe, qui est donné *a priori*. Dans ce cas, la
333 forme générique du modèle incrémental :math:`O` est légèrement modifiée comme
336 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
338 où :math:`\mathbf{u}` est le contrôle sur l'incrément d'état. En effet,
339 l'opérateur direct doit être appliqué à une paire de variables :math:`(X,U)`.
340 Schématiquement, l'opérateur :math:`O` doit être construit comme une fonction
341 applicable sur une paire:math:`\mathbf{(X, U)}` suit::
343 def DirectOperator( paire = (X, U) ):
344 """ Opérateur direct de simulation non-linéaire """
349 return quelque chose comme X(n+1) (évolution) ou Y(n+1) (observation)
351 Les opérateurs tangent et adjoint ont la même signature que précédemment, en
352 notant que les dérivées doivent être faites seulement partiellement par rapport
353 à :math:`\mathbf{x}`. Dans un tel cas de contrôle explicite, seule la deuxième
354 forme fonctionnelle (en utilisant "*ScriptWithFunctions*") et la troisième
355 forme fonctionnelle (en utilisant "*ScriptWithSwitch*") peuvent être utilisées.
357 .. _section_ref_operator_dimensionless:
359 Remarques complémentaires sur l'adimensionnement des opérateurs
360 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
362 .. index:: single: Adimensionnement
363 .. index:: single: Sans dimension
365 Il est fréquent que les grandeurs physiques, en entrée ou en sortie des
366 opérateurs, présentent des différences notables d'ordre de grandeur ou de taux
367 de variation. Une manière d'éviter des difficultés numériques est d'utiliser,
368 ou d'établir, un adimensionnement des calculs menés dans les opérateurs
369 [WikipediaND]_. Par principe, dans la mesure où la simulation de la physique
370 devrait être la plus adimensionnée possible, il est en premier lieu recommandé
371 d'utiliser les capacités existantes d'adimensionnement du code de calcul.
373 Néanmoins, dans le cas courant où l'on ne peut en disposer, il est souvent
374 utile d'environner le calcul pour l'adimensionner en entrée ou en sortie. Une
375 manière simple de faire cela en entrée consiste à transformer les paramètres
376 :math:`\mathbf{x}` en argument d'une fonction comme "*DirectOperator*". On
377 utilise le plus souvent comme référence les valeurs par défaut
378 :math:`\mathbf{x}^b` (ébauche, ou valeur nominale). Pourvu que chaque
379 composante de :math:`\mathbf{x}^b` soit non nulle, on peut ensuite procéder par
380 correction multiplicative. Pour cela, on peut par exemple poser:
382 .. math:: \mathbf{x} = \mathbf{\alpha}\mathbf{x}^b
384 et optimiser ensuite le paramètre multiplicatif :math:`\mathbf{\alpha}`. Ce
385 paramètre a pour valeur par défaut (ou pour ébauche) un vecteur de 1. De
386 manière similaire, on peut procéder par correction additive si c'est plus
387 judicieux pour la physique sous-jacente. Ainsi, dans ce cas, on peut poser:
389 .. math:: \mathbf{x} =\mathbf{x}^b + \mathbf{\alpha}
391 et optimiser ensuite le paramètre additif :math:`\mathbf{\alpha}`. Cette fois,
392 ce paramètre a pour valeur d'ébauche un vecteur de 0.
394 Attention, l'application d'une démarche d'adimensionnement nécessite aussi la
395 modification des covariances d'erreurs associées dans la formulation globale du
396 problème d'optimisation.
398 Une telle démarche suffit rarement à éviter tous les problèmes numériques, mais
399 permet souvent d'améliorer beaucoup le conditionnement numérique de
402 Gestion explicite de fonctions "multiples"
403 ++++++++++++++++++++++++++++++++++++++++++
407 il est fortement recommandé de ne pas utiliser cette gestion explicite de
408 fonctions "multiples" sans une très solide raison informatique pour le faire.
409 Cette gestion est déjà effectuée par défaut dans ADAO pour l'amélioration des
410 performances. Seul l'utilisateur très averti, cherchant à gérer des cas
411 particulièrement difficiles, peut s'intéresser à cette extension. En dépit de
412 sa simplicité, c'est au risque explicite de dégrader notablement les
415 Il est possible, lorsque l'on fournit des fonctions d'opérateurs, de les
416 définir comme des fonctions qui traitent non pas un seul argument, mais une
417 série d'arguments, pour restituer en sortie la série des valeurs
418 correspondantes. En pseudo-code, la fonction "multiple", ici nommée
419 ``MultiFunctionO``, représentant l'opérateur classique :math:`O` nommé
420 "*DirectOperator*", effectue::
422 def MultiFunctionO( Inputs ):
426 Y = DirectOperator( X )
430 La longueur de la sortie (c'est-à-dire le nombre de valeurs calculées) est
431 égale à la longueur de l'entrée (c'est-à-dire le nombre d'états dont on veut
432 calculer la valeur par l'opérateur).
434 Cette possibilité n'est disponible que dans l'interface textuelle d'ADAO. Pour
435 cela, lors de la définition d'une fonction d'opérateur, en même temps que l'on
436 définit de manière habituelle la fonction ou le script externe, il suffit
437 d'indiquer en plus en argument par un booléen "*InputFunctionAsMulti*" que la
438 définition est celle d'une fonction "multiple".