Salome HOME
Documentation de toutes les options du 3DVAR
[modules/adao.git] / doc / examples.rst
1 .. _section_examples:
2
3 ================================================================================
4 Examples on using the ADAO module
5 ================================================================================
6
7 .. |eficas_new| image:: images/eficas_new.png
8    :align: middle
9    :scale: 50%
10 .. |eficas_save| image:: images/eficas_save.png
11    :align: middle
12    :scale: 50%
13 .. |eficas_yacs| image:: images/eficas_yacs.png
14    :align: middle
15    :scale: 50%
16
17 This section presents some examples on using the ADAO module in SALOME. The
18 first one shows how to build a simple data assimilation case defining
19 explicitly all the required data through the GUI. The second one shows, on the
20 same case, how to define data using external sources through scripts.
21
22 Building a simple estimation case with explicit data definition
23 ---------------------------------------------------------------
24
25 This simple example is a demonstration one, and describes how to set a BLUE
26 estimation framework in order to get *weighted least square estimated state* of
27 a system from an observation of the state and from an *a priori* knowledge (or
28 background) of this state. In other words, we look for the weighted middle
29 between the observation and the background vectors. All the numerical values of
30 this example are arbitrary.
31
32 Experimental set up
33 +++++++++++++++++++
34
35 We choose to operate in a 3-dimensional space. 3D is chosen in order to restrict
36 the size of numerical object to explicitly enter by the user, but the problem is
37 not dependant of the dimension and can be set in dimension 1000... The
38 observation :math:`\mathbf{y}^o` is of value 1 in each direction, so:
39
40     ``Yo = [1 1 1]``
41
42 The background state :math:`\mathbf{x}^b`, which represent some *a priori*
43 knowledge or a regularization, is of value of 0 in each direction, which is:
44
45     ``Xb = [0 0 0]``
46
47 Data assimilation requires information on errors covariances :math:`\mathbf{R}`
48 and :math:`\mathbf{B}` respectively for observation and background variables. We
49 choose here to have uncorrelated errors (that is, diagonal matrices) and to have
50 the same variance of 1 for all variables (that is, identity matrices). We get:
51
52     ``B = R = [1 0 0 ; 0 1 0 ; 0 0 1]``
53
54 Last, we need an observation operator :math:`\mathbf{H}` to convert the
55 background value in the space of observation value. Here, because the space
56 dimensions are the same, we can choose the identity  as the observation
57 operator:
58
59     ``H = [1 0 0 ; 0 1 0 ; 0 0 1]``
60
61 With such choices, the Best Linear Unbiased Estimator (BLUE) will be the average
62 vector between :math:`\mathbf{y}^o` and :math:`\mathbf{x}^b`, named the
63 *analysis* and denoted by :math:`\mathbf{x}^a`:
64
65     ``Xa = [0.5 0.5 0.5]``
66
67 As en extension of this example, one can change the variances for
68 :math:`\mathbf{B}` or :math:`\mathbf{R}` independently, and the analysis will
69 move to :math:`\mathbf{y}^o` or to :math:`\mathbf{x}^b` in inverse proportion of
70 the variances in :math:`\mathbf{B}` and :math:`\mathbf{R}`. It is also
71 equivalent to search for the analysis thought a BLUE algorithm or a 3DVAR one.
72
73 Using the GUI to build the ADAO case
74 ++++++++++++++++++++++++++++++++++++
75
76 First, you have to activate the ADAO module by choosing the appropriate module
77 button or menu of SALOME, and you will see:
78
79   .. _adao_activate2:
80   .. image:: images/adao_activate.png
81     :align: center
82     :width: 100%
83   .. centered::
84     **Activating the module ADAO in SALOME**
85
86 Choose the "*New*" button in this window. You will directly get the EFICAS
87 interface for variables definition, along with the "*Object browser*". You can
88 then click on the "*New*" button |eficas_new| to create a new ADAO case, and you
89 will see:
90
91   .. _adao_viewer:
92   .. image:: images/adao_viewer.png
93     :align: center
94     :width: 100%
95   .. centered::
96     **The EFICAS viewer for cases definition in module ADAO**
97
98 Then fill in the variables to build the ADAO case by using the experimental set
99 up described above. All the technical information given above will be directly
100 inserted in the ADAO case definition, by using the *String* type for all the
101 variables. When the case definition is ready, save it to a "*JDC (\*.comm)*"
102 native file somewhere in your path. Remember that other files will be also
103 created near this first one, so it is better to make a specific directory for
104 your case, and to save the file inside. The name of the file will appear in the
105 "*Object browser*" window, under the "*ADAO*" menu. The final case definition
106 looks like this:
107
108   .. _adao_jdcexample01:
109   .. image:: images/adao_jdcexample01.png
110     :align: center
111     :width: 100%
112   .. centered::
113     **Definition of the experimental set up chosen for the ADAO case**
114
115 To go further, we need now to generate the YACS scheme from the ADAO case
116 definition. In order to do that, right click on the name of the file case in the
117 "*Object browser*" window, and choose the "*Export to YACS*" sub-menu (or the
118 "*Export to YACS*" button |eficas_yacs|) as below:
119
120   .. _adao_exporttoyacs:
121   .. image:: images/adao_exporttoyacs.png
122     :align: center
123     :scale: 75%
124   .. centered::
125     **"Export to YACS" sub-menu to generate the YACS scheme from the ADAO case**
126
127 This command will generate the YACS scheme, activate YACS module in SALOME, and
128 open the new scheme in the GUI of the YACS module [#]_. After reordering the
129 nodes by using the "*arrange local node*" sub-menu of the YACS graphical view of
130 the scheme, you get the following representation of the generated ADAO scheme:
131
132   .. _yacs_generatedscheme:
133   .. image:: images/yacs_generatedscheme.png
134     :align: center
135     :width: 100%
136   .. centered::
137     **YACS generated scheme from the ADAO case**
138
139 After that point, all the modifications, executions and post-processing of the
140 data assimilation scheme will be done in YACS. In order to check the result in a
141 simple way, we create here a new YACS node by using the "*in-line script node*"
142 sub-menu of the YACS graphical view, and we name it "*PostProcessing*".
143
144 This script will retrieve the data assimilation analysis from the
145 "*algoResults*" output port of the computation bloc (which gives access to a
146 SALOME Python Object), and will print it on the standard output. 
147
148 To obtain this, the in-line script node need to have an input port of type
149 "*pyobj*" named "*results*" for example, that have to be linked graphically to
150 the "*algoResults*" output port of the computation bloc. Then the code to fill
151 in the script node is::
152
153     Xa = results.ADD.get("Analysis").valueserie(-1)
154
155     print
156     print "Analysis =",Xa
157     print
158
159 The augmented YACS scheme can be saved (overwriting the generated scheme if the
160 simple "*Save*" command or button are used, or with a new name). Then,
161 classically in YACS, it have to be prepared for run, and then executed. After
162 completion, the printing on standard output is available in the "*YACS Container
163 Log*", obtained through the right click menu of the "*proc*" window in the YACS
164 scheme as shown below:
165
166   .. _yacs_containerlog:
167   .. image:: images/yacs_containerlog.png
168     :align: center
169     :width: 100%
170   .. centered::
171     **YACS menu for Container Log, and dialog window showing the log**
172
173 We verify that the result is correct by checking that the log dialog window
174 contains the following line::
175
176     Analysis = [0.5, 0.5, 0.5]
177
178 as shown in the image above.
179
180 As a simple extension of this example, one can notice that the same problem
181 solved with a 3DVAR algorithm gives the same result. This algorithm can be
182 chosen at the ADAO case building step, before entering in YACS step. The
183 ADAO 3DVAR case will look completely similar to the BLUE algorithmic case, as
184 shown by the following figure:
185
186   .. _adao_jdcexample02:
187   .. image:: images/adao_jdcexample02.png
188     :align: center
189     :width: 100%
190   .. centered::
191     **Defining an ADAO 3DVAR case looks completely similar to a BLUE case**
192
193 There is only one command changing, with "*3DVAR*" value instead of "*Blue*".
194
195 Building a simple estimation case with external data definition by scripts
196 --------------------------------------------------------------------------
197
198 It is useful to get parts or all of the data from external definition, using
199 Python script files to provide access to the data. As an example, we build here
200 an ADAO case representing the same experimental set up as in the above example
201 `Building a simple estimation case with explicit data definition`_, but using
202 data form a single one external Python script file.
203
204 First, we write the following script file, using conventional names for the
205 desired variables. Here, all the input variables are defined in the script, but
206 the user can choose to split the file in several ones, or to mix explicit data
207 definition in the ADAO GUI and implicit data definition by external files. The
208 present script looks like::
209
210     #-*-coding:iso-8859-1-*-
211     #
212     import numpy
213     #
214     # Definition of the Background as a vector
215     # ----------------------------------------
216     Background = [0, 0, 0]
217     #
218     # Definition of the Observation as a vector
219     # -----------------------------------------
220     Observation = "1 1 1"
221     #
222     # Definition of the Background Error covariance as a matrix
223     # ---------------------------------------------------------
224     BackgroundError = numpy.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
225     #
226     # Definition of the Observation Error covariance as a matrix
227     # ----------------------------------------------------------
228     ObservationError = numpy.matrix("1 0 0 ; 0 1 0 ; 0 0 1")
229     #
230     # Definition of the Observation Operator as a matrix
231     # --------------------------------------------------
232     ObservationOperator = numpy.identity(3)
233
234 The names of the Python variables above are mandatory, in order to define the
235 right variables, but the Python script can be bigger and define classes,
236 functions, etc. with other names. It shows different ways to define arrays and
237 matrices, using list, string (as in Numpy or Octave), Numpy array type or Numpy
238 matrix type, and Numpy special functions. All of these syntaxes are valid.
239
240 After saving this script somewhere in your path (named here "*script.py*" for
241 the example), we use the GUI to build the ADAO case. The procedure to fill in
242 the case is similar except that, instead of selecting the "*String*" option for
243 the "*FROM*" keyword, we select the "*Script*" one. This leads to a
244 "*SCRIPT_DATA/SCRIPT_FILE*" entry in the tree, allowing to choose a file as:
245
246   .. _adao_scriptentry01:
247   .. image:: images/adao_scriptentry01.png
248     :align: center
249     :width: 100%
250   .. centered::
251     **Defining an input value using an external script file**
252
253 Other steps and results are exactly the same as in the `Building a simple
254 estimation case with explicit data definition`_ previous example.
255
256 In fact, this script methodology allows to retrieve data from in-line or previous
257 calculations, from static files, from database or from stream, all of them
258 outside of SALOME. It allows also to modify easily some input data, for example
259 for debug purpose or for repetitive execution process, and it is the most
260 versatile method in order to parametrize the input data. **But be careful,
261 script methodology is not a "safe" procedure, in the sense that erroneous
262 data, or errors in calculations, can be directly injected into the YACS scheme
263 execution.**
264
265 Adding parameters to control the data assimilation algorithm
266 ------------------------------------------------------------
267
268 One can add some optional parameters to control the data assimilation algorithm
269 calculation. This is done by using the "*AlgorithmParameters*" keyword in the
270 definition of the ADAO case, which is an keyword of the ASSIMILATION_STUDY. This
271 keyword requires a Python dictionary, containing some key/value pairs.
272
273 For example, with a 3DVAR algorithm, the possible keys are "*Minimizer*",
274 "*MaximumNumberOfSteps*", "ProjectedGradientTolerance", "GradientNormTolerance"
275 and "*Bounds*":
276
277 #.   The "*Minimizer*" key allows to choose the optimisation minimizer. The
278      default choice is "LBFGSB", and the possible ones are "LBFGSB" (nonlinear
279      constrained minimizer, see [Byrd95] and [Zhu97]), "TNC" (nonlinear
280      constrained minimizer), "CG" (nonlinear unconstrained minimizer), "BFGS"
281      (nonlinear unconstrained minimizer), "NCG" (Newton CG minimizer).
282 #.   The "*MaximumNumberOfSteps*" key indicates the maximum number of iterations
283      allowed for iterative optimisation. The default is 15000, which very
284      similar of no limit on iterations. It is then recommended to adapt this
285      parameter to the needs on real problems.
286 #.   The "ProjectedGradientTolerance" key indicates a limit value, leading to
287      stop successfully the iterative optimisation process when all the components
288      of the projected gradient are under this limit.
289 #.   The "GradientNormTolerance" key indicates a limit value, leading to stop
290      successfully the iterative optimisation process when the norm of the
291      gradient is under this limit.
292 #.   The "*Bounds*" key allows to define upper and lower bounds for every
293      control variable being optimized. Bounds can be given by a list of list of
294      pairs of lower/upper bounds for each variable, with possibly ``None`` every
295      time there is no bound. The bounds can always be specified, but they are
296      taken into account only by the constrained minimizers.
297
298 If no bounds at all are required on the control variables, then one can choose
299 the "BFGS" or "CG" minimisation algorithm for the 3DVAR algorithm. For
300 constrained optimisation, the minimizer "LBFGSB" is often more robust, but the
301 "TNC" is always more performant.
302
303 This dictionary has to be defined, for example, in an external Python script
304 file, using the mandatory variable name "*AlgorithmParameters*" for the
305 dictionary. All the keys inside the dictionary are optional, they all have
306 default values, and can exist without being used. For example::
307
308     #-*-coding:iso-8859-1-*-
309     #
310     AlgorithmParameters = {
311         "Minimizer" : "CG", # Possible choice : "LBFGSB", "TNC", "CG", "BFGS"
312         "MaximumNumberOfSteps" : 10,
313         }
314
315 Then the script can be added to the ADAO case, in a file entry describing the
316 "*AlgorithmParameters*" keyword, as follows:
317
318   .. _adao_scriptentry02:
319   .. image:: images/adao_scriptentry02.png
320     :align: center
321     :width: 100%
322   .. centered::
323     **Adding parameters to control the algorithm**
324
325 Other steps and results are exactly the same as in the `Building a simple
326 estimation case with explicit data definition`_ previous example. The dictionary
327 can also be directly given in the input field associated with the keyword.
328
329 Building a complex case with external data definition by scripts
330 ----------------------------------------------------------------
331
332 This more complex and complete example has to been considered as a framework for
333 user inputs, that need to be tailored for each real application. Nevertheless,
334 the file skeletons are sufficiently general to have been used for various
335 applications in neutronic, fluid mechanics... Here, we will not focus on the
336 results, but more on the user control of inputs and outputs in an ADAO case. As
337 previously, all the numerical values of this example are arbitrary.
338
339 The objective is to set up the input and output definitions of a physical case
340 by external python scripts, using a general non-linear operator, adding control
341 on parameters and so on... The complete framework scripts can be found in the
342 ADAO examples standard directory.
343
344 Experimental set up
345 +++++++++++++++++++
346
347 We continue to operate in a 3-dimensional space, in order to restrict
348 the size of numerical object shown in the scripts, but the problem is
349 not dependant of the dimension. 
350
351 We choose a twin experiment context, using a known true state
352 :math:`\mathbf{x}^t` of arbitrary values:
353
354     ``Xt = [1 2 3]``
355
356 The background state :math:`\mathbf{x}^b`, which represent some *a priori*
357 knowledge of the true state, is build as a normal random perturbation of 20% the
358 true state :math:`\mathbf{x}^t` for each component, which is:
359
360     ``Xb = Xt + normal(0,20%*Xt)``
361
362 To describe the background error covariances matrix :math:`\mathbf{B}`, we make
363 as previously the hypothesis of uncorrelated errors (that is, a diagonal matrix,
364 of size 3x3 because :math:`\mathbf{x}^b` is of lenght 3) and to have the same
365 variance of 0.1 for all variables. We get:
366
367     ``B = 0.1 * diagonal( lenght(Xb) )``
368
369 We suppose that there exist an observation operator :math:`\mathbf{H}`, which
370 can be non linear. In real calibration procedure or inverse problems, the
371 physical simulation codes are embedded in the observation operator. We need also
372 to know its gradient with respect to each calibrated variable, which is a rarely
373 known information with industrial codes. But we will see later how to obtain an
374 approximated gradient in this case.
375
376 Being in twin experiments, the observation :math:`\mathbf{y}^o` and its error
377 covariances matrix :math:`\mathbf{R}` are generated by using the true state
378 :math:`\mathbf{x}^t` and the observation operator :math:`\mathbf{H}`:
379
380     ``Yo = H( Xt )``
381
382 and, with an arbitrary standard deviation of 1% on each error component:
383
384     ``R = 0.0001 * diagonal( lenght(Yo) )``
385
386 All the required data assimilation informations are then defined.
387
388 Skeletons of the scripts describing the setup
389 +++++++++++++++++++++++++++++++++++++++++++++
390
391 We give here the essential parts of each script used afterwards to build the ADAO
392 case. Remember that using these scripts in real Python files requires to
393 correctly define the path to imported modules or codes (even if the module is in
394 the same directory that the importing Python file ; we indicate the path
395 adjustment using the mention ``"# INSERT PHYSICAL SCRIPT PATH"``), the encoding
396 if necessary, etc. The indicated file names for the following scripts are
397 arbitrary. Examples of complete file scripts are available in the ADAO examples
398 standard directory.
399
400 We first define the true state :math:`\mathbf{x}^t` and some convenient matrix
401 building function, in a Python script file named
402 ``Physical_data_and_covariance_matrices.py``::
403
404     #-*-coding:iso-8859-1-*-
405     #
406     import numpy
407     #
408     def True_state():
409         """
410         Arbitrary values and names, as a tuple of two series of same length
411         """
412         return (numpy.array([1, 2, 3]), ['Para1', 'Para2', 'Para3'])
413     #
414     def Simple_Matrix( size, diagonal=None ):
415         """
416         Diagonal matrix, with either 1 or a given vector on the diagonal
417         """
418         if diagonal is not None:
419             S = numpy.diag( diagonal )
420         else:
421             S = numpy.matrix(numpy.identity(int(size)))
422         return S
423
424 We can then define the background state :math:`\mathbf{x}^b` as a random
425 perturbation of the true state, adding at the end of the script the definition
426 of a *required ADAO variable* in order to export the defined value. It is done
427 in a Python script file named ``Script_Background_xb.py``::
428
429     #-*-coding:iso-8859-1-*-
430     #
431     import sys ; sys.path.insert(0,"# INSERT PHYSICAL SCRIPT PATH")
432     from Physical_data_and_covariance_matrices import True_state
433     import numpy
434     #
435     xt, names = True_state()
436     #
437     Standard_deviation = 0.2*xt # 20% for each variable
438     #
439     xb = xt + abs(numpy.random.normal(0.,Standard_deviation,size=(len(xt),)))
440     #
441     # Creating the required ADAO variable
442     # -----------------------------------
443     Background = list(xb)
444
445 In the same way, we define the background error covariance matrix
446 :math:`\mathbf{B}` as a diagonal matrix of the same diagonal length as the
447 background of the true state, using the convenient function already defined. It
448 is done in a Python script file named ``Script_BackgroundError_B.py``::
449
450     #-*-coding:iso-8859-1-*-
451     #
452     import sys ; sys.path.insert(0,"# INSERT PHYSICAL SCRIPT PATH")
453     from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
454     #
455     xt, names = True_state()
456     #
457     B = 0.1 * Simple_Matrix( size = len(xt) )
458     #
459     # Creating the required ADAO variable
460     # -----------------------------------
461     BackgroundError = B
462
463 To continue, we need the observation operator :math:`\mathbf{H}` as a function
464 of the state. It is here defined in an external file named
465 ``"Physical_simulation_functions.py"``, which should contain functions
466 conveniently named here ``"FunctionH"`` and ``"AdjointH"``. These functions are
467 user ones, representing as programming functions the :math:`\mathbf{H}` operator
468 and its adjoint. We suppose these functions are given by the user (a simple
469 skeleton is given in the Python script file ``Physical_simulation_functions.py``
470 of the ADAO examples standard directory, not reproduced here).
471
472 To operates in ADAO, it is required to define different types of operators: the
473 (potentially non-linear) standard observation operator, named ``"Direct"``, its
474 linearised approximation, named ``"Tangent"``, and the adjoint operator named
475 ``"Adjoint"``. The Python script have to retrieve an input parameter, found
476 under the key "value", in a variable named ``"specificParameters"`` of the
477 SALOME input data and parameters ``"computation"`` dictionary variable. If the
478 operator is already linear, the ``"Direct"`` and ``"Tangent"`` functions are the
479 same, as it is supposed here. The following example Python script file named
480 ``Script_ObservationOperator_H.py``, illustrates the case::
481
482     #-*-coding:iso-8859-1-*-
483     #
484     import sys ; sys.path.insert(0,"# INSERT PHYSICAL SCRIPT PATH")
485     import Physical_simulation_functions
486     import numpy, logging
487     #
488     # -----------------------------------------------------------------------
489     # SALOME input data and parameters: all information are the required input
490     # variable "computation", containing for example:
491     #      {'inputValues': [[[[0.0, 0.0, 0.0]]]],
492     #       'inputVarList': ['adao_default'],
493     #       'outputVarList': ['adao_default'],
494     #       'specificParameters': [{'name': 'method', 'value': 'Direct'}]}
495     # -----------------------------------------------------------------------
496     #
497     # Recovering the type of computation: "Direct", "Tangent" or "Adjoint"
498     # --------------------------------------------------------------------
499     method = ""
500     for param in computation["specificParameters"]:
501         if param["name"] == "method":
502             method = param["value"]
503     logging.info("ComputationFunctionNode: Found method is \'%s\'"%method)
504     #
505     # Loading the H operator functions from external definitions
506     # ----------------------------------------------------------
507     logging.info("ComputationFunctionNode: Loading operator functions")
508     FunctionH = Physical_simulation_functions.FunctionH
509     AdjointH  = Physical_simulation_functions.AdjointH
510     #
511     # Executing the possible computations
512     # -----------------------------------
513     if method == "Direct":
514         logging.info("ComputationFunctionNode: Direct computation")
515         Xcurrent = computation["inputValues"][0][0][0]
516         data = FunctionH(numpy.matrix( Xcurrent ).T)
517     #
518     if method == "Tangent":
519         logging.info("ComputationFunctionNode: Tangent computation")
520         Xcurrent = computation["inputValues"][0][0][0]
521         data = FunctionH(numpy.matrix( Xcurrent ).T)
522     #
523     if method == "Adjoint":
524         logging.info("ComputationFunctionNode: Adjoint computation")
525         Xcurrent = computation["inputValues"][0][0][0]
526         Ycurrent = computation["inputValues"][0][0][1]
527         data = AdjointH((numpy.matrix( Xcurrent ).T, numpy.matrix( Ycurrent ).T))
528     #
529     # Formatting the output
530     # ---------------------
531     logging.info("ComputationFunctionNode: Formatting the output")
532     it = data.flat
533     outputValues = [[[[]]]]
534     for val in it:
535       outputValues[0][0][0].append(val)
536     #
537     # Creating the required ADAO variable
538     # -----------------------------------
539     result = {}
540     result["outputValues"]        = outputValues
541     result["specificOutputInfos"] = []
542     result["returnCode"]          = 0
543     result["errorMessage"]        = ""
544
545 As output, this script has to define a nested list variable, as shown above with
546 the ``"outputValues"`` variable, where the nested levels describe the different
547 variables included in the state, then the different possible states at the same
548 time, then the different time steps. In this case, because there is only one
549 time step and one state, and all the variables are stored together, we only set
550 the most inner level of the lists.
551
552 In this twin experiments framework, the observation :math:`\mathbf{y}^o` and its
553 error covariances matrix :math:`\mathbf{R}` can be generated. It is done in two
554 Python script files, the first one being named ``Script_Observation_yo.py``::
555
556     #-*-coding:iso-8859-1-*-
557     #
558     import sys ; sys.path.insert(0,"# INSERT PHYSICAL SCRIPT PATH")
559     from Physical_data_and_covariance_matrices import True_state
560     from Physical_simulation_functions import FunctionH
561     #
562     xt, noms = True_state()
563     #
564     yo = FunctionH( xt )
565     #
566     # Creating the required ADAO variable
567     # -----------------------------------
568     Observation = list(yo)
569
570 and the second one named ``Script_ObservationError_R.py``::
571
572     #-*-coding:iso-8859-1-*-
573     #
574     import sys ; sys.path.insert(0,"# INSERT PHYSICAL SCRIPT PATH")
575     from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
576     from Physical_simulation_functions import FunctionH
577     #
578     xt, names = True_state()
579     #
580     yo = FunctionH( xt )
581     #
582     R  = 0.0001 * Simple_Matrix( size = len(yo) )
583     #
584     # Creating the required ADAO variable
585     # -----------------------------------
586     ObservationError = R
587
588 As in previous examples, it can be useful to define some parameters for the data
589 assimilation algorithm. For example, if we use the standard 3DVAR algorithm, the
590 following parameters can be defined in a Python script file named
591 ``Script_AlgorithmParameters.py``::
592
593     #-*-coding:iso-8859-1-*-
594     #
595     # Creating the required ADAO variable
596     # -----------------------------------
597     AlgorithmParameters = {
598         "Minimizer" : "TNC",         # Possible : "LBFGSB", "TNC", "CG", "BFGS"
599         "MaximumNumberOfSteps" : 15, # Number of global iterative steps
600         "Bounds" : [
601             [ None, None ],          # Bound on the first parameter
602             [ 0., 4. ],              # Bound on the second parameter
603             [ 0., None ],            # Bound on the third parameter
604             ],
605     }
606
607 Finally, it is common to post-process the results, retrieving them after the
608 data assimilation phase in order to analyse, print or show them. It requires to
609 use a intermediary Python script file in order to extract these results. The
610 following example Python script file named ``Script_UserPostAnalysis.py``,
611 illustrates the fact::
612
613     #-*-coding:iso-8859-1-*-
614     #
615     import sys ; sys.path.insert(0,"# INSERT PHYSICAL SCRIPT PATH")
616     from Physical_data_and_covariance_matrices import True_state
617     import numpy
618     #
619     xt, names   = True_state()
620     xa          = ADD.get("Analysis").valueserie(-1)
621     x_series    = ADD.get("CurrentState").valueserie()
622     J           = ADD.get("CostFunctionJ").valueserie()
623     #
624     # Verifying the results by printing
625     # ---------------------------------
626     print
627     print "xt = %s"%xt
628     print "xa = %s"%numpy.array(xa)
629     print
630     for i in range( len(x_series) ):
631         print "Step %2i : J = %.5e  et  X = %s"%(i, J[i], x_series[i])
632     print
633
634 At the end, we get a description of the whole case setup through a set of files
635 listed here:
636
637 #.      ``Physical_data_and_covariance_matrices.py``
638 #.      ``Physical_simulation_functions.py``
639 #.      ``Script_AlgorithmParameters.py``
640 #.      ``Script_BackgroundError_B.py``
641 #.      ``Script_Background_xb.py``
642 #.      ``Script_ObservationError_R.py``
643 #.      ``Script_ObservationOperator_H.py``
644 #.      ``Script_Observation_yo.py``
645 #.      ``Script_UserPostAnalysis.py``
646
647 We insist here that all these scripts are written by the user and can not be
648 automatically tested. So the user is required to verify the scripts (and in
649 particular their input/output) in order to limit the difficulty of debug. We
650 recall: **script methodology is not a "safe" procedure, in the sense that
651 erroneous data, or errors in calculations, can be directly injected into the
652 YACS scheme execution.**
653
654 Building the case with external data definition by scripts
655 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
656
657 All these scripts can then be used to define the ADAO case with external data
658 definition by Python script files. It is entirely similar to the method
659 described in the `Building a simple estimation case with external data
660 definition by scripts`_ previous section. For each variable to be defined, we
661 select the "*Script*" option of the "*FROM*" keyword, which leads to a
662 "*SCRIPT_DATA/SCRIPT_FILE*" entry in the tree.
663
664 The other steps to build the ADAO case are exactly the same as in the `Building
665 a simple estimation case with explicit data definition`_ previous section.
666
667 Using the simple linear operator :math:`\mathbf{H}` from the Python script file
668 ``Physical_simulation_functions.py`` in the ADAO examples standard directory,
669 the results will look like::
670
671     xt = [1 2 3]
672     xa = [ 1.000014    2.000458  3.000390]
673
674     Step  0 : J = 1.81750e+03  et  X = [1.014011, 2.459175, 3.390462]
675     Step  1 : J = 1.81750e+03  et  X = [1.014011, 2.459175, 3.390462]
676     Step  2 : J = 1.79734e+01  et  X = [1.010771, 2.040342, 2.961378]
677     Step  3 : J = 1.79734e+01  et  X = [1.010771, 2.040342, 2.961378]
678     Step  4 : J = 1.81909e+00  et  X = [1.000826, 2.000352, 3.000487]
679     Step  5 : J = 1.81909e+00  et  X = [1.000826, 2.000352, 3.000487]
680     Step  6 : J = 1.81641e+00  et  X = [1.000247, 2.000651, 3.000156]
681     Step  7 : J = 1.81641e+00  et  X = [1.000247, 2.000651, 3.000156]
682     Step  8 : J = 1.81569e+00  et  X = [1.000015, 2.000432, 3.000364]
683     Step  9 : J = 1.81569e+00  et  X = [1.000015, 2.000432, 3.000364]
684     Step 10 : J = 1.81568e+00  et  X = [1.000013, 2.000458, 3.000390]
685     ...
686
687 The state at the first step is the randomly generated background state
688 :math:`\mathbf{x}^b`. After completion, these printing on standard output is
689 available in the "*YACS Container Log*", obtained through the right click menu
690 of the "*proc*" window in the YACS scheme.
691
692 .. [#] For more information on YACS, see the the *YACS module User's Guide* available in the main "*Help*" menu of SALOME GUI.