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
26 ================================================================================
27 **[DocU]** Tutorials on using the ADAO module
28 ================================================================================
30 .. |eficas_new| image:: images/eficas_new.png
33 .. |eficas_save| image:: images/eficas_save.png
36 .. |eficas_saveas| image:: images/eficas_saveas.png
39 .. |eficas_yacs| image:: images/eficas_yacs.png
43 This section presents some examples on using the ADAO module in SALOME. The
44 first one shows how to build a simple data assimilation case defining explicitly
45 all the required input data through the GUI. The second one shows, on the same
46 case, how to define input data using external sources through scripts. We
47 describe here always Python scripts because they can be directly inserted in
48 YACS script nodes, but external files can use other languages.
50 The mathematical notations used afterward are explained in the section
51 :ref:`section_theory`.
53 Building an estimation case with explicit data definition
54 ---------------------------------------------------------
56 This simple example is a demonstration one, and describes how to set a BLUE
57 estimation framework in order to get the *fully weighted least square estimated
58 state* of a system from an observation of the state and from an *a priori*
59 knowledge (or background) of this state. In other words, we look for the
60 weighted middle between the observation and the background vectors. All the
61 numerical values of this example are arbitrary.
66 We choose to operate in a 3-dimensional space. 3D is chosen in order to restrict
67 the size of numerical object to explicitly enter by the user, but the problem is
68 not dependent of the dimension and can be set in dimension 10, 100, 1000... The
69 observation :math:`\mathbf{y}^o` is of value 1 in each direction, so::
73 The background state :math:`\mathbf{x}^b`, which represent some *a priori*
74 knowledge or a mathematical regularization, is of value of 0 in each direction,
79 Data assimilation requires information on errors covariances :math:`\mathbf{R}`
80 and :math:`\mathbf{B}`, respectively for observation and background variables.
81 We choose here to have uncorrelated errors (that is, diagonal matrices) and to
82 have the same variance of 1 for all variables (that is, identity matrices). We
85 B = R = [1 0 0 ; 0 1 0 ; 0 0 1]
87 Last, we need an observation operator :math:`\mathbf{H}` to convert the
88 background value in the space of observation values. Here, because the space
89 dimensions are the same, we can choose the identity as the observation
92 H = [1 0 0 ; 0 1 0 ; 0 0 1]
94 With such choices, the "Best Linear Unbiased Estimator" (BLUE) will be the
95 average vector between :math:`\mathbf{y}^o` and :math:`\mathbf{x}^b`, named the
96 *analysis*, denoted by :math:`\mathbf{x}^a`, and its value is::
100 As an extension of this example, one can change the variances represented by
101 :math:`\mathbf{B}` or :math:`\mathbf{R}` independently, and the analysis
102 :math:`\mathbf{x}^a` will move to :math:`\mathbf{y}^o` or to
103 :math:`\mathbf{x}^b`, in inverse proportion of the variances in
104 :math:`\mathbf{B}` and :math:`\mathbf{R}`. As an other extension, it is also
105 equivalent to search for the analysis thought a BLUE algorithm or a 3DVAR one.
107 Using the graphical interface (GUI) to build the ADAO case
108 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
110 First, you have to activate the ADAO module by choosing the appropriate module
111 button or menu of SALOME, and you will see:
114 .. image:: images/adao_activate.png
118 **Activating the module ADAO in SALOME**
120 Choose the "*New*" button in this window. You will directly get the embedded
121 case editor interface for variables definition, along with the SALOME "*Object
122 browser*". You can then click on the "*New*" button |eficas_new| to create a new
123 ADAO case, and you will see:
126 .. image:: images/adao_viewer.png
130 **The embedded editor for cases definition in module ADAO**
132 Then, fill in the variables to build the ADAO case by using the experimental set
133 up described above. All the technical information given above will be directly
134 inserted in the ADAO case definition, by using the *String* type for all the
135 variables. When the case definition is ready, save it to a "*JDC (\*.comm)*"
136 native file somewhere in your path. Remember that other files will be also
137 created near this first one, so it is better to make a specific directory for
138 your case, and to save the file inside. The name of the file will appear in the
139 "*Object browser*" window, under the "*ADAO*" menu. The final case definition
142 .. _adao_jdcexample01:
143 .. image:: images/adao_jdcexample01.png
147 **Definition of the experimental setup chosen for the ADAO case**
149 To go further, we need now to generate the YACS scheme from the ADAO case
150 definition. In order to do that, right click on the name of the file case in the
151 "*Object browser*" window, and choose the "*Export to YACS*" sub-menu (or the
152 "*Export to YACS*" button |eficas_yacs|) as below:
154 .. _adao_exporttoyacs00:
155 .. image:: images/adao_exporttoyacs.png
159 **"Export to YACS" sub-menu to generate the YACS scheme from the ADAO case**
161 This command will generate the YACS scheme, activate YACS module in SALOME, and
162 open the new scheme in the GUI of the YACS module [#]_. After eventually
163 reordering the nodes by using the "*arrange local nodes*" sub-menu of the YACS
164 graphical view of the scheme, you get the following representation of the
165 generated ADAO scheme:
167 .. _yacs_generatedscheme:
168 .. image:: images/yacs_generatedscheme.png
172 **YACS generated scheme from the ADAO case**
174 After that point, all the modifications, executions and post-processing of the
175 data assimilation scheme will be done in the YACS module. In order to check the
176 result in a simple way, we use the "*UserPostAnalysis*" node (or we create here
177 a new YACS node by using the "*in-line script node*" sub-menu of the YACS
180 This script node will retrieve the data assimilation analysis from the
181 "*algoResults*" output port of the computation bloc (which gives access to a
182 SALOME Python Object), and will print it on the standard output.
184 To obtain this, the in-line script node need to have an input port of type
185 "*pyobj*", named "*Study*" for example, that have to be linked graphically to
186 the "*algoResults*" output port of the computation bloc. Then, the code to fill
187 in the script node is::
189 Xa = Study.getResults().get("Analysis")[-1]
192 print "Analysis =",Xa
195 The (initial or augmented) YACS scheme can be saved (overwriting the generated
196 scheme if the "*Save*" command or button are used, or with a new name through
197 the "*Save as*" command). Ideally, the implementation of such post-processing
198 procedure can be done in YACS to test, and then entirely saved in one Python
199 script that can be integrated in the ADAO case by using the keyword
200 "*UserPostAnalysis*".
202 Then, classically in YACS, the scheme have to be compiled for run, and then
203 executed. After completion, the printing on standard output is available in the
204 "*YACS Container Log*", obtained through the right click menu of the "*proc*"
205 window in the YACS scheme as shown below:
207 .. _yacs_containerlog:
208 .. image:: images/yacs_containerlog.png
212 **YACS menu for Container Log, and dialog window showing the log**
214 We verify that the result is correct by checking that the log dialog window
215 contains the following line::
217 Analysis = [0.5, 0.5, 0.5]
219 as shown in the image above.
221 As a simple extension of this example, one can notice that the same problem
222 solved with a 3DVAR algorithm gives the same result. This algorithm can be
223 chosen at the ADAO case building step, before entering in YACS step. The
224 ADAO 3DVAR case will look completely similar to the BLUE algorithmic case, as
225 shown by the following figure:
227 .. _adao_jdcexample02:
228 .. image:: images/adao_jdcexample02.png
232 **Defining an ADAO 3DVAR case looks completely similar to a BLUE case**
234 There is only one command changing, with "*3DVAR*" value in the "*Algorithm*"
235 field instead of "*Blue*".
237 Building an estimation case with external data definition by scripts
238 --------------------------------------------------------------------
240 It is useful to get parts or all of the data from external definition, using
241 Python script files to provide access to the data. As an example, we build here
242 an ADAO case representing the same experimental setup as in the above example
243 `Building an estimation case with explicit data definition`_, but using data
244 from a single one external Python script file.
246 First, we write the following script file, using conventional names for the
247 required variables. Here, all the input variables are defined in the same
248 script, but the user can choose to split the file in several ones, or to mix
249 explicit data definition in the ADAO GUI and implicit data definition by
250 external files. The present script file looks like::
254 # Definition of the Background as a vector
255 # ----------------------------------------
256 Background = [0, 0, 0]
258 # Definition of the Observation as a vector
259 # -----------------------------------------
260 Observation = "1 1 1"
262 # Definition of the Background Error covariance as a matrix
263 # ---------------------------------------------------------
264 BackgroundError = numpy.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
266 # Definition of the Observation Error covariance as a matrix
267 # ----------------------------------------------------------
268 ObservationError = numpy.matrix("1 0 0 ; 0 1 0 ; 0 0 1")
270 # Definition of the Observation Operator as a matrix
271 # --------------------------------------------------
272 ObservationOperator = numpy.identity(3)
274 The names of the Python variables above are mandatory, in order to define the
275 right case variables, but the Python script can be bigger and define classes,
276 functions, file or database access, etc. with other names. Moreover, the above
277 script shows different ways to define arrays and matrices, using list, string
278 (as in Numpy or Octave), Numpy array type or Numpy matrix type, and Numpy
279 special functions. All of these syntax are valid.
281 After saving this script in a file (named here "*script.py*" for the example)
282 somewhere in your path, we use the graphical interface (GUI) to build the ADAO
283 case. The procedure to fill in the case is similar to the previous example
284 except that, instead of selecting the "*String*" option for the "*FROM*" keyword
285 of each variable, we select the "*Script*" one. This leads to a
286 "*SCRIPT_DATA/SCRIPT_FILE*" entry in the graphical tree, allowing to choose a
289 .. _adao_scriptentry01:
290 .. image:: images/adao_scriptentry01.png
294 **Defining an input value using an external script file**
296 Other steps and results are exactly the same as in the `Building an estimation
297 case with explicit data definition`_ previous example.
299 In fact, this script methodology is the easiest way to retrieve data from
300 in-line or previous calculations, from static files, from database or from
301 stream, all of them inside or outside of SALOME. It allows also to modify easily
302 some input data, for example for debug purpose or for repetitive execution
303 process, and it is the most versatile method in order to parametrize the input
304 data. **But be careful, script methodology is not a "safe" procedure, in the
305 sense that erroneous data, or errors in calculations, can be directly injected
306 into the YACS scheme execution. The user have to carefully verify the content of
309 Adding parameters to control the data assimilation algorithm
310 ------------------------------------------------------------
312 One can add some optional parameters to control the data assimilation algorithm
313 calculation. This is done by using optional parameters in the
314 "*AlgorithmParameters*" command of the ADAO case definition, which is a keyword
315 of the "*ASSIMILATION_STUDY*" general command. This keyword requires an explicit
316 definition of the values from default ones, or from a Python dictionary,
317 containing some key/value pairs. The list of possible optional parameters are
318 given in the section :ref:`section_reference` and its subsections. The
319 recommandation is to use the explicit definition of values from the default list
320 of optionnal parameters, as here with the "*MaximumNumberOfSteps*":
322 .. _adao_scriptentry02:
323 .. image:: images/adao_scriptentry02.png
327 **Adding parameters to control the algorithm and the outputs**
329 This dictionary can be defined, for example, in an external Python script
330 file, using the mandatory variable name "*AlgorithmParameters*" for the
331 dictionary. All the keys inside the dictionary are optional, they all have
332 default values, and can exist without being used. For example::
334 AlgorithmParameters = {
335 "Minimizer" : "LBFGSB", # Recommended
336 "MaximumNumberOfSteps" : 10,
339 If no bounds at all are required on the control variables, then one can choose
340 the "*BFGS*" or "*CG*" minimization algorithm for all the variational data
341 assimilation or optimization algorithms. For constrained optimization, the
342 minimizer "*LBFGSB*" is often more robust, but the "*TNC*" is sometimes more
343 effective. In a general way, the "*LBFGSB*" algorithm choice is recommended.
344 Then the script can be added to the ADAO case, in a file entry describing the
345 "*Parameters*" keyword.
347 Other steps and results are exactly the same as in the `Building an estimation
348 case with explicit data definition`_ previous example. The dictionary can also
349 be directly given in the input field of string type associated for the keyword.
351 Building a complex case with external data definition by scripts
352 ----------------------------------------------------------------
354 This more complex and complete example has to been considered as a framework for
355 user inputs treatment, that need to be tailored for each real application.
356 Nevertheless, the file skeletons are sufficiently general to have been used for
357 various applications in neutronic, fluid mechanics... Here, we will not focus on
358 the results, but more on the user control of inputs and outputs in an ADAO case.
359 As previously, all the numerical values of this example are arbitrary.
361 The objective is to setup the input and output definitions of a physical
362 estimation case by external python scripts, using a general non-linear operator,
363 adding control on parameters and so on... The complete framework scripts can be
364 found in the ADAO skeletons examples directory under the name
365 "*External_data_definition_by_scripts*".
370 We continue to operate in a 3-dimensional space, in order to restrict
371 the size of numerical object shown in the scripts, but the problem is
372 not dependent of the dimension.
374 We choose a twin experiment context, using a known true state
375 :math:`\mathbf{x}^t` but of arbitrary value::
379 The background state :math:`\mathbf{x}^b`, which represent some *a priori*
380 knowledge of the true state, is build as a normal random perturbation of 20% of
381 the true state :math:`\mathbf{x}^t` for each component, which is::
383 Xb = Xt + normal(0, 20%*Xt)
385 To describe the background error covariances matrix :math:`\mathbf{B}`, we make
386 as previously the hypothesis of uncorrelated errors (that is, a diagonal matrix,
387 of size 3x3 because :math:`\mathbf{x}^b` is of lenght 3) and to have the same
388 variance of 0.1 for all variables. We get::
390 B = 0.1 * diagonal( length(Xb) )
392 We suppose that there exist an observation operator :math:`\mathbf{H}`, which
393 can be non linear. In real calibration procedure or inverse problems, the
394 physical simulation codes are embedded in the observation operator. We need also
395 to know its gradient with respect to each calibrated variable, which is a rarely
396 known information with industrial codes. But we will see later how to obtain an
397 approximated gradient in this case.
399 Being in twin experiments, the observation :math:`\mathbf{y}^o` and its error
400 covariances matrix :math:`\mathbf{R}` are generated by using the true state
401 :math:`\mathbf{x}^t` and the observation operator :math:`\mathbf{H}`::
405 and, with an arbitrary standard deviation of 1% on each error component::
407 R = 0.0001 * diagonal( length(Yo) )
409 All the information required for estimation by data assimilation are then
412 Skeletons of the scripts describing the setup
413 +++++++++++++++++++++++++++++++++++++++++++++
415 We give here the essential parts of each script used afterward to build the
416 ADAO case. Remember that using these scripts in real Python files requires to
417 correctly define the path to imported modules or codes (even if the module is in
418 the same directory that the importing Python file. We indicate the path
419 adjustment using the mention ``"# INSERT PHYSICAL SCRIPT PATH"``), the encoding
420 if necessary, etc. The indicated file names for the following scripts are
421 arbitrary. Examples of complete file scripts are available in the ADAO examples
424 We first define the true state :math:`\mathbf{x}^t` and some convenient matrix
425 building function, in a Python script file named
426 ``Physical_data_and_covariance_matrices.py``::
432 Arbitrary values and names, as a tuple of two series of same length
434 return (numpy.array([1, 2, 3]), ['Para1', 'Para2', 'Para3'])
436 def Simple_Matrix( size, diagonal=None ):
438 Diagonal matrix, with either 1 or a given vector on the diagonal
440 if diagonal is not None:
441 S = numpy.diag( diagonal )
443 S = numpy.matrix(numpy.identity(int(size)))
446 We can then define the background state :math:`\mathbf{x}^b` as a random
447 perturbation of the true state, adding a *required ADAO variable* at the end of
448 the script the definition, in order to export the defined value. It is done in a
449 Python script file named ``Script_Background_xb.py``::
451 from Physical_data_and_covariance_matrices import True_state
454 xt, names = True_state()
456 Standard_deviation = 0.2*xt # 20% for each variable
458 xb = xt + abs(numpy.random.normal(0.,Standard_deviation,size=(len(xt),)))
460 # Creating the required ADAO variable
461 # -----------------------------------
462 Background = list(xb)
464 In the same way, we define the background error covariance matrix
465 :math:`\mathbf{B}` as a diagonal matrix, of the same diagonal length as the
466 background of the true state, using the convenient function already defined. It
467 is done in a Python script file named ``Script_BackgroundError_B.py``::
469 from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
471 xt, names = True_state()
473 B = 0.1 * Simple_Matrix( size = len(xt) )
475 # Creating the required ADAO variable
476 # -----------------------------------
479 To continue, we need the observation operator :math:`\mathbf{H}` as a function
480 of the state. It is here defined in an external file named
481 ``"Physical_simulation_functions.py"``, which should contain one function
482 conveniently named here ``"DirectOperator"``. This function is user one,
483 representing as programming function the :math:`\mathbf{H}` operator. We suppose
484 this function is then given by the user. A simple skeleton is given here for
487 def DirectOperator( XX ):
488 """ Direct non-linear simulation operator """
490 # --------------------------------------> EXAMPLE TO BE REMOVED
491 if type(XX) is type(numpy.matrix([])): # EXAMPLE TO BE REMOVED
492 HX = XX.A1.tolist() # EXAMPLE TO BE REMOVED
493 elif type(XX) is type(numpy.array([])): # EXAMPLE TO BE REMOVED
494 HX = numpy.matrix(XX).A1.tolist() # EXAMPLE TO BE REMOVED
495 else: # EXAMPLE TO BE REMOVED
496 HX = XX # EXAMPLE TO BE REMOVED
497 # --------------------------------------> EXAMPLE TO BE REMOVED
499 return numpy.array( HX )
501 We does not need the linear companion operators ``"TangentOperator"`` and
502 ``"AdjointOperator"`` because they will be approximated using ADAO capabilities.
504 We insist on the fact that these non-linear operator ``"DirectOperator"``,
505 tangent operator ``"TangentOperator"`` and adjoint operator
506 ``"AdjointOperator"`` come from the physical knowledge, include the reference
507 physical simulation code, and have to be carefully setup by the data
508 assimilation or optimization user. The simulation errors or missuses of the
509 operators can not be detected or corrected by the data assimilation and
510 optimization ADAO framework alone.
512 In this twin experiments framework, the observation :math:`\mathbf{y}^o` and its
513 error covariances matrix :math:`\mathbf{R}` can be generated. It is done in two
514 Python script files, the first one being named ``Script_Observation_yo.py``::
516 from Physical_data_and_covariance_matrices import True_state
517 from Physical_simulation_functions import DirectOperator
519 xt, noms = True_state()
521 yo = DirectOperator( xt )
523 # Creating the required ADAO variable
524 # -----------------------------------
525 Observation = list(yo)
527 and the second one named ``Script_ObservationError_R.py``::
529 from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
530 from Physical_simulation_functions import DirectOperator
532 xt, names = True_state()
534 yo = DirectOperator( xt )
536 R = 0.0001 * Simple_Matrix( size = len(yo) )
538 # Creating the required ADAO variable
539 # -----------------------------------
542 As in previous examples, it can be useful to define some parameters for the data
543 assimilation algorithm. For example, if we use the standard "*3DVAR*" algorithm,
544 the following parameters can be defined in a Python script file named
545 ``Script_AlgorithmParameters.py``::
547 # Creating the required ADAO variable
548 # -----------------------------------
549 AlgorithmParameters = {
550 "Minimizer" : "LBFGSB", # Recommended
551 "MaximumNumberOfSteps" : 15, # Number of global iterative steps
553 [ None, None ], # Bound on the first parameter
554 [ 0., 4. ], # Bound on the second parameter
555 [ 0., None ], # Bound on the third parameter
559 Finally, it is common to post-process the results, retrieving them after the
560 data assimilation phase in order to analyze, print or show them. It requires to
561 use a intermediary Python script file in order to extract these results at the
562 end of the a data assimilation or optimization process. The following example
563 Python script file, named ``Script_UserPostAnalysis.py``, illustrates the fact::
565 from Physical_data_and_covariance_matrices import True_state
568 xt, names = True_state()
569 xa = ADD.get("Analysis")[-1]
570 x_series = ADD.get("CurrentState")[:]
571 J = ADD.get("CostFunctionJ")[:]
573 # Verifying the results by printing
574 # ---------------------------------
577 print "xa = %s"%numpy.array(xa)
579 for i in range( len(x_series) ):
580 print "Step %2i : J = %.5e and X = %s"%(i, J[i], x_series[i])
583 At the end, we get a description of the whole case setup through a set of files
586 #. ``Physical_data_and_covariance_matrices.py``
587 #. ``Physical_simulation_functions.py``
588 #. ``Script_AlgorithmParameters.py``
589 #. ``Script_BackgroundError_B.py``
590 #. ``Script_Background_xb.py``
591 #. ``Script_ObservationError_R.py``
592 #. ``Script_Observation_yo.py``
593 #. ``Script_UserPostAnalysis.py``
595 We insist here that all these scripts are written by the user and can not be
596 automatically tested by ADAO. So the user is required to verify the scripts (and
597 in particular their input/output) in order to limit the difficulty of debug. We
598 recall: **script methodology is not a "safe" procedure, in the sense that
599 erroneous data, or errors in calculations, can be directly injected into the
600 YACS scheme execution.**
602 Building the case with external data definition by scripts
603 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
605 All these scripts can then be used to define the ADAO case with external data
606 definition by Python script files. It is entirely similar to the method
607 described in the `Building an estimation case with external data definition by
608 scripts`_ previous section. For each variable to be defined, we select the
609 "*Script*" option of the "*FROM*" keyword, which leads to a
610 "*SCRIPT_DATA/SCRIPT_FILE*" entry in the tree. For the "*ObservationOperator*"
611 keyword, we choose the "*ScriptWithOneFunction*" form and keep the default
612 differential increment.
614 The other steps to build the ADAO case are exactly the same as in the `Building
615 an estimation case with explicit data definition`_ previous section.
617 Using the simple linear operator :math:`\mathbf{H}` from the Python script file
618 ``Physical_simulation_functions.py`` in the ADAO examples standard directory,
619 the results will look like::
622 xa = [ 1.000014 2.000458 3.000390]
624 Step 0 : J = 1.81750e+03 and X = [1.014011, 2.459175, 3.390462]
625 Step 1 : J = 1.81750e+03 and X = [1.014011, 2.459175, 3.390462]
626 Step 2 : J = 1.79734e+01 and X = [1.010771, 2.040342, 2.961378]
627 Step 3 : J = 1.79734e+01 and X = [1.010771, 2.040342, 2.961378]
628 Step 4 : J = 1.81909e+00 and X = [1.000826, 2.000352, 3.000487]
629 Step 5 : J = 1.81909e+00 and X = [1.000826, 2.000352, 3.000487]
630 Step 6 : J = 1.81641e+00 and X = [1.000247, 2.000651, 3.000156]
631 Step 7 : J = 1.81641e+00 and X = [1.000247, 2.000651, 3.000156]
632 Step 8 : J = 1.81569e+00 and X = [1.000015, 2.000432, 3.000364]
633 Step 9 : J = 1.81569e+00 and X = [1.000015, 2.000432, 3.000364]
634 Step 10 : J = 1.81568e+00 and X = [1.000013, 2.000458, 3.000390]
637 The state at the first step is the randomly generated background state
638 :math:`\mathbf{x}^b`. During calculation, these printings on standard output are
639 available in the "*YACS Container Log*" window, obtained through the right click
640 menu of the "*proc*" window in the YACS executed scheme.
642 .. [#] For more information on YACS, see the *YACS module* and its integrated help available from the main menu *Help* of the SALOME platform.