2 Copyright (C) 2008-2017 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. In
106 all cases, an internal cache mechanism is used to restrict the number of
107 operator evaluations at the minimum possible in a sequential or parallel
108 execution scheme for numerical approximations of the tangent and adjoint
109 operators, to avoid redundant calculations.
111 This first operator definition form allows easily to test the functional form
112 before its use in an ADAO case, greatly reducing the complexity of operator
113 implementation. One can then use the "*FunctionTest*" ADAO checking algorithm
114 (see the section on the :ref:`section_ref_algorithm_FunctionTest`) for this
117 **Important warning:** the name "*DirectOperator*" is mandatory, and the type of
118 the ``X`` argument can be either a list, a Numpy array or a Numpy 1D-matrix. The
119 user function has to accept and treat all these cases.
121 Second functional form: using "*ScriptWithFunctions*"
122 +++++++++++++++++++++++++++++++++++++++++++++++++++++
124 .. index:: single: ScriptWithFunctions
125 .. index:: single: DirectOperator
126 .. index:: single: TangentOperator
127 .. index:: single: AdjointOperator
129 **In general, it is recommended to use the first functional form rather than
130 the second one. A small performance improvement is not a good reason to use a
131 detailed implementation as this second functional form.**
133 The second one consist in providing directly the three associated operators
134 :math:`O`, :math:`\mathbf{O}` and :math:`\mathbf{O}^*`. This is done by using
135 the keyword "*ScriptWithFunctions*" for the description of the chosen operator
136 in the ADAO GUI. The user have to provide three functions in one script, with
137 three mandatory names "*DirectOperator*", "*TangentOperator*" and
138 "*AdjointOperator*". For example, the script can follow the template::
140 def DirectOperator( X ):
141 """ Direct non-linear simulation operator """
145 return something like Y
147 def TangentOperator( (X, dX) ):
148 """ Tangent linear operator, around X, applied to dX """
152 return something like Y
154 def AdjointOperator( (X, Y) ):
155 """ Adjoint operator, around X, applied to Y """
159 return something like X
161 Another time, this second operator definition allow easily to test the
162 functional forms before their use in an ADAO case, reducing the complexity of
163 operator implementation.
165 For some algorithms, it is required that the tangent and adjoint functions can
166 return the matrix equivalent to the linear operator. In this case, when
167 respectively the ``dX`` or the ``Y`` arguments are ``None``, the user has to
168 return the associated matrix.
170 **Important warning:** the names "*DirectOperator*", "*TangentOperator*" and
171 "*AdjointOperator*" are mandatory, and the type of the ``X``, Y``, ``dX``
172 arguments can be either a python list, a Numpy array or a Numpy 1D-matrix. The
173 user has to treat these cases in his script.
175 Third functional form: using "*ScriptWithSwitch*"
176 +++++++++++++++++++++++++++++++++++++++++++++++++
178 .. index:: single: ScriptWithSwitch
179 .. index:: single: DirectOperator
180 .. index:: single: TangentOperator
181 .. index:: single: AdjointOperator
183 **It is recommended not to use this third functional form without a solid
184 numerical or physical reason. A performance improvement is not a good reason to
185 use the implementation complexity of this third functional form. Only an
186 inability to use the first or second forms justifies the use of the third.**
188 This third form give more possibilities to control the execution of the three
189 functions representing the operator, allowing advanced usage and control over
190 each execution of the simulation code. This is done by using the keyword
191 "*ScriptWithSwitch*" for the description of the chosen operator in the ADAO GUI.
192 The user have to provide a switch in one script to control the execution of the
193 direct, tangent and adjoint forms of its simulation code. The user can then, for
194 example, use other approximations for the tangent and adjoint codes, or
195 introduce more complexity in the argument treatment of the functions. But it
196 will be far more complicated to implement and debug.
198 If, however, you want to use this third form, we recommend using the following
199 template for the switch. It requires an external script or code named here
200 "*Physical_simulation_functions.py*", containing three functions named
201 "*DirectOperator*", "*TangentOperator*" and "*AdjointOperator*" as previously.
202 Here is the switch template::
204 import Physical_simulation_functions
205 import numpy, logging
208 for param in computation["specificParameters"]:
209 if param["name"] == "method":
210 method = param["value"]
211 if method not in ["Direct", "Tangent", "Adjoint"]:
212 raise ValueError("No valid computation method is given")
213 logging.info("Found method is \'%s\'"%method)
215 logging.info("Loading operator functions")
216 Function = Physical_simulation_functions.DirectOperator
217 Tangent = Physical_simulation_functions.TangentOperator
218 Adjoint = Physical_simulation_functions.AdjointOperator
220 logging.info("Executing the possible computations")
222 if method == "Direct":
223 logging.info("Direct computation")
224 Xcurrent = computation["inputValues"][0][0][0]
225 data = Function(numpy.matrix( Xcurrent ).T)
226 if method == "Tangent":
227 logging.info("Tangent computation")
228 Xcurrent = computation["inputValues"][0][0][0]
229 dXcurrent = computation["inputValues"][0][0][1]
230 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
231 if method == "Adjoint":
232 logging.info("Adjoint computation")
233 Xcurrent = computation["inputValues"][0][0][0]
234 Ycurrent = computation["inputValues"][0][0][1]
235 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
237 logging.info("Formatting the output")
238 it = numpy.ravel(data)
239 outputValues = [[[[]]]]
241 outputValues[0][0][0].append(val)
244 result["outputValues"] = outputValues
245 result["specificOutputInfos"] = []
246 result["returnCode"] = 0
247 result["errorMessage"] = ""
249 All various modifications could be done from this template hypothesis.
251 .. _section_ref_operator_control:
253 Special case of controled evolution or observation operator
254 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
256 In some cases, the evolution or the observation operator is required to be
257 controlled by an external input control, given *a priori*. In this case, the
258 generic form of the incremental model is slightly modified as follows:
260 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
262 where :math:`\mathbf{u}` is the control over one state increment. In fact, the
263 direct operator has to be applied to a pair of variables :math:`(X,U)`.
264 Schematically, the operator has to be set as::
266 def DirectOperator( (X, U) ):
267 """ Direct non-linear simulation operator """
271 return something like X(n+1) (evolution) or Y(n+1) (observation)
273 The tangent and adjoint operators have the same signature as previously, noting
274 that the derivatives has to be done only partially against :math:`\mathbf{x}`.
275 In such a case with explicit control, only the second functional form (using
276 "*ScriptWithFunctions*") and third functional form (using "*ScriptWithSwitch*")
279 Additional notes on dimensionless transformation of operators
280 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
282 .. index:: single: Nondimensionalization
283 .. index:: single: Dimensionless
285 It is common that physical quantities, in input or output of the operators, have
286 significant differences in magnitude or rate of change. One way to avoid
287 numerical difficulties is to use, or to set, a dimensionless version of
288 calculations carried out in operators [WikipediaND]_. In principle, since
289 physical simulation should be as dimensionless as possible, it is at first
290 recommended to use the existing dimensionless capacity of the calculation code.
292 However, in the common case where we can not dispose of it, it is often useful
293 to surround the calculation to remove dimension for input or output. A simple
294 way to do this is to convert the input parameters :math:`\mathbf{x}` which are
295 arguments of a function like "*DirectOperator*". One mostly use the default
296 values :math:`\mathbf{x}^b` (background, or nominal value). Provided that each
297 component of :math:`\mathbf{x}^b` is non zero, one can indeed put:
299 .. math:: \mathbf{x} = \mathbf{\alpha}\mathbf{x}^b
301 and then optimize the multiplicative parameter :math:`\mathbf{\alpha}`. This
302 parameter has as default value (or as background) a vector of 1. Be careful,
303 applying a process of dimensionless transformation also requires changing the
304 associated error covariances in an ADAO formulation of the optimization problem.
306 Such a process is rarely enough to avoid all the numerical problems, but it
307 often improves a lot the numeric conditioning of the optimization.