2 Copyright (C) 2008-2023 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 La disponibilité des opérateurs d'observation et parfois d'évolution sont
34 nécessaires pour mettre en oeuvre les procédures d'assimilation de données ou
35 d'optimisation. Comme l'opérateur d'évolution est considéré dans sa forme
36 incrémentale, qui représente la transition entre deux états successifs, il est
37 alors formellement similaire à l'opérateur d'observation et la manière de les
40 Ces opérateurs comprennent la **simulation physique par des calculs
41 numériques**. Mais ils comprennent aussi **le filtrage, la projection ou la
42 restriction** des grandeurs simulées, qui sont nécessaires pour comparer la
43 simulation à l'observation.
45 Schématiquement, un opérateur :math:`O` a pour objet de restituer une
46 simulation ou une solution pour des paramètres d'entrée spécifiés. Une partie
47 des paramètres d'entrée peut être modifiée au cours de la procédure
48 d'optimisation. Ainsi, la représentation mathématique d'un tel processus est
49 une fonction. Il a été brièvement décrit dans la section :ref:`section_theory`.
50 Il est généralisé ici par la relation:
52 .. math:: \mathbf{y} = O( \mathbf{x} )
54 entre les pseudo-observations en sortie :math:`\mathbf{y}` et les paramètres
55 d'entrée :math:`\mathbf{x}` en utilisant l'opérateur :math:`O` d'observation ou
56 d'évolution. La même représentation fonctionnelle peut être utilisée
57 pour le modèle linéaire tangent :math:`\mathbf{O}` de :math:`O` et son adjoint
58 :math:`\mathbf{O}^*` qui sont aussi requis par certains algorithmes
59 d'assimilation de données ou d'optimisation.
61 En entrée et en sortie de ces opérateurs, les variables :math:`\mathbf{x}` et
62 :math:`\mathbf{y}`, ou leurs incréments, sont mathématiquement des vecteurs, et
63 ils peuvent donc être donnés par l'utilisateur comme des vecteurs non-orientés
64 (de type liste ou vecteur Numpy) ou orientés (de type matrice Numpy).
66 Ainsi, **pour décrire de manière complète un opérateur, l'utilisateur n'a qu'à
67 fournir une fonction qui réalise complètement et uniquement l'opération
70 Cette fonction est généralement donnée comme une **fonction ou un script
71 Python**, qui peuvent en particulier être exécuté comme une fonction Python
72 indépendante ou dans un noeud YACS. Cette fonction ou ce script peuvent, sans
73 différences, lancer des codes externes ou utiliser des appels et des méthodes
74 internes Python ou SALOME. Si l'algorithme nécessite les 3 aspects de
75 l'opérateur (forme directe, forme tangente et forme adjointe), l'utilisateur
76 doit donner les 3 fonctions ou les approximer grâce à ADAO.
78 Il existe pour l'utilisateur 3 méthodes effectives de fournir une représentation
79 fonctionnelle de l'opérateur, qui diffèrent selon le type d'argument choisi:
81 - :ref:`section_ref_operator_one`
82 - :ref:`section_ref_operator_funcs`
83 - :ref:`section_ref_operator_switch`
85 Dans le cas de l'interface textuelle d'ADAO (TUI), seules les deux premières
86 sont nécessaires car la troisième est incluse dans la seconde. Dans le cas de
87 l'interface graphique EFICAS d'ADAO, ces méthodes sont choisies dans le champ
88 "*FROM*" de chaque opérateur ayant une valeur "*Function*" comme
89 "*INPUT_TYPE*", comme le montre la figure suivante :
91 .. eficas_operator_function:
92 .. image:: images/eficas_operator_function.png
96 **Choisir graphiquement une représentation fonctionnelle de l'opérateur**
98 En interface textuelle d'ADAO (TUI), dans le cas précis illustré ci-dessus, on
99 réalise la même démarche en écrivant :
103 case.set( 'ObservationOperator',
105 Script = 'scripts_for_JDC.py'
109 .. _section_ref_operator_one:
111 Première forme fonctionnelle : un seul opérateur direct
112 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
114 .. index:: single: OneFunction
115 .. index:: single: ScriptWithOneFunction
116 .. index:: single: DirectOperator
117 .. index:: single: DifferentialIncrement
118 .. index:: single: CenteredFiniteDifference
120 La première consiste à ne fournir qu'une seule fonction, potentiellement non
121 linéaire, et à approximer les opérateurs tangent et adjoint associés.
123 Ceci est fait dans ADAO en utilisant, dans l'interface graphique EFICAS d'ADAO,
124 le mot-clé "*ScriptWithOneFunction*" pour la description par un script. Dans
125 l'interface textuelle, c'est le mot-clé "*OneFunction*", éventuellement combiné
126 avec le mot-clé "*Script*" selon que c'est une fonction ou un script. Si c'est
127 par script externe, l'utilisateur doit fournir un fichier contenant une
128 fonction qui porte le nom obligatoire "*DirectOperator*". Par exemple, un
129 script externe peut suivre le modèle générique suivant::
131 def DirectOperator( X ):
132 """ Opérateur direct de simulation non-linéaire """
136 # Résultat : Y = O(X)
137 return "un vecteur similaire à Y"
139 Dans ce cas, l'utilisateur doit aussi fournir une valeur pour l'incrément
140 différentiel ou conserver la valeur par défaut. Cela se réalise en utilisant
141 dans l'interface graphique (GUI) ou textuelle (TUI) le mot-clé
142 "*DifferentialIncrement*" comme paramètre, qui a une valeur par défaut de 1%.
143 Ce coefficient est utilisé dans l'approximation différences finies pour
144 construire les opérateurs tangent et adjoint. L'ordre de l'approximation
145 différences finies peut aussi être choisi à travers l'interface, en utilisant
146 le mot-clé "*CenteredFiniteDifference*", avec ``False`` ou 0 pour un schéma non
147 centré du premier ordre (qui est la valeur par défaut), et avec ``True`` ou 1
148 pour un schéma centré du second ordre (et qui coûte numériquement deux fois
149 plus cher que le premier ordre). Si nécessaire et si possible, on peut
150 :ref:`subsection_ref_parallel_df`. Dans tous les cas, un mécanisme de cache
151 interne permet de limiter le nombre d'évaluations de l'opérateur au minimum
152 possible du point de vue de l'exécution séquentielle ou parallèle des
153 approximations numériques des opérateurs tangent et adjoint, pour éviter des
154 calculs redondants. On se reportera à la partie permettant de
155 :ref:`subsection_iterative_convergence_control` pour connaître l'interaction
156 avec les paramètres relatifs à la convergence.
158 Cette première forme de définition de l'opérateur permet aisément de tester la
159 forme fonctionnelle avant son usage dans un cas ADAO, réduisant notablement la
160 complexité de l'implémentation de l'opérateur. On peut ainsi utiliser
161 l'algorithme ADAO de vérification "*FunctionTest*" (voir la section sur
162 l':ref:`section_ref_algorithm_FunctionTest`) spécifiquement prévu pour ce test.
164 **Important :** le nom "*DirectOperator*" est obligatoire lorsque l'on utilise
165 un script Python indépendant. Le type de l'argument ``X`` en entrée peut être
166 une liste de valeurs réelles, un vecteur Numpy ou une matrice Numpy, et la
167 fonction utilisateur doit accepter et traiter tous ces cas. Le type de
168 l'argument ``Y`` en sortie doit aussi être équivalent à une liste de valeurs
171 Des formes variées d'opérateurs sont disponibles dans les divers scripts inclus
172 dans les :ref:`section_docu_examples`.
174 .. _section_ref_operator_funcs:
176 Seconde forme fonctionnelle : trois opérateurs direct, tangent et adjoint
177 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
179 .. index:: single: ThreeFunctions
180 .. index:: single: ScriptWithFunctions
181 .. index:: single: DirectOperator
182 .. index:: single: TangentOperator
183 .. index:: single: AdjointOperator
187 En général, il est recommandé d'utiliser la première forme fonctionnelle
188 plutôt que la seconde. Un petit accroissement de performances n'est pas une
189 bonne raison pour utiliser l'implémentation détaillée de cette seconde forme
192 La seconde consiste à fournir directement les trois opérateurs liés :math:`O`,
193 :math:`\mathbf{O}` et :math:`\mathbf{O}^*`. C'est effectué en utilisant le
194 mot-clé "*ScriptWithFunctions*" pour la description de l'opérateur choisi dans
195 l'interface graphique EFICAS d'ADAO. Dans l'interface textuelle, c'est le
196 mot-clé "*ThreeFunctions*", éventuellement combiné avec le mot-clé "*Script*"
197 selon que c'est une fonction ou un script. L'utilisateur doit fournir dans un
198 script trois fonctions, avec les trois noms obligatoires "*DirectOperator*",
199 "*TangentOperator*" et "*AdjointOperator*". Par exemple, le script externe peut
200 suivre le squelette suivant::
202 def DirectOperator( X ):
203 """ Opérateur direct de simulation non-linéaire """
207 return "un vecteur similaire à Y"
209 def TangentOperator( paire = (X, dX) ):
210 """ Opérateur linéaire tangent, autour de X, appliqué à dX """
215 return "un vecteur similaire à Y"
217 def AdjointOperator( paire = (X, Y) ):
218 """ Opérateur adjoint, autour de X, appliqué à Y """
223 return "un vecteur similaire à X"
225 Un nouvelle fois, cette seconde définition d'opérateur permet aisément de tester
226 les formes fonctionnelles avant de les utiliser dans le cas ADAO, réduisant la
227 complexité de l'implémentation de l'opérateur.
229 Pour certains algorithmes (en particulier les filtres non ensemblistes), il
230 faut que les fonctions tangente et adjointe puisse renvoyer les matrices
231 équivalentes à l'opérateur linéaire. Dans ce cas, lorsque, respectivement, les
232 arguments ``dX`` ou ``Y`` valent ``None``, le script de l'utilisateur doit
233 renvoyer la matrice associée. Les squelettes des fonctions "*TangentOperator*"
234 et "*AdjointOperator*" deviennent alors les suivants::
236 def TangentOperator( paire = (X, dX) ):
237 """ Opérateur linéaire tangent, autour de X, appliqué à dX """
242 if dX is None or len(dX) == 0:
243 return "la matrice de l'opérateur linéaire tangent"
245 return "un vecteur similaire à Y"
247 def AdjointOperator( paire = (X, Y) ):
248 """ Opérateur adjoint, autour de X, appliqué à Y """
253 if Y is None or len(Y) == 0:
254 return "la matrice de l'opérateur linéaire adjoint"
256 return "un vecteur similaire à X"
258 **Important :** les noms "*DirectOperator*", "*TangentOperator*" et
259 "*AdjointOperator*" sont obligatoires lorsque l'on utilise un script Python
260 indépendant. Le type des arguments en entrée ou en sortie ``X``, ``Y``, ``dX``
261 peut être une liste de valeur réelles, un vecteur Numpy ou une matrice Numpy.
262 La fonction utilisateur doit accepter et traiter tous ces cas.
264 .. _section_ref_operator_switch:
266 Troisième forme fonctionnelle : trois opérateurs avec un branchement
267 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
269 .. index:: single: ScriptWithSwitch
270 .. index:: single: DirectOperator
271 .. index:: single: TangentOperator
272 .. index:: single: AdjointOperator
276 Il est recommandé de ne pas utiliser cette troisième forme fonctionnelle sans
277 une solide raison numérique ou physique. Un accroissement de performances
278 n'est pas une bonne raison pour utiliser la complexité de cette troisième
279 forme fonctionnelle. Seule une impossibilité à utiliser les première ou
280 seconde formes justifie l'usage de la troisième.
282 La troisième forme donne de plus grandes possibilités de contrôle de
283 l'exécution des trois fonctions représentant l'opérateur, permettant un usage
284 et un contrôle avancés sur chaque exécution du code de simulation. C'est
285 réalisable en utilisant le mot-clé "*ScriptWithSwitch*" pour la description de
286 l'opérateur à travers l'interface graphique EFICAS d'ADAO. Dans l'interface
287 textuelle, il suffit d'utiliser le mot-clé "*ThreeFunctions*" précédent pour
288 définir aussi ce cas, en indiquant les fonctions adéquates. L'utilisateur doit
289 fournir un script unique aiguillant, selon un contrôle, l'exécution des formes
290 directe, tangente et adjointe du code de simulation. L'utilisateur peut alors,
291 par exemple, utiliser des approximations pour les codes tangent et adjoint, ou
292 introduire une plus grande complexité du traitement des arguments des
293 fonctions. Mais cette démarche sera plus difficile à implémenter et à déboguer.
295 Toutefois, si vous souhaitez utiliser cette troisième forme, on recommande de
296 se baser sur le modèle suivant pour le script d'aiguillage. Il nécessite un
297 fichier script ou un code externe nommé ici
298 "*Physical_simulation_functions.py*", contenant trois fonctions nommées
299 "*DirectOperator*", "*TangentOperator*" et "*AdjointOperator*" comme
300 précédemment. Voici le squelette d'aiguillage:
303 import Physical_simulation_functions
304 import numpy, logging, codecs, pickle
306 return pickle.loads(codecs.decode(data.encode(), "base64"))
309 for param in computation["specificParameters"]:
310 if param["name"] == "method":
311 method = loads(param["value"])
312 if method not in ["Direct", "Tangent", "Adjoint"]:
313 raise ValueError("No valid computation method is given")
314 logging.info("Found method is \'%s\'"%method)
316 logging.info("Loading operator functions")
317 Function = Physical_simulation_functions.DirectOperator
318 Tangent = Physical_simulation_functions.TangentOperator
319 Adjoint = Physical_simulation_functions.AdjointOperator
321 logging.info("Executing the possible computations")
323 if method == "Direct":
324 logging.info("Direct computation")
325 Xcurrent = computation["inputValues"][0][0][0]
326 data = Function(numpy.matrix( Xcurrent ).T)
327 if method == "Tangent":
328 logging.info("Tangent computation")
329 Xcurrent = computation["inputValues"][0][0][0]
330 dXcurrent = computation["inputValues"][0][0][1]
331 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
332 if method == "Adjoint":
333 logging.info("Adjoint computation")
334 Xcurrent = computation["inputValues"][0][0][0]
335 Ycurrent = computation["inputValues"][0][0][1]
336 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
338 logging.info("Formatting the output")
339 it = numpy.ravel(data)
340 outputValues = [[[[]]]]
342 outputValues[0][0][0].append(val)
345 result["outputValues"] = outputValues
346 result["specificOutputInfos"] = []
347 result["returnCode"] = 0
348 result["errorMessage"] = ""
350 Toutes les modifications envisageables peuvent être faites à partir de cette
351 hypothèse de squelette.
353 .. _section_ref_operator_control:
355 Cas spécial d'un opérateur d'évolution avec contrôle
356 ++++++++++++++++++++++++++++++++++++++++++++++++++++
358 Dans certains cas, l'opérateur d'évolution ou d'observation doit être contrôlé
359 par un contrôle d'entrée externe, qui est donné *a priori*. Dans ce cas, la
360 forme générique du modèle incrémental :math:`O` est légèrement modifiée comme
363 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
365 où :math:`\mathbf{u}` est le contrôle sur l'incrément d'état. En effet,
366 l'opérateur direct doit être appliqué à une paire de variables :math:`(X,U)`.
367 Schématiquement, l'opérateur :math:`O` doit être construit comme une fonction
368 applicable sur une paire :math:`\mathbf{(X, U)}` comme suit :
371 def DirectOperator( paire = (X, U) ):
372 """ Opérateur direct de simulation non-linéaire """
377 return quelque chose comme X(n+1) (évolution) ou Y(n+1) (observation)
379 Les opérateurs tangent et adjoint ont la même signature que précédemment, en
380 notant que les dérivées doivent être faites seulement partiellement par rapport
381 à :math:`\mathbf{x}`. Dans un tel cas de contrôle explicite, seule la deuxième
382 forme fonctionnelle (en utilisant "*ScriptWithFunctions*") et la troisième
383 forme fonctionnelle (en utilisant "*ScriptWithSwitch*") peuvent être utilisées.
385 .. _section_ref_operator_dimensionless:
387 Remarques complémentaires sur l'adimensionnement des opérateurs
388 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
390 .. index:: single: Adimensionnement
391 .. index:: single: Sans dimension
393 Il est fréquent que les grandeurs physiques, en entrée ou en sortie des
394 opérateurs, présentent des différences notables d'ordre de grandeur ou de taux
395 de variation. Une manière d'éviter des difficultés numériques est d'utiliser,
396 ou d'établir, un adimensionnement des calculs menés dans les opérateurs
397 [WikipediaND]_. Par principe, dans la mesure où la simulation de la physique
398 devrait être la plus adimensionnée possible, il est en premier lieu recommandé
399 d'utiliser les capacités existantes d'adimensionnement du code de calcul.
401 Néanmoins, dans le cas courant où l'on ne peut en disposer, il est souvent
402 utile d'environner le calcul pour l'adimensionner en entrée ou en sortie. Une
403 manière simple de faire cela en entrée consiste à transformer les paramètres
404 :math:`\mathbf{x}` en argument d'une fonction comme "*DirectOperator*". On
405 utilise le plus souvent comme référence les valeurs par défaut
406 :math:`\mathbf{x}^b` (ébauche, ou valeur nominale). Pourvu que chaque
407 composante de :math:`\mathbf{x}^b` soit non nulle, on peut ensuite procéder par
408 correction multiplicative. Pour cela, on peut par exemple poser :
410 .. math:: \mathbf{x} = \mathbf{\alpha}\mathbf{x}^b
412 et optimiser ensuite le paramètre multiplicatif :math:`\mathbf{\alpha}`. Ce
413 paramètre a pour valeur par défaut (ou pour ébauche) un vecteur de 1. De
414 manière similaire, on peut procéder par correction additive si c'est plus
415 judicieux pour la physique sous-jacente. Ainsi, dans ce cas, on peut poser :
417 .. math:: \mathbf{x} =\mathbf{x}^b + \mathbf{\alpha}
419 et optimiser ensuite le paramètre additif :math:`\mathbf{\alpha}`. Cette fois,
420 ce paramètre a pour valeur d'ébauche un vecteur de 0.
422 Attention, l'application d'une démarche d'adimensionnement nécessite aussi la
423 modification des covariances d'erreurs associées dans la formulation globale du
424 problème d'optimisation.
426 Une telle démarche suffit rarement à éviter tous les problèmes numériques, mais
427 permet souvent d'améliorer beaucoup le conditionnement numérique de
430 .. index:: single: InputFunctionAsMulti
432 Gestion explicite de fonctions "multiples"
433 ++++++++++++++++++++++++++++++++++++++++++
437 Il est fortement recommandé de ne pas utiliser cette gestion explicite de
438 fonctions "multiples" sans une très solide raison informatique pour le faire.
439 Cette gestion est déjà effectuée par défaut dans ADAO pour l'amélioration des
440 performances. Seul l'utilisateur très averti, cherchant à gérer des cas
441 particulièrement difficiles, peut s'intéresser à cette extension. En dépit de
442 sa simplicité, c'est au risque explicite de dégrader notablement les
443 performances, ou d'avoir des erreurs d'exécution étranges.
445 Il est possible, lorsque l'on fournit des fonctions d'opérateurs, de les
446 définir comme des fonctions qui traitent non pas un seul argument, mais une
447 série d'arguments, pour restituer en sortie la série des valeurs
448 correspondantes. En pseudo-code, la fonction "multiple", ici nommée
449 ``MultiFunctionO``, représentant l'opérateur classique :math:`O` nommé
450 "*DirectOperator*", effectue :
453 def MultiFunctionO( Inputs ):
457 Y = DirectOperator( X )
461 La longueur de la sortie (c'est-à-dire le nombre de valeurs calculées) est
462 égale à la longueur de l'entrée (c'est-à-dire le nombre d'états dont on veut
463 calculer la valeur par l'opérateur).
465 Cette possibilité n'est disponible que dans l'interface textuelle TUI d'ADAO.
466 Pour cela, lors de la définition d'une fonction d'opérateur, en même temps que
467 l'on définit de manière habituelle la fonction ou le script externe, il suffit
468 d'indiquer en plus en argument par un booléen supplémentaire
469 "*InputFunctionAsMulti*" que la définition est celle d'une fonction "multiple".
470 Par exemple, si c'est l'opérateur d'observation que l'on définit de cette
471 manière, il faut écrire (sachant que toutes les autres commandes optionnelles
472 restent inchangées) :
475 case.set( 'ObservationOperator',
476 OneFunction = MultiFunctionO,
478 InputFunctionAsMulti = True,