2 Copyright (C) 2008-2014 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 Requirements for functions describing an operator
27 -------------------------------------------------
29 The operators for observation and evolution are required to implement the data
30 assimilation or optimization procedures. They include the physical simulation by
31 numerical calculations, but also the filtering and restriction to compare the
32 simulation to observation. The evolution operator is considered here in its
33 incremental form, representing the transition between two successive states, and
34 is then similar to the observation operator.
36 Schematically, an operator has to give a output solution given the input
37 parameters. Part of the input parameters can be modified during the optimization
38 procedure. So the mathematical representation of such a process is a function.
39 It was briefly described in the section :ref:`section_theory` and is generalized
42 .. math:: \mathbf{y} = O( \mathbf{x} )
44 between the pseudo-observations :math:`\mathbf{y}` and the parameters
45 :math:`\mathbf{x}` using the observation or evolution operator :math:`O`. The
46 same functional representation can be used for the linear tangent model
47 :math:`\mathbf{O}` of :math:`O` and its adjoint :math:`\mathbf{O}^*`, also
48 required by some data assimilation or optimization algorithms.
50 On input and output of these operators, the :math:`\mathbf{x}` and
51 :math:`\mathbf{y}` variables or their increments are mathematically vectors,
52 and they are given as non-oriented vectors (of type list or Numpy array) or
53 oriented ones (of type Numpy matrix).
55 Then, **to describe completely an operator, the user has only to provide a
56 function that fully and only realize the functional operation**.
58 This function is usually given as a script that can be executed in a YACS node.
59 This script can without difference launch external codes or use internal SALOME
60 calls and methods. If the algorithm requires the 3 aspects of the operator
61 (direct form, tangent form and adjoint form), the user has to give the 3
62 functions or to approximate them.
64 There are 3 practical methods for the user to provide an operator functional
65 representation. These methods are chosen in the "*FROM*" field of each operator
66 having a "*Function*" value as "*INPUT_TYPE*", as shown by the following figure:
68 .. eficas_operator_function:
69 .. image:: images/eficas_operator_function.png
73 **Choosing an operator functional representation**
75 First functional form: using "*ScriptWithOneFunction*"
76 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
78 .. index:: single: ScriptWithOneFunction
79 .. index:: single: DirectOperator
80 .. index:: single: DifferentialIncrement
81 .. index:: single: CenteredFiniteDifference
83 The first one consist in providing only one potentially non-linear function, and
84 to approximate the tangent and the adjoint operators. This is done by using the
85 keyword "*ScriptWithOneFunction*" for the description of the chosen operator in
86 the ADAO GUI. The user have to provide the function in a script, with a
87 mandatory name "*DirectOperator*". For example, the script can follow the
90 def DirectOperator( X ):
91 """ Direct non-linear simulation operator """
97 In this case, the user has also provide a value for the differential increment
98 (or keep the default value), using through the GUI the keyword
99 "*DifferentialIncrement*", which has a default value of 1%. This coefficient
100 will be used in the finite differences approximation to build the tangent and
101 adjoint operators. The finite differences approximation order can also be chosen
102 through the GUI, using the keyword "*CenteredFiniteDifference*", with 0 for an
103 uncentered schema of first order (which is the default value), and with 1 for a
104 centered schema of second order (of twice the first order computational cost).
106 This first operator definition form allows easily to test the functional form
107 before its use in an ADAO case, greatly reducing the complexity of
108 operator implementation.
110 **Important warning:** the name "*DirectOperator*" is mandatory, and the type of
111 the ``X`` argument can be either a list, a numpy array or a numpy 1D-matrix. The
112 user has to treat these cases in his function.
114 Second functional form: using "*ScriptWithFunctions*"
115 +++++++++++++++++++++++++++++++++++++++++++++++++++++
117 .. index:: single: ScriptWithFunctions
118 .. index:: single: DirectOperator
119 .. index:: single: TangentOperator
120 .. index:: single: AdjointOperator
122 **In general, it is recommended to use the first functional form rather than
123 the second one. A small performance improvement is not a good reason to use a
124 detailed implementation as this second functional form.**
126 The second one consist in providing directly the three associated operators
127 :math:`O`, :math:`\mathbf{O}` and :math:`\mathbf{O}^*`. This is done by using
128 the keyword "*ScriptWithFunctions*" for the description of the chosen operator
129 in the ADAO GUI. The user have to provide three functions in one script, with
130 three mandatory names "*DirectOperator*", "*TangentOperator*" and
131 "*AdjointOperator*". For example, the script can follow the template::
133 def DirectOperator( X ):
134 """ Direct non-linear simulation operator """
138 return something like Y
140 def TangentOperator( (X, dX) ):
141 """ Tangent linear operator, around X, applied to dX """
145 return something like Y
147 def AdjointOperator( (X, Y) ):
148 """ Adjoint operator, around X, applied to Y """
152 return something like X
154 Another time, this second operator definition allow easily to test the
155 functional forms before their use in an ADAO case, reducing the complexity of
156 operator implementation.
158 For some algorithms, it is required that the tangent and adjoint functions can
159 return the matrix equivalent to the linear operator. In this case, when
160 respectively the ``dX`` or the ``Y`` arguments are ``None``, the user has to
161 return the associated matrix.
163 **Important warning:** the names "*DirectOperator*", "*TangentOperator*" and
164 "*AdjointOperator*" are mandatory, and the type of the ``X``, Y``, ``dX``
165 arguments can be either a python list, a numpy array or a numpy 1D-matrix. The
166 user has to treat these cases in his script.
168 Third functional form: using "*ScriptWithSwitch*"
169 +++++++++++++++++++++++++++++++++++++++++++++++++
171 .. index:: single: ScriptWithSwitch
172 .. index:: single: DirectOperator
173 .. index:: single: TangentOperator
174 .. index:: single: AdjointOperator
176 **It is recommended not to use this third functional form without a solid
177 numerical or physical reason. A performance improvement is not a good reason to
178 use the implementation complexity of this third functional form. Only an
179 inability to use the first or second forms justifies the use of the third.**
181 This third form give more possibilities to control the execution of the three
182 functions representing the operator, allowing advanced usage and control over
183 each execution of the simulation code. This is done by using the keyword
184 "*ScriptWithSwitch*" for the description of the chosen operator in the ADAO GUI.
185 The user have to provide a switch in one script to control the execution of the
186 direct, tangent and adjoint forms of its simulation code. The user can then, for
187 example, use other approximations for the tangent and adjoint codes, or
188 introduce more complexity in the argument treatment of the functions. But it
189 will be far more complicated to implement and debug.
191 If, however, you want to use this third form, we recommend using the following
192 template for the switch. It requires an external script or code named here
193 "*Physical_simulation_functions.py*", containing three functions named
194 "*DirectOperator*", "*TangentOperator*" and "*AdjointOperator*" as previously.
195 Here is the switch template::
197 import Physical_simulation_functions
198 import numpy, logging
201 for param in computation["specificParameters"]:
202 if param["name"] == "method":
203 method = param["value"]
204 if method not in ["Direct", "Tangent", "Adjoint"]:
205 raise ValueError("No valid computation method is given")
206 logging.info("Found method is \'%s\'"%method)
208 logging.info("Loading operator functions")
209 Function = Physical_simulation_functions.DirectOperator
210 Tangent = Physical_simulation_functions.TangentOperator
211 Adjoint = Physical_simulation_functions.AdjointOperator
213 logging.info("Executing the possible computations")
215 if method == "Direct":
216 logging.info("Direct computation")
217 Xcurrent = computation["inputValues"][0][0][0]
218 data = Function(numpy.matrix( Xcurrent ).T)
219 if method == "Tangent":
220 logging.info("Tangent computation")
221 Xcurrent = computation["inputValues"][0][0][0]
222 dXcurrent = computation["inputValues"][0][0][1]
223 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
224 if method == "Adjoint":
225 logging.info("Adjoint computation")
226 Xcurrent = computation["inputValues"][0][0][0]
227 Ycurrent = computation["inputValues"][0][0][1]
228 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
230 logging.info("Formatting the output")
231 it = numpy.ravel(data)
232 outputValues = [[[[]]]]
234 outputValues[0][0][0].append(val)
237 result["outputValues"] = outputValues
238 result["specificOutputInfos"] = []
239 result["returnCode"] = 0
240 result["errorMessage"] = ""
242 All various modifications could be done from this template hypothesis.
244 Special case of controlled evolution or observation operator
245 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
247 In some cases, the evolution or the observation operator is required to be
248 controlled by an external input control, given *a priori*. In this case, the
249 generic form of the incremental model is slightly modified as follows:
251 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
253 where :math:`\mathbf{u}` is the control over one state increment. In this case,
254 the direct operator has to be applied to a pair of variables :math:`(X,U)`.
255 Schematically, the operator has to be set as::
257 def DirectOperator( (X, U) ):
258 """ Direct non-linear simulation operator """
262 return something like X(n+1) (evolution) or Y(n+1) (observation)
264 The tangent and adjoint operators have the same signature as previously, noting
265 that the derivatives has to be done only partially against :math:`\mathbf{x}`.
266 In such a case with explicit control, only the second functional form (using
267 "*ScriptWithFunctions*") and third functional form (using "*ScriptWithSwitch*")