3 ================================================================================
4 Tutorials on using the ADAO module
5 ================================================================================
7 .. |eficas_new| image:: images/eficas_new.png
10 .. |eficas_save| image:: images/eficas_save.png
13 .. |eficas_saveas| image:: images/eficas_saveas.png
16 .. |eficas_yacs| image:: images/eficas_yacs.png
20 This section presents some examples on using the ADAO module in SALOME. The
21 first one shows how to build a simple data assimilation case defining
22 explicitly all the required data through the GUI. The second one shows, on the
23 same case, how to define data using external sources through scripts.
25 Building a simple estimation case with explicit data definition
26 ---------------------------------------------------------------
28 This simple example is a demonstration one, and describes how to set a BLUE
29 estimation framework in order to get *ponderated (or fully weighted) least
30 square estimated state* of a system from an observation of the state and from an
31 *a priori* knowledge (or background) of this state. In other words, we look for
32 the weighted middle between the observation and the background vectors. All the
33 numerical values of this example are arbitrary.
38 We choose to operate in a 3-dimensional space. 3D is chosen in order to restrict
39 the size of numerical object to explicitly enter by the user, but the problem is
40 not dependant of the dimension and can be set in dimension 1000... The
41 observation :math:`\mathbf{y}^o` is of value 1 in each direction, so:
45 The background state :math:`\mathbf{x}^b`, which represent some *a priori*
46 knowledge or a regularization, is of value of 0 in each direction, which is:
50 Data assimilation requires information on errors covariances :math:`\mathbf{R}`
51 and :math:`\mathbf{B}` respectively for observation and background variables. We
52 choose here to have uncorrelated errors (that is, diagonal matrices) and to have
53 the same variance of 1 for all variables (that is, identity matrices). We get:
55 ``B = R = [1 0 0 ; 0 1 0 ; 0 0 1]``
57 Last, we need an observation operator :math:`\mathbf{H}` to convert the
58 background value in the space of observation value. Here, because the space
59 dimensions are the same, we can choose the identity as the observation
62 ``H = [1 0 0 ; 0 1 0 ; 0 0 1]``
64 With such choices, the Best Linear Unbiased Estimator (BLUE) will be the average
65 vector between :math:`\mathbf{y}^o` and :math:`\mathbf{x}^b`, named the
66 *analysis* and denoted by :math:`\mathbf{x}^a`:
68 ``Xa = [0.5 0.5 0.5]``
70 As en extension of this example, one can change the variances for
71 :math:`\mathbf{B}` or :math:`\mathbf{R}` independently, and the analysis will
72 move to :math:`\mathbf{y}^o` or to :math:`\mathbf{x}^b` in inverse proportion of
73 the variances in :math:`\mathbf{B}` and :math:`\mathbf{R}`. It is also
74 equivalent to search for the analysis thought a BLUE algorithm or a 3DVAR one.
76 Using the GUI to build the ADAO case
77 ++++++++++++++++++++++++++++++++++++
79 First, you have to activate the ADAO module by choosing the appropriate module
80 button or menu of SALOME, and you will see:
83 .. image:: images/adao_activate.png
87 **Activating the module ADAO in SALOME**
89 Choose the "*New*" button in this window. You will directly get the EFICAS
90 interface for variables definition, along with the "*Object browser*". You can
91 then click on the "*New*" button |eficas_new| to create a new ADAO case, and you
95 .. image:: images/adao_viewer.png
99 **The EFICAS viewer for cases definition in module ADAO**
101 Then fill in the variables to build the ADAO case by using the experimental set
102 up described above. All the technical information given above will be directly
103 inserted in the ADAO case definition, by using the *String* type for all the
104 variables. When the case definition is ready, save it to a "*JDC (\*.comm)*"
105 native file somewhere in your path. Remember that other files will be also
106 created near this first one, so it is better to make a specific directory for
107 your case, and to save the file inside. The name of the file will appear in the
108 "*Object browser*" window, under the "*ADAO*" menu. The final case definition
111 .. _adao_jdcexample01:
112 .. image:: images/adao_jdcexample01.png
116 **Definition of the experimental set up chosen for the ADAO case**
118 To go further, we need now to generate the YACS scheme from the ADAO case
119 definition. In order to do that, right click on the name of the file case in the
120 "*Object browser*" window, and choose the "*Export to YACS*" sub-menu (or the
121 "*Export to YACS*" button |eficas_yacs|) as below:
123 .. _adao_exporttoyacs00:
124 .. image:: images/adao_exporttoyacs.png
128 **"Export to YACS" sub-menu to generate the YACS scheme from the ADAO case**
130 This command will generate the YACS scheme, activate YACS module in SALOME, and
131 open the new scheme in the GUI of the YACS module [#]_. After reordering the
132 nodes by using the "*arrange local node*" sub-menu of the YACS graphical view of
133 the scheme, you get the following representation of the generated ADAO scheme:
135 .. _yacs_generatedscheme:
136 .. image:: images/yacs_generatedscheme.png
140 **YACS generated scheme from the ADAO case**
142 After that point, all the modifications, executions and post-processing of the
143 data assimilation scheme will be done in YACS. In order to check the result in a
144 simple way, we create here a new YACS node by using the "*in-line script node*"
145 sub-menu of the YACS graphical view, and we name it "*PostProcessing*".
147 This script will retrieve the data assimilation analysis from the
148 "*algoResults*" output port of the computation bloc (which gives access to a
149 SALOME Python Object), and will print it on the standard output.
151 To obtain this, the in-line script node need to have an input port of type
152 "*pyobj*" named "*results*" for example, that have to be linked graphically to
153 the "*algoResults*" output port of the computation bloc. Then the code to fill
154 in the script node is::
156 Xa = results.ADD.get("Analysis")[-1]
159 print "Analysis =",Xa
162 The augmented YACS scheme can be saved (overwriting the generated scheme if the
163 simple "*Save*" command or button are used, or with a new name). Ideally, the
164 implementation of such post-processing procedure can be done in YACS to test,
165 and then entirely saved in one script that can be integrated in the ADAO case by
166 using the keyword "*UserPostAnalysis*".
168 Then, classically in YACS, it have to be prepared for run, and then executed.
169 After completion, the printing on standard output is available in the "*YACS
170 Container Log*", obtained through the right click menu of the "*proc*" window in
171 the YACS scheme as shown below:
173 .. _yacs_containerlog:
174 .. image:: images/yacs_containerlog.png
178 **YACS menu for Container Log, and dialog window showing the log**
180 We verify that the result is correct by checking that the log dialog window
181 contains the following line::
183 Analysis = [0.5, 0.5, 0.5]
185 as shown in the image above.
187 As a simple extension of this example, one can notice that the same problem
188 solved with a 3DVAR algorithm gives the same result. This algorithm can be
189 chosen at the ADAO case building step, before entering in YACS step. The
190 ADAO 3DVAR case will look completely similar to the BLUE algorithmic case, as
191 shown by the following figure:
193 .. _adao_jdcexample02:
194 .. image:: images/adao_jdcexample02.png
198 **Defining an ADAO 3DVAR case looks completely similar to a BLUE case**
200 There is only one command changing, with "*3DVAR*" value instead of "*Blue*".
202 Building a simple estimation case with external data definition by scripts
203 --------------------------------------------------------------------------
205 It is useful to get parts or all of the data from external definition, using
206 Python script files to provide access to the data. As an example, we build here
207 an ADAO case representing the same experimental set up as in the above example
208 `Building a simple estimation case with explicit data definition`_, but using
209 data form a single one external Python script file.
211 First, we write the following script file, using conventional names for the
212 desired variables. Here, all the input variables are defined in the script, but
213 the user can choose to split the file in several ones, or to mix explicit data
214 definition in the ADAO GUI and implicit data definition by external files. The
215 present script looks like::
219 # Definition of the Background as a vector
220 # ----------------------------------------
221 Background = [0, 0, 0]
223 # Definition of the Observation as a vector
224 # -----------------------------------------
225 Observation = "1 1 1"
227 # Definition of the Background Error covariance as a matrix
228 # ---------------------------------------------------------
229 BackgroundError = numpy.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
231 # Definition of the Observation Error covariance as a matrix
232 # ----------------------------------------------------------
233 ObservationError = numpy.matrix("1 0 0 ; 0 1 0 ; 0 0 1")
235 # Definition of the Observation Operator as a matrix
236 # --------------------------------------------------
237 ObservationOperator = numpy.identity(3)
239 The names of the Python variables above are mandatory, in order to define the
240 right variables, but the Python script can be bigger and define classes,
241 functions, etc. with other names. It shows different ways to define arrays and
242 matrices, using list, string (as in Numpy or Octave), Numpy array type or Numpy
243 matrix type, and Numpy special functions. All of these syntax are valid.
245 After saving this script somewhere in your path (named here "*script.py*" for
246 the example), we use the GUI to build the ADAO case. The procedure to fill in
247 the case is similar except that, instead of selecting the "*String*" option for
248 the "*FROM*" keyword, we select the "*Script*" one. This leads to a
249 "*SCRIPT_DATA/SCRIPT_FILE*" entry in the tree, allowing to choose a file as:
251 .. _adao_scriptentry01:
252 .. image:: images/adao_scriptentry01.png
256 **Defining an input value using an external script file**
258 Other steps and results are exactly the same as in the `Building a simple
259 estimation case with explicit data definition`_ previous example.
261 In fact, this script methodology allows to retrieve data from in-line or previous
262 calculations, from static files, from database or from stream, all of them
263 outside of SALOME. It allows also to modify easily some input data, for example
264 for debug purpose or for repetitive execution process, and it is the most
265 versatile method in order to parametrize the input data. **But be careful,
266 script methodology is not a "safe" procedure, in the sense that erroneous
267 data, or errors in calculations, can be directly injected into the YACS scheme
270 Adding parameters to control the data assimilation algorithm
271 ------------------------------------------------------------
273 One can add some optional parameters to control the data assimilation algorithm
274 calculation. This is done by using the "*AlgorithmParameters*" keyword in the
275 definition of the ADAO case, which is an keyword of the ASSIMILATION_STUDY. This
276 keyword requires a Python dictionary, containing some key/value pairs. The list
277 of possible optional parameters are given in the subsection
278 :ref:`subsection_algo_options`.
280 If no bounds at all are required on the control variables, then one can choose
281 the "BFGS" or "CG" minimisation algorithm for the 3DVAR algorithm. For
282 constrained optimization, the minimizer "LBFGSB" is often more robust, but the
283 "TNC" is sometimes more performant.
285 This dictionary has to be defined, for example, in an external Python script
286 file, using the mandatory variable name "*AlgorithmParameters*" for the
287 dictionary. All the keys inside the dictionary are optional, they all have
288 default values, and can exist without being used. For example::
290 AlgorithmParameters = {
291 "Minimizer" : "CG", # Possible choice : "LBFGSB", "TNC", "CG", "BFGS"
292 "MaximumNumberOfSteps" : 10,
295 Then the script can be added to the ADAO case, in a file entry describing the
296 "*AlgorithmParameters*" keyword, as follows:
298 .. _adao_scriptentry02:
299 .. image:: images/adao_scriptentry02.png
303 **Adding parameters to control the algorithm**
305 Other steps and results are exactly the same as in the `Building a simple
306 estimation case with explicit data definition`_ previous example. The dictionary
307 can also be directly given in the input field associated with the keyword.
309 Building a complex case with external data definition by scripts
310 ----------------------------------------------------------------
312 This more complex and complete example has to been considered as a framework for
313 user inputs, that need to be tailored for each real application. Nevertheless,
314 the file skeletons are sufficiently general to have been used for various
315 applications in neutronic, fluid mechanics... Here, we will not focus on the
316 results, but more on the user control of inputs and outputs in an ADAO case. As
317 previously, all the numerical values of this example are arbitrary.
319 The objective is to set up the input and output definitions of a physical case
320 by external python scripts, using a general non-linear operator, adding control
321 on parameters and so on... The complete framework scripts can be found in the
322 ADAO skeletons examples directory under the name
323 "*External_data_definition_by_scripts*".
328 We continue to operate in a 3-dimensional space, in order to restrict
329 the size of numerical object shown in the scripts, but the problem is
330 not dependant of the dimension.
332 We choose a twin experiment context, using a known true state
333 :math:`\mathbf{x}^t` of arbitrary values:
337 The background state :math:`\mathbf{x}^b`, which represent some *a priori*
338 knowledge of the true state, is build as a normal random perturbation of 20% the
339 true state :math:`\mathbf{x}^t` for each component, which is:
341 ``Xb = Xt + normal(0, 20%*Xt)``
343 To describe the background error covariances matrix :math:`\mathbf{B}`, we make
344 as previously the hypothesis of uncorrelated errors (that is, a diagonal matrix,
345 of size 3x3 because :math:`\mathbf{x}^b` is of lenght 3) and to have the same
346 variance of 0.1 for all variables. We get:
348 ``B = 0.1 * diagonal( length(Xb) )``
350 We suppose that there exist an observation operator :math:`\mathbf{H}`, which
351 can be non linear. In real calibration procedure or inverse problems, the
352 physical simulation codes are embedded in the observation operator. We need also
353 to know its gradient with respect to each calibrated variable, which is a rarely
354 known information with industrial codes. But we will see later how to obtain an
355 approximated gradient in this case.
357 Being in twin experiments, the observation :math:`\mathbf{y}^o` and its error
358 covariances matrix :math:`\mathbf{R}` are generated by using the true state
359 :math:`\mathbf{x}^t` and the observation operator :math:`\mathbf{H}`:
363 and, with an arbitrary standard deviation of 1% on each error component:
365 ``R = 0.0001 * diagonal( lenght(Yo) )``
367 All the required data assimilation informations are then defined.
369 Skeletons of the scripts describing the setup
370 +++++++++++++++++++++++++++++++++++++++++++++
372 We give here the essential parts of each script used afterwards to build the ADAO
373 case. Remember that using these scripts in real Python files requires to
374 correctly define the path to imported modules or codes (even if the module is in
375 the same directory that the importing Python file ; we indicate the path
376 adjustment using the mention ``"# INSERT PHYSICAL SCRIPT PATH"``), the encoding
377 if necessary, etc. The indicated file names for the following scripts are
378 arbitrary. Examples of complete file scripts are available in the ADAO examples
381 We first define the true state :math:`\mathbf{x}^t` and some convenient matrix
382 building function, in a Python script file named
383 ``Physical_data_and_covariance_matrices.py``::
389 Arbitrary values and names, as a tuple of two series of same length
391 return (numpy.array([1, 2, 3]), ['Para1', 'Para2', 'Para3'])
393 def Simple_Matrix( size, diagonal=None ):
395 Diagonal matrix, with either 1 or a given vector on the diagonal
397 if diagonal is not None:
398 S = numpy.diag( diagonal )
400 S = numpy.matrix(numpy.identity(int(size)))
403 We can then define the background state :math:`\mathbf{x}^b` as a random
404 perturbation of the true state, adding at the end of the script the definition
405 of a *required ADAO variable* in order to export the defined value. It is done
406 in a Python script file named ``Script_Background_xb.py``::
408 from Physical_data_and_covariance_matrices import True_state
411 xt, names = True_state()
413 Standard_deviation = 0.2*xt # 20% for each variable
415 xb = xt + abs(numpy.random.normal(0.,Standard_deviation,size=(len(xt),)))
417 # Creating the required ADAO variable
418 # -----------------------------------
419 Background = list(xb)
421 In the same way, we define the background error covariance matrix
422 :math:`\mathbf{B}` as a diagonal matrix of the same diagonal length as the
423 background of the true state, using the convenient function already defined. It
424 is done in a Python script file named ``Script_BackgroundError_B.py``::
426 from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
428 xt, names = True_state()
430 B = 0.1 * Simple_Matrix( size = len(xt) )
432 # Creating the required ADAO variable
433 # -----------------------------------
436 To continue, we need the observation operator :math:`\mathbf{H}` as a function
437 of the state. It is here defined in an external file named
438 ``"Physical_simulation_functions.py"``, which should contain functions
439 conveniently named here ``"FunctionH"`` and ``"AdjointH"``. These functions are
440 user ones, representing as programming functions the :math:`\mathbf{H}` operator
441 and its adjoint. We suppose these functions are given by the user. A simple
442 skeleton is given here for convenience::
445 """ Direct non-linear simulation operator """
447 # --------------------------------------> EXAMPLE TO BE REMOVED
448 if type(XX) is type(numpy.matrix([])): # EXAMPLE TO BE REMOVED
449 HX = XX.A1.tolist() # EXAMPLE TO BE REMOVED
450 elif type(XX) is type(numpy.array([])): # EXAMPLE TO BE REMOVED
451 HX = numpy.matrix(XX).A1.tolist() # EXAMPLE TO BE REMOVED
452 else: # EXAMPLE TO BE REMOVED
453 HX = XX # EXAMPLE TO BE REMOVED
454 # --------------------------------------> EXAMPLE TO BE REMOVED
456 return numpy.array( HX )
458 We does not need the operators ``"TangentH"`` and ``"AdjointH"`` because they
459 will be approximated using ADAO capabilities.
461 We insist on the fact that these non-linear operator ``"FunctionH"``, tangent
462 operator ``"TangentH"`` and adjoint operator ``"AdjointH"`` come from the
463 physical knowledge, include the reference physical simulation code and its
464 eventual adjoint, and have to be carefully set up by the data assimilation user.
465 The errors in or missuses of the operators can not be detected or corrected by
466 the data assimilation framework alone.
468 In this twin experiments framework, the observation :math:`\mathbf{y}^o` and its
469 error covariances matrix :math:`\mathbf{R}` can be generated. It is done in two
470 Python script files, the first one being named ``Script_Observation_yo.py``::
472 from Physical_data_and_covariance_matrices import True_state
473 from Physical_simulation_functions import FunctionH
475 xt, noms = True_state()
479 # Creating the required ADAO variable
480 # -----------------------------------
481 Observation = list(yo)
483 and the second one named ``Script_ObservationError_R.py``::
485 from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
486 from Physical_simulation_functions import FunctionH
488 xt, names = True_state()
492 R = 0.0001 * Simple_Matrix( size = len(yo) )
494 # Creating the required ADAO variable
495 # -----------------------------------
498 As in previous examples, it can be useful to define some parameters for the data
499 assimilation algorithm. For example, if we use the standard 3DVAR algorithm, the
500 following parameters can be defined in a Python script file named
501 ``Script_AlgorithmParameters.py``::
503 # Creating the required ADAO variable
504 # -----------------------------------
505 AlgorithmParameters = {
506 "Minimizer" : "TNC", # Possible : "LBFGSB", "TNC", "CG", "BFGS"
507 "MaximumNumberOfSteps" : 15, # Number of global iterative steps
509 [ None, None ], # Bound on the first parameter
510 [ 0., 4. ], # Bound on the second parameter
511 [ 0., None ], # Bound on the third parameter
515 Finally, it is common to post-process the results, retrieving them after the
516 data assimilation phase in order to analyze, print or show them. It requires to
517 use a intermediary Python script file in order to extract these results. The
518 following example Python script file named ``Script_UserPostAnalysis.py``,
519 illustrates the fact::
521 from Physical_data_and_covariance_matrices import True_state
524 xt, names = True_state()
525 xa = ADD.get("Analysis")[-1]
526 x_series = ADD.get("CurrentState")[:]
527 J = ADD.get("CostFunctionJ")[:]
529 # Verifying the results by printing
530 # ---------------------------------
533 print "xa = %s"%numpy.array(xa)
535 for i in range( len(x_series) ):
536 print "Step %2i : J = %.5e et X = %s"%(i, J[i], x_series[i])
539 At the end, we get a description of the whole case setup through a set of files
542 #. ``Physical_data_and_covariance_matrices.py``
543 #. ``Physical_simulation_functions.py``
544 #. ``Script_AlgorithmParameters.py``
545 #. ``Script_BackgroundError_B.py``
546 #. ``Script_Background_xb.py``
547 #. ``Script_ObservationError_R.py``
548 #. ``Script_Observation_yo.py``
549 #. ``Script_UserPostAnalysis.py``
551 We insist here that all these scripts are written by the user and can not be
552 automatically tested. So the user is required to verify the scripts (and in
553 particular their input/output) in order to limit the difficulty of debug. We
554 recall: **script methodology is not a "safe" procedure, in the sense that
555 erroneous data, or errors in calculations, can be directly injected into the
556 YACS scheme execution.**
558 Building the case with external data definition by scripts
559 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
561 All these scripts can then be used to define the ADAO case with external data
562 definition by Python script files. It is entirely similar to the method
563 described in the `Building a simple estimation case with external data
564 definition by scripts`_ previous section. For each variable to be defined, we
565 select the "*Script*" option of the "*FROM*" keyword, which leads to a
566 "*SCRIPT_DATA/SCRIPT_FILE*" entry in the tree. For the "*ObservationOperator*"
567 keyword, we choose the "*ScriptWithOneFunction*" form and keep the default
568 differential increment.
570 The other steps to build the ADAO case are exactly the same as in the `Building
571 a simple estimation case with explicit data definition`_ previous section.
573 Using the simple linear operator :math:`\mathbf{H}` from the Python script file
574 ``Physical_simulation_functions.py`` in the ADAO examples standard directory,
575 the results will look like::
578 xa = [ 1.000014 2.000458 3.000390]
580 Step 0 : J = 1.81750e+03 et X = [1.014011, 2.459175, 3.390462]
581 Step 1 : J = 1.81750e+03 et X = [1.014011, 2.459175, 3.390462]
582 Step 2 : J = 1.79734e+01 et X = [1.010771, 2.040342, 2.961378]
583 Step 3 : J = 1.79734e+01 et X = [1.010771, 2.040342, 2.961378]
584 Step 4 : J = 1.81909e+00 et X = [1.000826, 2.000352, 3.000487]
585 Step 5 : J = 1.81909e+00 et X = [1.000826, 2.000352, 3.000487]
586 Step 6 : J = 1.81641e+00 et X = [1.000247, 2.000651, 3.000156]
587 Step 7 : J = 1.81641e+00 et X = [1.000247, 2.000651, 3.000156]
588 Step 8 : J = 1.81569e+00 et X = [1.000015, 2.000432, 3.000364]
589 Step 9 : J = 1.81569e+00 et X = [1.000015, 2.000432, 3.000364]
590 Step 10 : J = 1.81568e+00 et X = [1.000013, 2.000458, 3.000390]
593 The state at the first step is the randomly generated background state
594 :math:`\mathbf{x}^b`. After completion, these printing on standard output is
595 available in the "*YACS Container Log*", obtained through the right click menu
596 of the "*proc*" window in the YACS scheme.
598 .. [#] For more information on YACS, see the the *YACS module User's Guide* available in the main "*Help*" menu of SALOME GUI.