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 explicitly
22 all the required input data through the GUI. The second one shows, on the same
23 case, how to define input data using external sources through scripts. We
24 describe here always Python scripts because they can be directly inserted in
25 YACS script nodes, but external files can use other langages.
27 The mathematical notations used afterward are explained in the section
28 :ref:`section_theory`.
30 Building an estimation case with explicit data definition
31 ---------------------------------------------------------
33 This simple example is a demonstration one, and describes how to set a BLUE
34 estimation framework in order to get *ponderated (or fully weighted) least
35 square estimated state* of a system from an observation of the state and from an
36 *a priori* knowledge (or background) of this state. In other words, we look for
37 the weighted middle between the observation and the background vectors. All the
38 numerical values of this example are arbitrary.
43 We choose to operate in a 3-dimensional space. 3D is chosen in order to restrict
44 the size of numerical object to explicitly enter by the user, but the problem is
45 not dependant of the dimension and can be set in dimension 10, 100, 1000... The
46 observation :math:`\mathbf{y}^o` is of value 1 in each direction, so::
50 The background state :math:`\mathbf{x}^b`, which represent some *a priori*
51 knowledge or a mathematical regularization, is of value of 0 in each direction,
56 Data assimilation requires information on errors covariances :math:`\mathbf{R}`
57 and :math:`\mathbf{B}`, respectively for observation and background variables.
58 We choose here to have uncorrelated errors (that is, diagonal matrices) and to
59 have the same variance of 1 for all variables (that is, identity matrices). We
62 B = R = [1 0 0 ; 0 1 0 ; 0 0 1]
64 Last, we need an observation operator :math:`\mathbf{H}` to convert the
65 background value in the space of observation values. Here, because the space
66 dimensions are the same, we can choose the identity as the observation
69 H = [1 0 0 ; 0 1 0 ; 0 0 1]
71 With such choices, the "Best Linear Unbiased Estimator" (BLUE) will be the
72 average vector between :math:`\mathbf{y}^o` and :math:`\mathbf{x}^b`, named the
73 *analysis*, denoted by :math:`\mathbf{x}^a`, and its value is::
77 As an extension of this example, one can change the variances represented by
78 :math:`\mathbf{B}` or :math:`\mathbf{R}` independently, and the analysis
79 :math:`\mathbf{x}^a` will move to :math:`\mathbf{y}^o` or to
80 :math:`\mathbf{x}^b`, in inverse proportion of the variances in
81 :math:`\mathbf{B}` and :math:`\mathbf{R}`. As an other extension, it is also
82 equivalent to search for the analysis thought a BLUE algorithm or a 3DVAR one.
84 Using the graphical interface (GUI) to build the ADAO case
85 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 First, you have to activate the ADAO module by choosing the appropriate module
88 button or menu of SALOME, and you will see:
91 .. image:: images/adao_activate.png
95 **Activating the module ADAO in SALOME**
97 Choose the "*New*" button in this window. You will directly get the EFICAS
98 interface for variables definition, along with the SALOME "*Object browser*".
99 You can then click on the "*New*" button |eficas_new| to create a new ADAO case,
103 .. image:: images/adao_viewer.png
107 **The EFICAS editor for cases definition in module ADAO**
109 Then, fill in the variables to build the ADAO case by using the experimental set
110 up described above. All the technical information given above will be directly
111 inserted in the ADAO case definition, by using the *String* type for all the
112 variables. When the case definition is ready, save it to a "*JDC (\*.comm)*"
113 native file somewhere in your path. Remember that other files will be also
114 created near this first one, so it is better to make a specific directory for
115 your case, and to save the file inside. The name of the file will appear in the
116 "*Object browser*" window, under the "*ADAO*" menu. The final case definition
119 .. _adao_jdcexample01:
120 .. image:: images/adao_jdcexample01.png
124 **Definition of the experimental setup chosen for the ADAO case**
126 To go further, we need now to generate the YACS scheme from the ADAO case
127 definition. In order to do that, right click on the name of the file case in the
128 "*Object browser*" window, and choose the "*Export to YACS*" sub-menu (or the
129 "*Export to YACS*" button |eficas_yacs|) as below:
131 .. _adao_exporttoyacs00:
132 .. image:: images/adao_exporttoyacs.png
136 **"Export to YACS" sub-menu to generate the YACS scheme from the ADAO case**
138 This command will generate the YACS scheme, activate YACS module in SALOME, and
139 open the new scheme in the GUI of the YACS module [#]_. After eventually
140 reordering the nodes by using the "*arrange local nodes*" sub-menu of the YACS
141 graphical view of the scheme, you get the following representation of the
142 generated ADAO scheme:
144 .. _yacs_generatedscheme:
145 .. image:: images/yacs_generatedscheme.png
149 **YACS generated scheme from the ADAO case**
151 After that point, all the modifications, executions and post-processing of the
152 data assimilation scheme will be done in the YACS module. In order to check the
153 result in a simple way, we create here a new YACS node by using the "*in-line
154 script node*" sub-menu of the YACS graphical view, and we name it
157 This script node will retrieve the data assimilation analysis from the
158 "*algoResults*" output port of the computation bloc (which gives access to a
159 SALOME Python Object), and will print it on the standard output.
161 To obtain this, the in-line script node need to have an input port of type
162 "*pyobj*", named "*results*" for example, that have to be linked graphically to
163 the "*algoResults*" output port of the computation bloc. Then, the code to fill
164 in the script node is::
166 Xa = results.ADD.get("Analysis")[-1]
169 print "Analysis =",Xa
172 The augmented YACS scheme can be saved (overwriting the generated scheme if the
173 "*Save*" command or button are used, or with a new name through the "*Save as*"
174 command). Ideally, the implementation of such post-processing procedure can be
175 done in YACS to test, and then entirely saved in one Python script that can be
176 integrated in the ADAO case by using the keyword "*UserPostAnalysis*".
178 Then, classically in YACS, the scheme have to be compiled for run, and then
179 executed. After completion, the printing on standard output is available in the
180 "*YACS Container Log*", obtained through the right click menu of the "*proc*"
181 window in the YACS scheme as shown below:
183 .. _yacs_containerlog:
184 .. image:: images/yacs_containerlog.png
188 **YACS menu for Container Log, and dialog window showing the log**
190 We verify that the result is correct by checking that the log dialog window
191 contains the following line::
193 Analysis = [0.5, 0.5, 0.5]
195 as shown in the image above.
197 As a simple extension of this example, one can notice that the same problem
198 solved with a 3DVAR algorithm gives the same result. This algorithm can be
199 chosen at the ADAO case building step, before entering in YACS step. The
200 ADAO 3DVAR case will look completely similar to the BLUE algorithmic case, as
201 shown by the following figure:
203 .. _adao_jdcexample02:
204 .. image:: images/adao_jdcexample02.png
208 **Defining an ADAO 3DVAR case looks completely similar to a BLUE case**
210 There is only one command changing, with "*3DVAR*" value in the "*Algorithm*"
211 field instead of "*Blue*".
213 Building a estimation case with external data definition by scripts
214 -------------------------------------------------------------------
216 It is useful to get parts or all of the data from external definition, using
217 Python script files to provide access to the data. As an example, we build here
218 an ADAO case representing the same experimental setup as in the above example
219 `Building an estimation case with explicit data definition`_, but using data
220 from a single one external Python script file.
222 First, we write the following script file, using conventional names for the
223 required variables. Here, all the input variables are defined in the same
224 script, but the user can choose to split the file in several ones, or to mix
225 explicit data definition in the ADAO GUI and implicit data definition by
226 external files. The present script file looks like::
230 # Definition of the Background as a vector
231 # ----------------------------------------
232 Background = [0, 0, 0]
234 # Definition of the Observation as a vector
235 # -----------------------------------------
236 Observation = "1 1 1"
238 # Definition of the Background Error covariance as a matrix
239 # ---------------------------------------------------------
240 BackgroundError = numpy.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
242 # Definition of the Observation Error covariance as a matrix
243 # ----------------------------------------------------------
244 ObservationError = numpy.matrix("1 0 0 ; 0 1 0 ; 0 0 1")
246 # Definition of the Observation Operator as a matrix
247 # --------------------------------------------------
248 ObservationOperator = numpy.identity(3)
250 The names of the Python variables above are mandatory, in order to define the
251 right case variables, but the Python script can be bigger and define classes,
252 functions, file or database access, etc. with other names. Moreover, the above
253 script shows different ways to define arrays and matrices, using list, string
254 (as in Numpy or Octave), Numpy array type or Numpy matrix type, and Numpy
255 special functions. All of these syntax are valid.
257 After saving this script in a file (named here "*script.py*" for the example)
258 somewhere in your path, we use the graphical interface (GUI) to build the ADAO
259 case. The procedure to fill in the case is similar to the previous example
260 except that, instead of selecting the "*String*" option for the "*FROM*" keyword
261 of each variable, we select the "*Script*" one. This leads to a
262 "*SCRIPT_DATA/SCRIPT_FILE*" entry in the graphical tree, allowing to choose a
265 .. _adao_scriptentry01:
266 .. image:: images/adao_scriptentry01.png
270 **Defining an input value using an external script file**
272 Other steps and results are exactly the same as in the `Building an estimation
273 case with explicit data definition`_ previous example.
275 In fact, this script methodology is the easiest way to retrieve data from
276 in-line or previous calculations, from static files, from database or from
277 stream, all of them inside or outside of SALOME. It allows also to modify easily
278 some input data, for example for debug purpose or for repetitive execution
279 process, and it is the most versatile method in order to parametrize the input
280 data. **But be careful, script methodology is not a "safe" procedure, in the
281 sense that erroneous data, or errors in calculations, can be directly injected
282 into the YACS scheme execution.**
284 Adding parameters to control the data assimilation algorithm
285 ------------------------------------------------------------
287 One can add some optional parameters to control the data assimilation algorithm
288 calculation. This is done by using the "*AlgorithmParameters*" keyword in the
289 definition of the ADAO case, which is an keyword of the "*ASSIMILATION_STUDY*"
290 general command. This keyword requires a Python dictionary, containing some
291 key/value pairs. The list of possible optional parameters are given in the
292 subsection :ref:`section_reference`.
294 This dictionary has to be defined, for example, in an external Python script
295 file, using the mandatory variable name "*AlgorithmParameters*" for the
296 dictionary. All the keys inside the dictionary are optional, they all have
297 default values, and can exist without being used. For example::
299 AlgorithmParameters = {
300 "Minimizer" : "CG", # Possible choice : "LBFGSB", "TNC", "CG", "BFGS"
301 "MaximumNumberOfSteps" : 10,
304 If no bounds at all are required on the control variables, then one can choose
305 the "*BFGS*" or "*CG*" minimisation algorithm for all the variational data
306 assimilation or optimization algorithms. For constrained optimization, the
307 minimizer "*LBFGSB*" is often more robust, but the "*TNC*" is sometimes more
310 Then the script can be added to the ADAO case, in a file entry describing the
311 "*AlgorithmParameters*" keyword, as follows:
313 .. _adao_scriptentry02:
314 .. image:: images/adao_scriptentry02.png
318 **Adding parameters to control the algorithm and the outputs**
320 Other steps and results are exactly the same as in the `Building an estimation
321 case with explicit data definition`_ previous example. The dictionary can also
322 be directly given in the input field of string type associated for the keyword.
324 Building a complex case with external data definition by scripts
325 ----------------------------------------------------------------
327 This more complex and complete example has to been considered as a framework for
328 user inputs treatment, that need to be tailored for each real application.
329 Nevertheless, the file skeletons are sufficiently general to have been used for
330 various applications in neutronic, fluid mechanics... Here, we will not focus on
331 the results, but more on the user control of inputs and outputs in an ADAO case.
332 As previously, all the numerical values of this example are arbitrary.
334 The objective is to setup the input and output definitions of a physical
335 estimation case by external python scripts, using a general non-linear operator,
336 adding control on parameters and so on... The complete framework scripts can be
337 found in the ADAO skeletons examples directory under the name
338 "*External_data_definition_by_scripts*".
343 We continue to operate in a 3-dimensional space, in order to restrict
344 the size of numerical object shown in the scripts, but the problem is
345 not dependant of the dimension.
347 We choose a twin experiment context, using a known true state
348 :math:`\mathbf{x}^t` but of arbitrary value::
352 The background state :math:`\mathbf{x}^b`, which represent some *a priori*
353 knowledge of the true state, is build as a normal random perturbation of 20% of
354 the true state :math:`\mathbf{x}^t` for each component, which is::
356 Xb = Xt + normal(0, 20%*Xt)
358 To describe the background error covariances matrix :math:`\mathbf{B}`, we make
359 as previously the hypothesis of uncorrelated errors (that is, a diagonal matrix,
360 of size 3x3 because :math:`\mathbf{x}^b` is of lenght 3) and to have the same
361 variance of 0.1 for all variables. We get::
363 B = 0.1 * diagonal( length(Xb) )
365 We suppose that there exist an observation operator :math:`\mathbf{H}`, which
366 can be non linear. In real calibration procedure or inverse problems, the
367 physical simulation codes are embedded in the observation operator. We need also
368 to know its gradient with respect to each calibrated variable, which is a rarely
369 known information with industrial codes. But we will see later how to obtain an
370 approximated gradient in this case.
372 Being in twin experiments, the observation :math:`\mathbf{y}^o` and its error
373 covariances matrix :math:`\mathbf{R}` are generated by using the true state
374 :math:`\mathbf{x}^t` and the observation operator :math:`\mathbf{H}`::
378 and, with an arbitrary standard deviation of 1% on each error component::
380 R = 0.0001 * diagonal( lenght(Yo) )
382 All the informations required for estimation by data assimilation are then
385 Skeletons of the scripts describing the setup
386 +++++++++++++++++++++++++++++++++++++++++++++
388 We give here the essential parts of each script used afterwards to build the
389 ADAO case. Remember that using these scripts in real Python files requires to
390 correctly define the path to imported modules or codes (even if the module is in
391 the same directory that the importing Python file. We indicate the path
392 adjustment using the mention ``"# INSERT PHYSICAL SCRIPT PATH"``), the encoding
393 if necessary, etc. The indicated file names for the following scripts are
394 arbitrary. Examples of complete file scripts are available in the ADAO examples
397 We first define the true state :math:`\mathbf{x}^t` and some convenient matrix
398 building function, in a Python script file named
399 ``Physical_data_and_covariance_matrices.py``::
405 Arbitrary values and names, as a tuple of two series of same length
407 return (numpy.array([1, 2, 3]), ['Para1', 'Para2', 'Para3'])
409 def Simple_Matrix( size, diagonal=None ):
411 Diagonal matrix, with either 1 or a given vector on the diagonal
413 if diagonal is not None:
414 S = numpy.diag( diagonal )
416 S = numpy.matrix(numpy.identity(int(size)))
419 We can then define the background state :math:`\mathbf{x}^b` as a random
420 perturbation of the true state, adding a *required ADAO variable* at the end of
421 the script the definition, in order to export the defined value. It is done in a
422 Python script file named ``Script_Background_xb.py``::
424 from Physical_data_and_covariance_matrices import True_state
427 xt, names = True_state()
429 Standard_deviation = 0.2*xt # 20% for each variable
431 xb = xt + abs(numpy.random.normal(0.,Standard_deviation,size=(len(xt),)))
433 # Creating the required ADAO variable
434 # -----------------------------------
435 Background = list(xb)
437 In the same way, we define the background error covariance matrix
438 :math:`\mathbf{B}` as a diagonal matrix, of the same diagonal length as the
439 background of the true state, using the convenient function already defined. It
440 is done in a Python script file named ``Script_BackgroundError_B.py``::
442 from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
444 xt, names = True_state()
446 B = 0.1 * Simple_Matrix( size = len(xt) )
448 # Creating the required ADAO variable
449 # -----------------------------------
452 To continue, we need the observation operator :math:`\mathbf{H}` as a function
453 of the state. It is here defined in an external file named
454 ``"Physical_simulation_functions.py"``, which should contain one function
455 conveniently named here ``"DirectOperator"``. This function is user one,
456 representing as programming function the :math:`\mathbf{H}` operator. We suppose
457 this function is then given by the user. A simple skeleton is given here for
460 def DirectOperator( XX ):
461 """ Direct non-linear simulation operator """
463 # --------------------------------------> EXAMPLE TO BE REMOVED
464 if type(XX) is type(numpy.matrix([])): # EXAMPLE TO BE REMOVED
465 HX = XX.A1.tolist() # EXAMPLE TO BE REMOVED
466 elif type(XX) is type(numpy.array([])): # EXAMPLE TO BE REMOVED
467 HX = numpy.matrix(XX).A1.tolist() # EXAMPLE TO BE REMOVED
468 else: # EXAMPLE TO BE REMOVED
469 HX = XX # EXAMPLE TO BE REMOVED
470 # --------------------------------------> EXAMPLE TO BE REMOVED
472 return numpy.array( HX )
474 We does not need the linear compagnion operators ``"TangentOperator"`` and
475 ``"AdjointOperator"`` because they will be approximated using ADAO capabilities.
477 We insist on the fact that these non-linear operator ``"DirectOperator"``,
478 tangent operator ``"TangentOperator"`` and adjoint operator
479 ``"AdjointOperator"`` come from the physical knowledge, include the reference
480 physical simulation code, and have to be carefully setup by the data
481 assimilation or optimization user. The simulation errors or missuses of the
482 operators can not be detected or corrected by the data assimilation and
483 optimization ADAO framework alone.
485 In this twin experiments framework, the observation :math:`\mathbf{y}^o` and its
486 error covariances matrix :math:`\mathbf{R}` can be generated. It is done in two
487 Python script files, the first one being named ``Script_Observation_yo.py``::
489 from Physical_data_and_covariance_matrices import True_state
490 from Physical_simulation_functions import DirectOperator
492 xt, noms = True_state()
494 yo = DirectOperator( xt )
496 # Creating the required ADAO variable
497 # -----------------------------------
498 Observation = list(yo)
500 and the second one named ``Script_ObservationError_R.py``::
502 from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
503 from Physical_simulation_functions import DirectOperator
505 xt, names = True_state()
507 yo = DirectOperator( xt )
509 R = 0.0001 * Simple_Matrix( size = len(yo) )
511 # Creating the required ADAO variable
512 # -----------------------------------
515 As in previous examples, it can be useful to define some parameters for the data
516 assimilation algorithm. For example, if we use the standard "*3DVAR*" algorithm,
517 the following parameters can be defined in a Python script file named
518 ``Script_AlgorithmParameters.py``::
520 # Creating the required ADAO variable
521 # -----------------------------------
522 AlgorithmParameters = {
523 "Minimizer" : "TNC", # Possible : "LBFGSB", "TNC", "CG", "BFGS"
524 "MaximumNumberOfSteps" : 15, # Number of global iterative steps
526 [ None, None ], # Bound on the first parameter
527 [ 0., 4. ], # Bound on the second parameter
528 [ 0., None ], # Bound on the third parameter
532 Finally, it is common to post-process the results, retrieving them after the
533 data assimilation phase in order to analyze, print or show them. It requires to
534 use a intermediary Python script file in order to extract these results at the
535 end of the adata assimilation or optimization process. The following example
536 Python script file, named ``Script_UserPostAnalysis.py``, illustrates the fact::
538 from Physical_data_and_covariance_matrices import True_state
541 xt, names = True_state()
542 xa = ADD.get("Analysis")[-1]
543 x_series = ADD.get("CurrentState")[:]
544 J = ADD.get("CostFunctionJ")[:]
546 # Verifying the results by printing
547 # ---------------------------------
550 print "xa = %s"%numpy.array(xa)
552 for i in range( len(x_series) ):
553 print "Step %2i : J = %.5e and X = %s"%(i, J[i], x_series[i])
556 At the end, we get a description of the whole case setup through a set of files
559 #. ``Physical_data_and_covariance_matrices.py``
560 #. ``Physical_simulation_functions.py``
561 #. ``Script_AlgorithmParameters.py``
562 #. ``Script_BackgroundError_B.py``
563 #. ``Script_Background_xb.py``
564 #. ``Script_ObservationError_R.py``
565 #. ``Script_Observation_yo.py``
566 #. ``Script_UserPostAnalysis.py``
568 We insist here that all these scripts are written by the user and can not be
569 automatically tested by ADAO. So the user is required to verify the scripts (and
570 in particular their input/output) in order to limit the difficulty of debug. We
571 recall: **script methodology is not a "safe" procedure, in the sense that
572 erroneous data, or errors in calculations, can be directly injected into the
573 YACS scheme execution.**
575 Building the case with external data definition by scripts
576 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
578 All these scripts can then be used to define the ADAO case with external data
579 definition by Python script files. It is entirely similar to the method
580 described in the `Building a simple estimation case with external data
581 definition by scripts`_ previous section. For each variable to be defined, we
582 select the "*Script*" option of the "*FROM*" keyword, which leads to a
583 "*SCRIPT_DATA/SCRIPT_FILE*" entry in the tree. For the "*ObservationOperator*"
584 keyword, we choose the "*ScriptWithOneFunction*" form and keep the default
585 differential increment.
587 The other steps to build the ADAO case are exactly the same as in the `Building
588 an estimation case with explicit data definition`_ previous section.
590 Using the simple linear operator :math:`\mathbf{H}` from the Python script file
591 ``Physical_simulation_functions.py`` in the ADAO examples standard directory,
592 the results will look like::
595 xa = [ 1.000014 2.000458 3.000390]
597 Step 0 : J = 1.81750e+03 and X = [1.014011, 2.459175, 3.390462]
598 Step 1 : J = 1.81750e+03 and X = [1.014011, 2.459175, 3.390462]
599 Step 2 : J = 1.79734e+01 and X = [1.010771, 2.040342, 2.961378]
600 Step 3 : J = 1.79734e+01 and X = [1.010771, 2.040342, 2.961378]
601 Step 4 : J = 1.81909e+00 and X = [1.000826, 2.000352, 3.000487]
602 Step 5 : J = 1.81909e+00 and X = [1.000826, 2.000352, 3.000487]
603 Step 6 : J = 1.81641e+00 and X = [1.000247, 2.000651, 3.000156]
604 Step 7 : J = 1.81641e+00 and X = [1.000247, 2.000651, 3.000156]
605 Step 8 : J = 1.81569e+00 and X = [1.000015, 2.000432, 3.000364]
606 Step 9 : J = 1.81569e+00 and X = [1.000015, 2.000432, 3.000364]
607 Step 10 : J = 1.81568e+00 and X = [1.000013, 2.000458, 3.000390]
610 The state at the first step is the randomly generated background state
611 :math:`\mathbf{x}^b`. During calculation, these printings on standard output are
612 available in the "*YACS Container Log*" window, obtained through the right click
613 menu of the "*proc*" window in the YACS executed scheme.
615 .. [#] For more information on YACS, see the *YACS module User's Guide* available in the main "*Help*" menu of SALOME platform.