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 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::
287 import Physical_simulation_functions
288 import numpy, logging, codecs, pickle
290 return pickle.loads(codecs.decode(data.encode(), "base64"))
293 for param in computation["specificParameters"]:
294 if param["name"] == "method":
295 method = loads(param["value"])
296 if method not in ["Direct", "Tangent", "Adjoint"]:
297 raise ValueError("No valid computation method is given")
298 logging.info("Found method is \'%s\'"%method)
300 logging.info("Loading operator functions")
301 Function = Physical_simulation_functions.DirectOperator
302 Tangent = Physical_simulation_functions.TangentOperator
303 Adjoint = Physical_simulation_functions.AdjointOperator
305 logging.info("Executing the possible computations")
307 if method == "Direct":
308 logging.info("Direct computation")
309 Xcurrent = computation["inputValues"][0][0][0]
310 data = Function(numpy.matrix( Xcurrent ).T)
311 if method == "Tangent":
312 logging.info("Tangent computation")
313 Xcurrent = computation["inputValues"][0][0][0]
314 dXcurrent = computation["inputValues"][0][0][1]
315 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
316 if method == "Adjoint":
317 logging.info("Adjoint computation")
318 Xcurrent = computation["inputValues"][0][0][0]
319 Ycurrent = computation["inputValues"][0][0][1]
320 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
322 logging.info("Formatting the output")
323 it = numpy.ravel(data)
324 outputValues = [[[[]]]]
326 outputValues[0][0][0].append(val)
329 result["outputValues"] = outputValues
330 result["specificOutputInfos"] = []
331 result["returnCode"] = 0
332 result["errorMessage"] = ""
334 Toutes les modifications envisageables peuvent être faites à partir de cette
335 hypothèse de squelette.
337 .. _section_ref_operator_control:
339 Cas spécial d'un opérateur d'évolution avec contrôle
340 ++++++++++++++++++++++++++++++++++++++++++++++++++++
342 Dans certains cas, l'opérateur d'évolution ou d'observation doit être contrôlé
343 par un contrôle d'entrée externe, qui est donné *a priori*. Dans ce cas, la
344 forme générique du modèle incrémental :math:`O` est légèrement modifiée comme
347 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
349 où :math:`\mathbf{u}` est le contrôle sur l'incrément d'état. En effet,
350 l'opérateur direct doit être appliqué à une paire de variables :math:`(X,U)`.
351 Schématiquement, l'opérateur :math:`O` doit être construit comme une fonction
352 applicable sur une paire:math:`\mathbf{(X, U)}` suit::
354 def DirectOperator( paire = (X, U) ):
355 """ Opérateur direct de simulation non-linéaire """
360 return quelque chose comme X(n+1) (évolution) ou Y(n+1) (observation)
362 Les opérateurs tangent et adjoint ont la même signature que précédemment, en
363 notant que les dérivées doivent être faites seulement partiellement par rapport
364 à :math:`\mathbf{x}`. Dans un tel cas de contrôle explicite, seule la deuxième
365 forme fonctionnelle (en utilisant "*ScriptWithFunctions*") et la troisième
366 forme fonctionnelle (en utilisant "*ScriptWithSwitch*") peuvent être utilisées.
368 .. _section_ref_operator_dimensionless:
370 Remarques complémentaires sur l'adimensionnement des opérateurs
371 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
373 .. index:: single: Adimensionnement
374 .. index:: single: Sans dimension
376 Il est fréquent que les grandeurs physiques, en entrée ou en sortie des
377 opérateurs, présentent des différences notables d'ordre de grandeur ou de taux
378 de variation. Une manière d'éviter des difficultés numériques est d'utiliser,
379 ou d'établir, un adimensionnement des calculs menés dans les opérateurs
380 [WikipediaND]_. Par principe, dans la mesure où la simulation de la physique
381 devrait être la plus adimensionnée possible, il est en premier lieu recommandé
382 d'utiliser les capacités existantes d'adimensionnement du code de calcul.
384 Néanmoins, dans le cas courant où l'on ne peut en disposer, il est souvent
385 utile d'environner le calcul pour l'adimensionner en entrée ou en sortie. Une
386 manière simple de faire cela en entrée consiste à transformer les paramètres
387 :math:`\mathbf{x}` en argument d'une fonction comme "*DirectOperator*". On
388 utilise le plus souvent comme référence les valeurs par défaut
389 :math:`\mathbf{x}^b` (ébauche, ou valeur nominale). Pourvu que chaque
390 composante de :math:`\mathbf{x}^b` soit non nulle, on peut ensuite procéder par
391 correction multiplicative. Pour cela, on peut par exemple poser:
393 .. math:: \mathbf{x} = \mathbf{\alpha}\mathbf{x}^b
395 et optimiser ensuite le paramètre multiplicatif :math:`\mathbf{\alpha}`. Ce
396 paramètre a pour valeur par défaut (ou pour ébauche) un vecteur de 1. De
397 manière similaire, on peut procéder par correction additive si c'est plus
398 judicieux pour la physique sous-jacente. Ainsi, dans ce cas, on peut poser:
400 .. math:: \mathbf{x} =\mathbf{x}^b + \mathbf{\alpha}
402 et optimiser ensuite le paramètre additif :math:`\mathbf{\alpha}`. Cette fois,
403 ce paramètre a pour valeur d'ébauche un vecteur de 0.
405 Attention, l'application d'une démarche d'adimensionnement nécessite aussi la
406 modification des covariances d'erreurs associées dans la formulation globale du
407 problème d'optimisation.
409 Une telle démarche suffit rarement à éviter tous les problèmes numériques, mais
410 permet souvent d'améliorer beaucoup le conditionnement numérique de
413 Gestion explicite de fonctions "multiples"
414 ++++++++++++++++++++++++++++++++++++++++++
418 il est fortement recommandé de ne pas utiliser cette gestion explicite de
419 fonctions "multiples" sans une très solide raison informatique pour le faire.
420 Cette gestion est déjà effectuée par défaut dans ADAO pour l'amélioration des
421 performances. Seul l'utilisateur très averti, cherchant à gérer des cas
422 particulièrement difficiles, peut s'intéresser à cette extension. En dépit de
423 sa simplicité, c'est au risque explicite de dégrader notablement les
426 Il est possible, lorsque l'on fournit des fonctions d'opérateurs, de les
427 définir comme des fonctions qui traitent non pas un seul argument, mais une
428 série d'arguments, pour restituer en sortie la série des valeurs
429 correspondantes. En pseudo-code, la fonction "multiple", ici nommée
430 ``MultiFunctionO``, représentant l'opérateur classique :math:`O` nommé
431 "*DirectOperator*", effectue::
433 def MultiFunctionO( Inputs ):
437 Y = DirectOperator( X )
441 La longueur de la sortie (c'est-à-dire le nombre de valeurs calculées) est
442 égale à la longueur de l'entrée (c'est-à-dire le nombre d'états dont on veut
443 calculer la valeur par l'opérateur).
445 Cette possibilité n'est disponible que dans l'interface textuelle d'ADAO. Pour
446 cela, lors de la définition d'une fonction d'opérateur, en même temps que l'on
447 définit de manière habituelle la fonction ou le script externe, il suffit
448 d'indiquer en plus en argument par un booléen "*InputFunctionAsMulti*" que la
449 définition est celle d'une fonction "multiple".