2 Copyright (C) 2008-2015 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).
105 If necessary and if possible, :ref:`subsection_ref_parallel_df` can be used.
107 This first operator definition form allows easily to test the functional form
108 before its use in an ADAO case, greatly reducing the complexity of operator
109 implementation. One can then use the "*FunctionTest*" ADAO checking algorithm
110 (see the section on the :ref:`section_ref_algorithm_FunctionTest`) for this
113 **Important warning:** the name "*DirectOperator*" is mandatory, and the type of
114 the ``X`` argument can be either a list, a numpy array or a numpy 1D-matrix. The
115 user function has to accept and treat all these cases.
117 Second functional form: using "*ScriptWithFunctions*"
118 +++++++++++++++++++++++++++++++++++++++++++++++++++++
120 .. index:: single: ScriptWithFunctions
121 .. index:: single: DirectOperator
122 .. index:: single: TangentOperator
123 .. index:: single: AdjointOperator
125 **In general, it is recommended to use the first functional form rather than
126 the second one. A small performance improvement is not a good reason to use a
127 detailed implementation as this second functional form.**
129 The second one consist in providing directly the three associated operators
130 :math:`O`, :math:`\mathbf{O}` and :math:`\mathbf{O}^*`. This is done by using
131 the keyword "*ScriptWithFunctions*" for the description of the chosen operator
132 in the ADAO GUI. The user have to provide three functions in one script, with
133 three mandatory names "*DirectOperator*", "*TangentOperator*" and
134 "*AdjointOperator*". For example, the script can follow the template::
136 def DirectOperator( X ):
137 """ Direct non-linear simulation operator """
141 return something like Y
143 def TangentOperator( (X, dX) ):
144 """ Tangent linear operator, around X, applied to dX """
148 return something like Y
150 def AdjointOperator( (X, Y) ):
151 """ Adjoint operator, around X, applied to Y """
155 return something like X
157 Another time, this second operator definition allow easily to test the
158 functional forms before their use in an ADAO case, reducing the complexity of
159 operator implementation.
161 For some algorithms, it is required that the tangent and adjoint functions can
162 return the matrix equivalent to the linear operator. In this case, when
163 respectively the ``dX`` or the ``Y`` arguments are ``None``, the user has to
164 return the associated matrix.
166 **Important warning:** the names "*DirectOperator*", "*TangentOperator*" and
167 "*AdjointOperator*" are mandatory, and the type of the ``X``, Y``, ``dX``
168 arguments can be either a python list, a numpy array or a numpy 1D-matrix. The
169 user has to treat these cases in his script.
171 Third functional form: using "*ScriptWithSwitch*"
172 +++++++++++++++++++++++++++++++++++++++++++++++++
174 .. index:: single: ScriptWithSwitch
175 .. index:: single: DirectOperator
176 .. index:: single: TangentOperator
177 .. index:: single: AdjointOperator
179 **It is recommended not to use this third functional form without a solid
180 numerical or physical reason. A performance improvement is not a good reason to
181 use the implementation complexity of this third functional form. Only an
182 inability to use the first or second forms justifies the use of the third.**
184 This third form give more possibilities to control the execution of the three
185 functions representing the operator, allowing advanced usage and control over
186 each execution of the simulation code. This is done by using the keyword
187 "*ScriptWithSwitch*" for the description of the chosen operator in the ADAO GUI.
188 The user have to provide a switch in one script to control the execution of the
189 direct, tangent and adjoint forms of its simulation code. The user can then, for
190 example, use other approximations for the tangent and adjoint codes, or
191 introduce more complexity in the argument treatment of the functions. But it
192 will be far more complicated to implement and debug.
194 If, however, you want to use this third form, we recommend using the following
195 template for the switch. It requires an external script or code named here
196 "*Physical_simulation_functions.py*", containing three functions named
197 "*DirectOperator*", "*TangentOperator*" and "*AdjointOperator*" as previously.
198 Here is the switch template::
200 import Physical_simulation_functions
201 import numpy, logging
204 for param in computation["specificParameters"]:
205 if param["name"] == "method":
206 method = param["value"]
207 if method not in ["Direct", "Tangent", "Adjoint"]:
208 raise ValueError("No valid computation method is given")
209 logging.info("Found method is \'%s\'"%method)
211 logging.info("Loading operator functions")
212 Function = Physical_simulation_functions.DirectOperator
213 Tangent = Physical_simulation_functions.TangentOperator
214 Adjoint = Physical_simulation_functions.AdjointOperator
216 logging.info("Executing the possible computations")
218 if method == "Direct":
219 logging.info("Direct computation")
220 Xcurrent = computation["inputValues"][0][0][0]
221 data = Function(numpy.matrix( Xcurrent ).T)
222 if method == "Tangent":
223 logging.info("Tangent computation")
224 Xcurrent = computation["inputValues"][0][0][0]
225 dXcurrent = computation["inputValues"][0][0][1]
226 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
227 if method == "Adjoint":
228 logging.info("Adjoint computation")
229 Xcurrent = computation["inputValues"][0][0][0]
230 Ycurrent = computation["inputValues"][0][0][1]
231 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
233 logging.info("Formatting the output")
234 it = numpy.ravel(data)
235 outputValues = [[[[]]]]
237 outputValues[0][0][0].append(val)
240 result["outputValues"] = outputValues
241 result["specificOutputInfos"] = []
242 result["returnCode"] = 0
243 result["errorMessage"] = ""
245 All various modifications could be done from this template hypothesis.
247 .. _section_ref_operator_control:
249 Special case of controled evolution or observation operator
250 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
252 In some cases, the evolution or the observation operator is required to be
253 controlled by an external input control, given *a priori*. In this case, the
254 generic form of the incremental model is slightly modified as follows:
256 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
258 where :math:`\mathbf{u}` is the control over one state increment. In fact, the
259 direct operator has to be applied to a pair of variables :math:`(X,U)`.
260 Schematically, the operator has to be set as::
262 def DirectOperator( (X, U) ):
263 """ Direct non-linear simulation operator """
267 return something like X(n+1) (evolution) or Y(n+1) (observation)
269 The tangent and adjoint operators have the same signature as previously, noting
270 that the derivatives has to be done only partially against :math:`\mathbf{x}`.
271 In such a case with explicit control, only the second functional form (using
272 "*ScriptWithFunctions*") and third functional form (using "*ScriptWithSwitch*")