.. include:: snippets/Header2Algo09.rst
-.. include:: scripts/simple_3DVAR.rst
+.. --------- ..
+.. include:: scripts/simple_3DVAR1.rst
-.. literalinclude:: scripts/simple_3DVAR.py
+.. literalinclude:: scripts/simple_3DVAR1.py
.. include:: snippets/Header2Algo10.rst
-.. literalinclude:: scripts/simple_3DVAR.res
+.. literalinclude:: scripts/simple_3DVAR1.res
+ :language: none
.. include:: snippets/Header2Algo11.rst
-.. _simple_3DVAR:
-.. image:: scripts/simple_3DVAR.png
+.. _simple_3DVAR1:
+.. image:: scripts/simple_3DVAR1.png
+ :align: center
+ :width: 90%
+
+.. include:: scripts/simple_3DVAR1Plus.rst
+
+.. _simple_3DVAR1Plus:
+.. image:: scripts/simple_3DVAR1Plus.png
+ :align: center
+ :width: 90%
+
+.. --------- ..
+.. include:: scripts/simple_3DVAR2.rst
+
+.. literalinclude:: scripts/simple_3DVAR2.py
+
+.. include:: snippets/Header2Algo10.rst
+
+.. literalinclude:: scripts/simple_3DVAR2.res
+ :language: none
+
+.. include:: snippets/Header2Algo11.rst
+
+.. _simple_3DVAR2_state:
+.. image:: scripts/simple_3DVAR2_state.png
+ :align: center
+ :width: 90%
+
+.. simple_3DVAR2_variance:
+.. image:: scripts/simple_3DVAR2_variance.png
+ :align: center
+ :width: 90%
+
+.. --------- ..
+.. include:: scripts/simple_3DVAR3.rst
+
+.. literalinclude:: scripts/simple_3DVAR3.py
+
+.. include:: snippets/Header2Algo10.rst
+
+.. literalinclude:: scripts/simple_3DVAR3.res
+ :language: none
+
+.. include:: snippets/Header2Algo11.rst
+
+.. _simple_3DVAR3_state:
+.. image:: scripts/simple_3DVAR3_state.png
+ :align: center
+ :width: 90%
+
+.. _simple_3DVAR3_variance:
+.. image:: scripts/simple_3DVAR3_variance.png
:align: center
:width: 90%
- :ref:`section_ref_algorithm_Blue`
- :ref:`section_ref_algorithm_ExtendedBlue`
+- :ref:`section_ref_algorithm_KalmanFilter`
- :ref:`section_ref_algorithm_LinearityTest`
.. ------------------------------------ ..
+++ /dev/null
-# -*- coding: utf-8 -*-
-#
-from numpy import array, ravel
-def QuadFunction( coefficients ):
- """
- Quadratic simulation in x: y = a x^2 + b x + c
- """
- a, b, c = list(ravel(coefficients))
- x_points = (-5, 0, 1, 3, 10)
- y_points = []
- for x in x_points:
- y_points.append( a*x*x + b*x + c )
- return array(y_points)
-#
-Xb = array([1., 1., 1.])
-Yobs = array([57, 2, 3, 17, 192])
-#
-print("Resolution of the calibration problem")
-print("-------------------------------------")
-print("")
-from adao import adaoBuilder
-case = adaoBuilder.New()
-case.setBackground( Vector = Xb, Stored=True )
-case.setBackgroundError( ScalarSparseMatrix = 1.e6 )
-case.setObservation( Vector = Yobs, Stored=True )
-case.setObservationError( ScalarSparseMatrix = 1. )
-case.setObservationOperator( OneFunction = QuadFunction )
-case.setAlgorithmParameters(
- Algorithm='3DVAR',
- Parameters={
- 'MaximumNumberOfIterations': 100,
- 'StoreSupplementaryCalculations': [
- 'CurrentState',
- ],
- },
- )
-case.setObserver(
- Info=" Intermediate state at the current iteration:",
- Template='ValuePrinter',
- Variable='CurrentState',
- )
-case.execute()
-print("")
-#
-#-------------------------------------------------------------------------------
-#
-print("Calibration of %i coefficients in a 1D quadratic function on %i measures"%(
- len(case.get('Background')),
- len(case.get('Observation')),
- ))
-print("----------------------------------------------------------------------")
-print("")
-print("Observation vector.................:", ravel(case.get('Observation')))
-print("A priori background state..........:", ravel(case.get('Background')))
-print("")
-print("Expected theoretical coefficients..:", ravel((2,-1,2)))
-print("")
-print("Number of iterations...............:", len(case.get('CurrentState')))
-print("Number of simulations..............:", len(case.get('CurrentState'))*4)
-print("Calibration resulting coefficients.:", ravel(case.get('Analysis')[-1]))
-#
-Xa = case.get('Analysis')[-1]
-import matplotlib.pyplot as plt
-plt.rcParams['figure.figsize'] = (10, 4)
-#
-plt.figure()
-plt.plot((-5,0,1,3,10),QuadFunction(Xb),'b-',label="Simulation at background")
-plt.plot((-5,0,1,3,10),Yobs, 'kX',label='Observation',markersize=10)
-plt.plot((-5,0,1,3,10),QuadFunction(Xa),'r-',label="Simulation at optimum")
-plt.legend()
-plt.title('Coefficients calibration', fontweight='bold')
-plt.xlabel('Arbitrary coordinate')
-plt.ylabel('Observations')
-plt.savefig("simple_3DVAR.png")
+++ /dev/null
-Resolution of the calibration problem
--------------------------------------
-
- Intermediate state at the current iteration: [1. 1. 1.]
- Intermediate state at the current iteration: [1.99739508 1.07086406 1.01346638]
- Intermediate state at the current iteration: [1.83891966 1.04815981 1.01208385]
- Intermediate state at the current iteration: [1.8390702 1.03667176 1.01284797]
- Intermediate state at the current iteration: [1.83967236 0.99071957 1.01590445]
- Intermediate state at the current iteration: [1.84208099 0.8069108 1.02813037]
- Intermediate state at the current iteration: [ 1.93711599 -0.56383145 1.12097995]
- Intermediate state at the current iteration: [ 1.99838848 -1.00480573 1.1563713 ]
- Intermediate state at the current iteration: [ 2.0135905 -1.04815933 1.16155285]
- Intermediate state at the current iteration: [ 2.01385679 -1.03874809 1.16129657]
- Intermediate state at the current iteration: [ 2.01377856 -1.03700044 1.16157611]
- Intermediate state at the current iteration: [ 2.01338902 -1.02943736 1.16528951]
- Intermediate state at the current iteration: [ 2.01265633 -1.0170847 1.17793974]
- Intermediate state at the current iteration: [ 2.0112487 -0.99745509 1.21485091]
- Intermediate state at the current iteration: [ 2.00863696 -0.96943284 1.30917045]
- Intermediate state at the current iteration: [ 2.00453385 -0.94011716 1.51021882]
- Intermediate state at the current iteration: [ 2.00013539 -0.93313893 1.80539433]
- Intermediate state at the current iteration: [ 1.95437244 -0.76890418 2.04566856]
- Intermediate state at the current iteration: [ 1.99797362 -0.92538074 1.81674451]
- Intermediate state at the current iteration: [ 1.99760514 -0.95929669 2.01402091]
- Intermediate state at the current iteration: [ 1.99917565 -0.99152672 2.03171791]
- Intermediate state at the current iteration: [ 1.99990376 -0.99963123 2.00671578]
- Intermediate state at the current iteration: [ 1.99999841 -1.00005285 2.00039699]
- Intermediate state at the current iteration: [ 2.00000014 -1.00000307 2.00000221]
- Intermediate state at the current iteration: [ 2. -0.99999992 1.99999987]
-
-Calibration of 3 coefficients in a 1D quadratic function on 5 measures
-----------------------------------------------------------------------
-
-Observation vector.................: [ 57. 2. 3. 17. 192.]
-A priori background state..........: [1. 1. 1.]
-
-Expected theoretical coefficients..: [ 2 -1 2]
-
-Number of iterations...............: 25
-Number of simulations..............: 100
-Calibration resulting coefficients.: [ 2. -0.99999992 1.99999987]
+++ /dev/null
-.. index:: single: 3DVAR (example)
-
-This example describes the calibration of parameters :math:`\mathbf{x}` of a
-quadratic observation model :math:`H`. This model is here represented as a
-function named ``QuadFunction``. This function get as input the coefficients
-vector :math:`\mathbf{x}`, and return as output the evaluation vector
-:math:`\mathbf{y}` of the quadratic model at the predefined internal control
-points. The calibration is done using an initial coefficient set (background
-state specified by ``Xb`` in the code), and with the information
-:math:`\mathbf{y}^o` (specified by ``Yobs`` in the code) of 5 measures obtained
-in these same internal control points. We set twin experiments (see
-:ref:`section_methodology_twin`) and the measurements are supposed to be
-perfect. We choose to emphasize the observations versus the background by
-setting a great variance for the background error, here of :math:`10^{6}`.
-
-The adjustment is carried out by displaying intermediate results during
-iterative optimization.
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+from numpy import array, ravel
+def QuadFunction( coefficients ):
+ """
+ Quadratic simulation in x: y = a x^2 + b x + c
+ """
+ a, b, c = list(ravel(coefficients))
+ x_points = (-5, 0, 1, 3, 10)
+ y_points = []
+ for x in x_points:
+ y_points.append( a*x*x + b*x + c )
+ return array(y_points)
+#
+Xb = array([1., 1., 1.])
+Yobs = array([57, 2, 3, 17, 192])
+#
+print("Resolution of the calibration problem")
+print("-------------------------------------")
+print("")
+from adao import adaoBuilder
+case = adaoBuilder.New()
+case.setBackground( Vector = Xb, Stored=True )
+case.setBackgroundError( ScalarSparseMatrix = 1.e6 )
+case.setObservation( Vector = Yobs, Stored=True )
+case.setObservationError( ScalarSparseMatrix = 1. )
+case.setObservationOperator( OneFunction = QuadFunction )
+case.setAlgorithmParameters(
+ Algorithm='3DVAR',
+ Parameters={
+ 'MaximumNumberOfIterations': 100,
+ 'StoreSupplementaryCalculations': [
+ 'CurrentState',
+ ],
+ },
+ )
+case.setObserver(
+ Info=" Intermediate state at the current iteration:",
+ Template='ValuePrinter',
+ Variable='CurrentState',
+ )
+case.execute()
+print("")
+#
+#-------------------------------------------------------------------------------
+#
+print("Calibration of %i coefficients in a 1D quadratic function on %i measures"%(
+ len(case.get('Background')),
+ len(case.get('Observation')),
+ ))
+print("----------------------------------------------------------------------")
+print("")
+print("Observation vector.................:", ravel(case.get('Observation')))
+print("A priori background state..........:", ravel(case.get('Background')))
+print("")
+print("Expected theoretical coefficients..:", ravel((2,-1,2)))
+print("")
+print("Number of iterations...............:", len(case.get('CurrentState')))
+print("Number of simulations..............:", len(case.get('CurrentState'))*4)
+print("Calibration resulting coefficients.:", ravel(case.get('Analysis')[-1]))
+#
+Xa = case.get('Analysis')[-1]
+import matplotlib.pyplot as plt
+plt.rcParams['figure.figsize'] = (10, 4)
+#
+plt.figure()
+plt.plot((-5,0,1,3,10),QuadFunction(Xb),'b--',label="Simulation at background")
+plt.plot((-5,0,1,3,10),Yobs, 'kX', label='Observation',markersize=10)
+plt.plot((-5,0,1,3,10),QuadFunction(Xa),'r-', label="Simulation at optimum")
+plt.legend()
+plt.title('Coefficients calibration', fontweight='bold')
+plt.xlabel('Arbitrary coordinate')
+plt.ylabel('Observations')
+plt.savefig("simple_3DVAR1.png")
--- /dev/null
+Resolution of the calibration problem
+-------------------------------------
+
+ Intermediate state at the current iteration: [1. 1. 1.]
+ Intermediate state at the current iteration: [1.99739508 1.07086406 1.01346638]
+ Intermediate state at the current iteration: [1.83891966 1.04815981 1.01208385]
+ Intermediate state at the current iteration: [1.8390702 1.03667176 1.01284797]
+ Intermediate state at the current iteration: [1.83967236 0.99071957 1.01590445]
+ Intermediate state at the current iteration: [1.84208099 0.8069108 1.02813037]
+ Intermediate state at the current iteration: [ 1.93711599 -0.56383145 1.12097995]
+ Intermediate state at the current iteration: [ 1.99838848 -1.00480573 1.1563713 ]
+ Intermediate state at the current iteration: [ 2.0135905 -1.04815933 1.16155285]
+ Intermediate state at the current iteration: [ 2.01385679 -1.03874809 1.16129657]
+ Intermediate state at the current iteration: [ 2.01377856 -1.03700044 1.16157611]
+ Intermediate state at the current iteration: [ 2.01338902 -1.02943736 1.16528951]
+ Intermediate state at the current iteration: [ 2.01265633 -1.0170847 1.17793974]
+ Intermediate state at the current iteration: [ 2.0112487 -0.99745509 1.21485091]
+ Intermediate state at the current iteration: [ 2.00863696 -0.96943284 1.30917045]
+ Intermediate state at the current iteration: [ 2.00453385 -0.94011716 1.51021882]
+ Intermediate state at the current iteration: [ 2.00013539 -0.93313893 1.80539433]
+ Intermediate state at the current iteration: [ 1.95437244 -0.76890418 2.04566856]
+ Intermediate state at the current iteration: [ 1.99797362 -0.92538074 1.81674451]
+ Intermediate state at the current iteration: [ 1.99760514 -0.95929669 2.01402091]
+ Intermediate state at the current iteration: [ 1.99917565 -0.99152672 2.03171791]
+ Intermediate state at the current iteration: [ 1.99990376 -0.99963123 2.00671578]
+ Intermediate state at the current iteration: [ 1.99999841 -1.00005285 2.00039699]
+ Intermediate state at the current iteration: [ 2.00000014 -1.00000307 2.00000221]
+ Intermediate state at the current iteration: [ 2. -0.99999992 1.99999987]
+
+Calibration of 3 coefficients in a 1D quadratic function on 5 measures
+----------------------------------------------------------------------
+
+Observation vector.................: [ 57. 2. 3. 17. 192.]
+A priori background state..........: [1. 1. 1.]
+
+Expected theoretical coefficients..: [ 2 -1 2]
+
+Number of iterations...............: 25
+Number of simulations..............: 100
+Calibration resulting coefficients.: [ 2. -0.99999992 1.99999987]
--- /dev/null
+.. index:: single: 3DVAR (example)
+
+First example
+.............
+
+This example describes the calibration of parameters :math:`\mathbf{x}` of a
+quadratic observation model :math:`H`. This model :math:`H` is here represented
+as a function named ``QuadFunction`` for the purpose of this example. This
+function get as input the coefficients vector :math:`\mathbf{x}` of the
+quadratic form, and return as output the evaluation vector :math:`\mathbf{y}`
+of the quadratic model at the predefined internal control points, predefined in
+a static way in the model.
+
+The calibration is done using an initial coefficient set (background state
+specified by ``Xb`` in the example), and with the information
+:math:`\mathbf{y}^o` (specified by ``Yobs`` in the example) of 5 measures
+obtained in these same internal control points. We set twin experiments (see
+:ref:`section_methodology_twin`) and the measurements are supposed to be
+perfect. We choose to emphasize the observations, versus the background, by
+setting artificially a great variance for the background error, here of
+:math:`10^{6}`.
+
+The adjustment is carried out by displaying intermediate results using
+"*observer*" (for reference, see :ref:`section_advanced_observer`) during
+iterative optimization.
--- /dev/null
+In the figure, the lines simply join the simulated values at the control points
+and make them easier to read.
+
+We see, in this particularly simple case, that the optimum is identical to the
+values of coefficients ``[2, -1, 2]`` used to generate the pseudo-observations
+for the twin experiment. Naturally, the final simulation obtained with the
+adjusted coefficients coincides with the observations at the 5 control points.
+
+To illustrate more clearly the recalibration performed, using the same
+information, we can plot the continuous version of the quadratic curves
+available for the background state of the coefficients, and for the optimal
+state of the coefficients, together with the values obtained at the control
+points. In addition, to improve readability, the lines joining the simulations
+in the previous graph have been removed.
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+from numpy import array, random
+random.seed(1234567)
+Xtrue = -0.37727
+Yobs = []
+for i in range(51):
+ Yobs.append([random.normal(Xtrue, 0.1, size=(1,)),])
+#
+print("Variational estimation of a time trajectory")
+print("-------------------------------------------")
+print(" Noisy measurements acquired on %i time steps"%(len(Yobs)-1,))
+print("")
+from adao import adaoBuilder
+case = adaoBuilder.New()
+#
+case.setBackground (Vector = [0.])
+case.setBackgroundError (ScalarSparseMatrix = 0.1**2)
+#
+case.setObservationOperator(Matrix = [1.])
+case.setObservation (VectorSerie = Yobs)
+case.setObservationError (ScalarSparseMatrix = 0.3**2)
+#
+case.setEvolutionModel (Matrix = [1.])
+case.setEvolutionError (ScalarSparseMatrix = 1e-5)
+#
+case.setAlgorithmParameters(
+ Algorithm="3DVAR",
+ Parameters={
+ "EstimationOf":"State",
+ "StoreSupplementaryCalculations":[
+ "Analysis",
+ "APosterioriCovariance",
+ ],
+ },
+ )
+#
+case.execute()
+Xa = case.get("Analysis")
+Pa = case.get("APosterioriCovariance")
+#
+print(" Analyzed state at final observation:", Xa[-1])
+print("")
+print(" Final a posteriori variance:", Pa[-1])
+print("")
+#
+#-------------------------------------------------------------------------------
+#
+Observations = array([yo[0] for yo in Yobs])
+Estimates = array([xa[0] for xa in case.get("Analysis")])
+Variances = array([pa[0,0] for pa in case.get("APosterioriCovariance")])
+#
+import matplotlib.pyplot as plt
+plt.rcParams['figure.figsize'] = (10, 4)
+#
+plt.figure()
+plt.plot(Observations,'kx',label='Noisy measurements')
+plt.plot(Estimates,'r-',label='Estimated state')
+plt.axhline(Xtrue,color='b',label='Truth value')
+plt.legend()
+plt.title('Estimate of the state', fontweight='bold')
+plt.xlabel('Observation step')
+plt.ylabel('Voltage')
+plt.savefig("simple_3DVAR2_state.png")
+#
+plt.figure()
+iobs = range(1,len(Observations))
+plt.plot(iobs,Variances[iobs],label='A posteriori error variance')
+plt.title('Estimate of the a posteriori error variance', fontweight='bold')
+plt.xlabel('Observation step')
+plt.ylabel('$(Voltage)^2$')
+plt.setp(plt.gca(),'ylim',[0,.01])
+plt.savefig("simple_3DVAR2_variance.png")
--- /dev/null
+Variational estimation of a time trajectory
+-------------------------------------------
+ Noisy measurements acquired on 50 time steps
+
+ Analyzed state at final observation: [-0.37110687]
+
+ Final a posteriori variance: [[0.009]]
+
--- /dev/null
+Second example
+..............
+
+The 3DVAR can also be used for a **time analysis of the observations of a given
+dynamic model**. In this case, the analysis is performed iteratively, at the
+arrival of each observation. For this example, we use the same simple dynamic
+system [Welch06]_ that is analyzed in the Kalman Filter
+:ref:`section_ref_algorithm_KalmanFilter_examples`.
+
+At each step, the classical 3DVAR analysis updates only the state of the
+system. By modifying the *a priori* covariance values with respect to the
+initial assumptions of the filtering, this 3DVAR reanalysis allows to converge
+towards the true trajectory, as illustrated in the associated figure, in a
+slightly slower speed than with a Kalman Filter.
+
+.. note::
+
+ Note about *a posteriori* covariances: classically, the 3DVAR iterative
+ analysis updates only the state and not its covariance. As the assumptions
+ of operators and *a priori* covariance remain unchanged here during the
+ evolution, the *a posteriori* covariance is constant. The following plot of
+ this *a posteriori* covariance allows us to insist on this property, which
+ is entirely expected from the 3DVAR analysis. A more advanced hypothesis is
+ proposed in the forthcoming example.
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+from numpy import array, random
+random.seed(1234567)
+Xtrue = -0.37727
+Yobs = []
+for i in range(51):
+ Yobs.append([random.normal(Xtrue, 0.1, size=(1,)),])
+#
+print("Variational estimation of a time trajectory")
+print("-------------------------------------------")
+print(" Noisy measurements acquired on %i time steps"%(len(Yobs)-1,))
+print("")
+from adao import adaoBuilder
+case = adaoBuilder.New()
+#
+case.setBackground (Vector = [0.])
+#
+case.setObservationOperator(Matrix = [1.])
+case.setObservationError (ScalarSparseMatrix = 0.3**2)
+#
+case.setEvolutionModel (Matrix = [1.])
+case.setEvolutionError (ScalarSparseMatrix = 1e-5)
+#
+case.setAlgorithmParameters(
+ Algorithm="3DVAR",
+ Parameters={
+ "EstimationOf":"State",
+ "StoreSupplementaryCalculations":[
+ "Analysis",
+ "APosterioriCovariance",
+ ],
+ },
+ )
+#
+# Loop to obtain an analysis at each observation arrival
+#
+Ca = 1.
+for i in range(1,len(Yobs)):
+ case.setObservation(Vector = Yobs[i])
+ case.setBackgroundError(ScalarSparseMatrix = Ca) # Update of the covariance
+ case.execute( nextStep = True )
+ #
+ # Hypothesis of the a priori covariance decay of the background error
+ #
+ if float(Ca) <= 0.1**2:
+ Ca = 0.1**2
+ else:
+ Ca = Ca * 0.9**2
+#
+Xa = case.get("Analysis")
+Pa = case.get("APosterioriCovariance")
+#
+print(" Analyzed state at final observation:", Xa[-1])
+print("")
+print(" Final a posteriori variance:", Pa[-1])
+print("")
+#
+#-------------------------------------------------------------------------------
+#
+Observations = array([yo[0] for yo in Yobs])
+Estimates = array([xa[0] for xa in case.get("Analysis")])
+Variances = array([pa[0,0] for pa in case.get("APosterioriCovariance")])
+#
+import matplotlib.pyplot as plt
+plt.rcParams['figure.figsize'] = (10, 4)
+#
+plt.figure()
+plt.plot(Observations,'kx',label='Noisy measurements')
+plt.plot(Estimates,'r-',label='Estimated state')
+plt.axhline(Xtrue,color='b',label='Truth value')
+plt.legend()
+plt.title('Estimate of the state', fontweight='bold')
+plt.xlabel('Observation step')
+plt.ylabel('Voltage')
+plt.savefig("simple_3DVAR3_state.png")
+#
+plt.figure()
+iobs = range(1,len(Observations))
+plt.plot(iobs,Variances[iobs],label='A posteriori error variance')
+plt.title('Estimate of the a posteriori error variance', fontweight='bold')
+plt.xlabel('Observation step')
+plt.ylabel('$(Voltage)^2$')
+plt.savefig("simple_3DVAR3_variance.png")
--- /dev/null
+Variational estimation of a time trajectory
+-------------------------------------------
+ Noisy measurements acquired on 50 time steps
+
+ Analyzed state at final observation: [-0.37334336]
+
+ Final a posteriori variance: [[0.009]]
+
--- /dev/null
+Third example
+.............
+
+From the preceding example, if one wants to adapt the time convergence of the
+3DVAR, one can change, for example, the *a priori* covariance assumptions of
+the background errors during the iterations. This update is an **assumption**
+of the user, and there are multiple alternatives that will depend on the
+physics of the case. We illustrate one of them here.
+
+We choose, in an arbitrary way, to make the *a priori* covariance of the
+background errors to decrease by a constant factor :math:`0.9^2=0.81` as long
+as it remains above a limit value of :math:`0.1^2=0.01` (which is the fixed
+value of *a priori* covariance of the background errors of the previous
+example), knowing that it starts at the value `1` (which is the fixed value of
+*a priori* covariance of the background errors used for the first step of
+Kalman filtering). This value is updated at each step, by reinjecting it as the
+*a priori* covariance of the state which is used as a background in the next
+step of analysis, in an explicit loop.
+
+We notice in this case that the state estimation converges faster to the true
+value, and that the assimilation then behaves similarly to the examples for the
+Kalman Filter, or to the previous example with the manually adapted *a priori*
+covariances. Moreover, the *a posteriori* covariance decreases as long as we
+force the decrease of the *a priori* covariance.
+
+.. note::
+
+ We insist on the fact that the *a priori* covariance variations, which
+ determine the *a posteriori* covariance variations, are a **user
+ assumption** and not an obligation. This assumption must therefore be
+ **adapted to the physical case**.
Python (TUI) use examples
+++++++++++++++++++++++++
-Here is a very simple use of the given algorithm and its parameters, written in
-:ref:`section_tui`, and from which input information allow to define an
-equivalent case in graphical user interface.
+Here is one or more very simple examples of the proposed algorithm and its
+parameters, written in :ref:`section_tui`. Moreover, the information given as
+input also allows to define an equivalent case in :ref:`section_gui_in_salome`.
Python, 3.6.5, 3.10.9
Numpy, 1.14.3, 1.24.2
- Scipy, 1.1.0, 1.10.0
+ Scipy, 0.19.1, 1.10.1
MatplotLib, 2.2.2, 3.7.0
GnuplotPy, 1.8, 1.8
NLopt, 2.4.2, 2.7.1
.. include:: snippets/Header2Algo09.rst
-.. include:: scripts/simple_3DVAR.rst
+.. --------- ..
+.. include:: scripts/simple_3DVAR1.rst
-.. literalinclude:: scripts/simple_3DVAR.py
+.. literalinclude:: scripts/simple_3DVAR1.py
.. include:: snippets/Header2Algo10.rst
-.. literalinclude:: scripts/simple_3DVAR.res
+.. literalinclude:: scripts/simple_3DVAR1.res
:language: none
.. include:: snippets/Header2Algo11.rst
-.. _simple_3DVAR:
-.. image:: scripts/simple_3DVAR.png
+.. _simple_3DVAR1:
+.. image:: scripts/simple_3DVAR1.png
+ :align: center
+ :width: 90%
+
+.. include:: scripts/simple_3DVAR1Plus.rst
+
+.. _simple_3DVAR1Plus:
+.. image:: scripts/simple_3DVAR1Plus.png
+ :align: center
+ :width: 90%
+
+.. --------- ..
+.. include:: scripts/simple_3DVAR2.rst
+
+.. literalinclude:: scripts/simple_3DVAR2.py
+
+.. include:: snippets/Header2Algo10.rst
+
+.. literalinclude:: scripts/simple_3DVAR2.res
+ :language: none
+
+.. include:: snippets/Header2Algo11.rst
+
+.. _simple_3DVAR2_state:
+.. image:: scripts/simple_3DVAR2_state.png
+ :align: center
+ :width: 90%
+
+.. simple_3DVAR2_variance:
+.. image:: scripts/simple_3DVAR2_variance.png
+ :align: center
+ :width: 90%
+
+.. --------- ..
+.. include:: scripts/simple_3DVAR3.rst
+
+.. literalinclude:: scripts/simple_3DVAR3.py
+
+.. include:: snippets/Header2Algo10.rst
+
+.. literalinclude:: scripts/simple_3DVAR3.res
+ :language: none
+
+.. include:: snippets/Header2Algo11.rst
+
+.. _simple_3DVAR3_state:
+.. image:: scripts/simple_3DVAR3_state.png
+ :align: center
+ :width: 90%
+
+.. _simple_3DVAR3_variance:
+.. image:: scripts/simple_3DVAR3_variance.png
:align: center
:width: 90%
- :ref:`section_ref_algorithm_Blue`
- :ref:`section_ref_algorithm_ExtendedBlue`
+- :ref:`section_ref_algorithm_KalmanFilter`
- :ref:`section_ref_algorithm_LinearityTest`
.. ------------------------------------ ..
+++ /dev/null
-# -*- coding: utf-8 -*-
-#
-from numpy import array, ravel
-def QuadFunction( coefficients ):
- """
- Simulation quadratique aux points x : y = a x^2 + b x + c
- """
- a, b, c = list(ravel(coefficients))
- x_points = (-5, 0, 1, 3, 10)
- y_points = []
- for x in x_points:
- y_points.append( a*x*x + b*x + c )
- return array(y_points)
-#
-Xb = array([1., 1., 1.])
-Yobs = array([57, 2, 3, 17, 192])
-#
-print("Résolution du problème de calage")
-print("--------------------------------")
-print("")
-from adao import adaoBuilder
-case = adaoBuilder.New()
-case.setBackground( Vector = Xb, Stored=True )
-case.setBackgroundError( ScalarSparseMatrix = 1.e6 )
-case.setObservation( Vector = Yobs, Stored=True )
-case.setObservationError( ScalarSparseMatrix = 1. )
-case.setObservationOperator( OneFunction = QuadFunction )
-case.setAlgorithmParameters(
- Algorithm='3DVAR',
- Parameters={
- 'MaximumNumberOfIterations': 100,
- 'StoreSupplementaryCalculations': [
- 'CurrentState',
- ],
- },
- )
-case.setObserver(
- Info=" État intermédiaire en itération courante :",
- Template='ValuePrinter',
- Variable='CurrentState',
- )
-case.execute()
-print("")
-#
-#-------------------------------------------------------------------------------
-#
-print("Calage de %i coefficients pour une forme quadratique 1D sur %i mesures"%(
- len(case.get('Background')),
- len(case.get('Observation')),
- ))
-print("--------------------------------------------------------------------")
-print("")
-print("Vecteur d'observation.............:", ravel(case.get('Observation')))
-print("État d'ébauche a priori...........:", ravel(case.get('Background')))
-print("")
-print("Coefficients théoriques attendus..:", ravel((2,-1,2)))
-print("")
-print("Nombre d'itérations...............:", len(case.get('CurrentState')))
-print("Nombre de simulations.............:", len(case.get('CurrentState'))*4)
-print("Coefficients résultants du calage.:", ravel(case.get('Analysis')[-1]))
-#
-Xa = case.get('Analysis')[-1]
-import matplotlib.pyplot as plt
-plt.rcParams['figure.figsize'] = (10, 4)
-#
-plt.figure()
-plt.plot((-5,0,1,3,10),QuadFunction(Xb),'b-',label="Simulation à l'ébauche")
-plt.plot((-5,0,1,3,10),Yobs, 'kX',label='Observation',markersize=10)
-plt.plot((-5,0,1,3,10),QuadFunction(Xa),'r-',label="Simulation à l'optimum")
-plt.legend()
-plt.title('Calage de coefficients', fontweight='bold')
-plt.xlabel('Coordonnée arbitraire')
-plt.ylabel('Observations')
-plt.savefig("simple_3DVAR.png")
+++ /dev/null
-Résolution du problème de calage
---------------------------------
-
- État intermédiaire en itération courante : [1. 1. 1.]
- État intermédiaire en itération courante : [1.99739508 1.07086406 1.01346638]
- État intermédiaire en itération courante : [1.83891966 1.04815981 1.01208385]
- État intermédiaire en itération courante : [1.8390702 1.03667176 1.01284797]
- État intermédiaire en itération courante : [1.83967236 0.99071957 1.01590445]
- État intermédiaire en itération courante : [1.84208099 0.8069108 1.02813037]
- État intermédiaire en itération courante : [ 1.93711599 -0.56383145 1.12097995]
- État intermédiaire en itération courante : [ 1.99838848 -1.00480573 1.1563713 ]
- État intermédiaire en itération courante : [ 2.0135905 -1.04815933 1.16155285]
- État intermédiaire en itération courante : [ 2.01385679 -1.03874809 1.16129657]
- État intermédiaire en itération courante : [ 2.01377856 -1.03700044 1.16157611]
- État intermédiaire en itération courante : [ 2.01338902 -1.02943736 1.16528951]
- État intermédiaire en itération courante : [ 2.01265633 -1.0170847 1.17793974]
- État intermédiaire en itération courante : [ 2.0112487 -0.99745509 1.21485091]
- État intermédiaire en itération courante : [ 2.00863696 -0.96943284 1.30917045]
- État intermédiaire en itération courante : [ 2.00453385 -0.94011716 1.51021882]
- État intermédiaire en itération courante : [ 2.00013539 -0.93313893 1.80539433]
- État intermédiaire en itération courante : [ 1.95437244 -0.76890418 2.04566856]
- État intermédiaire en itération courante : [ 1.99797362 -0.92538074 1.81674451]
- État intermédiaire en itération courante : [ 1.99760514 -0.95929669 2.01402091]
- État intermédiaire en itération courante : [ 1.99917565 -0.99152672 2.03171791]
- État intermédiaire en itération courante : [ 1.99990376 -0.99963123 2.00671578]
- État intermédiaire en itération courante : [ 1.99999841 -1.00005285 2.00039699]
- État intermédiaire en itération courante : [ 2.00000014 -1.00000307 2.00000221]
- État intermédiaire en itération courante : [ 2. -0.99999992 1.99999987]
-
-Calage de 3 coefficients pour une forme quadratique 1D sur 5 mesures
---------------------------------------------------------------------
-
-Vecteur d'observation.............: [ 57. 2. 3. 17. 192.]
-État d'ébauche a priori...........: [1. 1. 1.]
-
-Coefficients théoriques attendus..: [ 2 -1 2]
-
-Nombre d'itérations...............: 25
-Nombre de simulations.............: 100
-Coefficients résultants du calage.: [ 2. -0.99999992 1.99999987]
+++ /dev/null
-.. index:: single: 3DVAR (exemple)
-
-Cet exemple décrit le recalage des paramètres :math:`\mathbf{x}` d'un modèle
-d'observation :math:`H` quadratique. Ce modèle est représenté ici comme une
-fonction nommée ``QuadFunction``. Cette fonction accepte en entrée le vecteur
-de coefficients :math:`\mathbf{x}`, et fournit en sortie le vecteur
-:math:`\mathbf{y}` d'évaluation du modèle quadratique aux points de contrôle
-internes prédéfinis dans le modèle. Le calage s'effectue sur la base d'un jeu
-initial de coefficients (état d'ébauche désigné par ``Xb`` dans l'exemple), et
-avec l'information :math:`\mathbf{y}^o` (désignée par ``Yobs`` dans l'exemple)
-de 5 mesures obtenues à ces mêmes points de contrôle internes. On se place en
-expériences jumelles (voir :ref:`section_methodology_twin`) et les mesures sont
-parfaites. On privilégie les observations au détriment de l'ébauche par
-l'indication d'une très importante variance d'erreur d'ébauche, ici de
-:math:`10^{6}`.
-
-L'ajustement s'effectue en affichant des résultats intermédiaires lors de
-l'optimisation itérative.
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+from numpy import array, ravel
+def QuadFunction( coefficients ):
+ """
+ Simulation quadratique aux points x : y = a x^2 + b x + c
+ """
+ a, b, c = list(ravel(coefficients))
+ x_points = (-5, 0, 1, 3, 10)
+ y_points = []
+ for x in x_points:
+ y_points.append( a*x*x + b*x + c )
+ return array(y_points)
+#
+Xb = array([1., 1., 1.])
+Yobs = array([57, 2, 3, 17, 192])
+#
+print("Résolution du problème de calage")
+print("--------------------------------")
+print("")
+from adao import adaoBuilder
+case = adaoBuilder.New()
+case.setBackground( Vector = Xb, Stored=True )
+case.setBackgroundError( ScalarSparseMatrix = 1.e6 )
+case.setObservation( Vector = Yobs, Stored=True )
+case.setObservationError( ScalarSparseMatrix = 1. )
+case.setObservationOperator( OneFunction = QuadFunction )
+case.setAlgorithmParameters(
+ Algorithm='3DVAR',
+ Parameters={
+ 'MaximumNumberOfIterations': 100,
+ 'StoreSupplementaryCalculations': [
+ 'CurrentState',
+ ],
+ },
+ )
+case.setObserver(
+ Info=" État intermédiaire en itération courante :",
+ Template='ValuePrinter',
+ Variable='CurrentState',
+ )
+case.execute()
+print("")
+#
+#-------------------------------------------------------------------------------
+#
+print("Calage de %i coefficients pour une forme quadratique 1D sur %i mesures"%(
+ len(case.get('Background')),
+ len(case.get('Observation')),
+ ))
+print("--------------------------------------------------------------------")
+print("")
+print("Vecteur d'observation.............:", ravel(case.get('Observation')))
+print("État d'ébauche a priori...........:", ravel(case.get('Background')))
+print("")
+print("Coefficients théoriques attendus..:", ravel((2,-1,2)))
+print("")
+print("Nombre d'itérations...............:", len(case.get('CurrentState')))
+print("Nombre de simulations.............:", len(case.get('CurrentState'))*4)
+print("Coefficients résultants du calage.:", ravel(case.get('Analysis')[-1]))
+#
+Xa = case.get('Analysis')[-1]
+import matplotlib.pyplot as plt
+plt.rcParams['figure.figsize'] = (10, 4)
+#
+plt.figure()
+plt.plot((-5,0,1,3,10),QuadFunction(Xb),'b--',label="Simulation à l'ébauche")
+plt.plot((-5,0,1,3,10),Yobs, 'kX', label='Observation',markersize=10)
+plt.plot((-5,0,1,3,10),QuadFunction(Xa),'r-', label="Simulation à l'optimum")
+plt.legend()
+plt.title('Calage de coefficients', fontweight='bold')
+plt.xlabel('Coordonnée arbitraire')
+plt.ylabel('Observations')
+plt.savefig("simple_3DVAR1.png")
--- /dev/null
+Résolution du problème de calage
+--------------------------------
+
+ État intermédiaire en itération courante : [1. 1. 1.]
+ État intermédiaire en itération courante : [1.99739508 1.07086406 1.01346638]
+ État intermédiaire en itération courante : [1.83891966 1.04815981 1.01208385]
+ État intermédiaire en itération courante : [1.8390702 1.03667176 1.01284797]
+ État intermédiaire en itération courante : [1.83967236 0.99071957 1.01590445]
+ État intermédiaire en itération courante : [1.84208099 0.8069108 1.02813037]
+ État intermédiaire en itération courante : [ 1.93711599 -0.56383145 1.12097995]
+ État intermédiaire en itération courante : [ 1.99838848 -1.00480573 1.1563713 ]
+ État intermédiaire en itération courante : [ 2.0135905 -1.04815933 1.16155285]
+ État intermédiaire en itération courante : [ 2.01385679 -1.03874809 1.16129657]
+ État intermédiaire en itération courante : [ 2.01377856 -1.03700044 1.16157611]
+ État intermédiaire en itération courante : [ 2.01338902 -1.02943736 1.16528951]
+ État intermédiaire en itération courante : [ 2.01265633 -1.0170847 1.17793974]
+ État intermédiaire en itération courante : [ 2.0112487 -0.99745509 1.21485091]
+ État intermédiaire en itération courante : [ 2.00863696 -0.96943284 1.30917045]
+ État intermédiaire en itération courante : [ 2.00453385 -0.94011716 1.51021882]
+ État intermédiaire en itération courante : [ 2.00013539 -0.93313893 1.80539433]
+ État intermédiaire en itération courante : [ 1.95437244 -0.76890418 2.04566856]
+ État intermédiaire en itération courante : [ 1.99797362 -0.92538074 1.81674451]
+ État intermédiaire en itération courante : [ 1.99760514 -0.95929669 2.01402091]
+ État intermédiaire en itération courante : [ 1.99917565 -0.99152672 2.03171791]
+ État intermédiaire en itération courante : [ 1.99990376 -0.99963123 2.00671578]
+ État intermédiaire en itération courante : [ 1.99999841 -1.00005285 2.00039699]
+ État intermédiaire en itération courante : [ 2.00000014 -1.00000307 2.00000221]
+ État intermédiaire en itération courante : [ 2. -0.99999992 1.99999987]
+
+Calage de 3 coefficients pour une forme quadratique 1D sur 5 mesures
+--------------------------------------------------------------------
+
+Vecteur d'observation.............: [ 57. 2. 3. 17. 192.]
+État d'ébauche a priori...........: [1. 1. 1.]
+
+Coefficients théoriques attendus..: [ 2 -1 2]
+
+Nombre d'itérations...............: 25
+Nombre de simulations.............: 100
+Coefficients résultants du calage.: [ 2. -0.99999992 1.99999987]
--- /dev/null
+.. index:: single: 3DVAR (exemple)
+
+Premier exemple
+...............
+
+Cet exemple décrit le **recalage des paramètres** :math:`\mathbf{x}` d'un
+modèle d'observation :math:`H` quadratique. Ce modèle :math:`H` est représenté
+ici comme une fonction nommée ``QuadFunction`` pour les besoins de l'exemple.
+Cette fonction accepte en entrée le vecteur de coefficients :math:`\mathbf{x}`
+de la forme quadratique, et fournit en sortie le vecteur :math:`\mathbf{y}`
+d'évaluation du modèle quadratique aux points de contrôle internes, prédéfinis
+de manière statique dans le modèle.
+
+Le recalage s'effectue sur la base d'un jeu initial de coefficients (état
+d'ébauche désigné par ``Xb`` dans l'exemple), et avec l'information
+:math:`\mathbf{y}^o` (désignée par ``Yobs`` dans l'exemple) de 5 mesures
+obtenues à aux points de contrôle internes. On se place ici en expériences
+jumelles (voir :ref:`section_methodology_twin`), et les mesures sont
+considérées comme parfaites. On choisit de privilégier les observations, au
+détriment de l'ébauche, par l'indication artificielle d'une très importante
+variance d'erreur d'ébauche, ici de :math:`10^{6}`.
+
+L'ajustement s'effectue en affichant des résultats intermédiaires lors de
+l'optimisation itérative grâce à des "*observer*" (pour mémoire, voir
+:ref:`section_advanced_observer`).
--- /dev/null
+Sur la figure, les traits permettent simplement de joindre les valeurs simulées
+aux points de contrôle et d'en faciliter la lisibilité.
+
+On constate, dans ce cas particulièrement simple, que l'optimum est identique
+aux valeurs de coefficients ``[2, -1, 2]`` utilisés pour générer les
+pseudo-observations de l'expérience jumelle. Naturellement, la simulation
+finale obtenue avec les coefficients recalés coïncide donc avec les
+observations aux 5 points de contrôle.
+
+Pour illustrer plus clairement le recalage effectué, en utilisant les mêmes
+informations, on peut tracer la version continue des courbes quadratiques dont
+on dispose pour l'état d'ébauche des coefficients, et pour l'état optimal des
+coefficients, en même temps que les valeurs obtenues aux points de contrôle. De plus,
+pour améliorer la lisibilité, on a retiré les traits qui joignent les
+simulations dans le graphique précédent.
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+from numpy import array, random
+random.seed(1234567)
+Xtrue = -0.37727
+Yobs = []
+for i in range(51):
+ Yobs.append([random.normal(Xtrue, 0.1, size=(1,)),])
+#
+print("Estimation variationnelle d'une trajectoire temporelle")
+print("------------------------------------------------------")
+print(" Observations bruitées acquises sur %i pas de temps"%(len(Yobs)-1,))
+print("")
+from adao import adaoBuilder
+case = adaoBuilder.New()
+#
+case.setBackground (Vector = [0.])
+case.setBackgroundError (ScalarSparseMatrix = 0.1**2)
+#
+case.setObservationOperator(Matrix = [1.])
+case.setObservation (VectorSerie = Yobs)
+case.setObservationError (ScalarSparseMatrix = 0.3**2)
+#
+case.setEvolutionModel (Matrix = [1.])
+case.setEvolutionError (ScalarSparseMatrix = 1e-5)
+#
+case.setAlgorithmParameters(
+ Algorithm="3DVAR",
+ Parameters={
+ "EstimationOf":"State",
+ "StoreSupplementaryCalculations":[
+ "Analysis",
+ "APosterioriCovariance",
+ ],
+ },
+ )
+#
+case.execute()
+Xa = case.get("Analysis")
+Pa = case.get("APosterioriCovariance")
+#
+print(" État analysé à l'observation finale :", Xa[-1])
+print("")
+print(" Variance a posteriori finale :", Pa[-1])
+print("")
+#
+#-------------------------------------------------------------------------------
+#
+Observations = array([yo[0] for yo in Yobs])
+Estimates = array([xa[0] for xa in case.get("Analysis")])
+Variances = array([pa[0,0] for pa in case.get("APosterioriCovariance")])
+#
+import matplotlib.pyplot as plt
+plt.rcParams['figure.figsize'] = (10, 4)
+#
+plt.figure()
+plt.plot(Observations,'kx',label='Mesures bruitées')
+plt.plot(Estimates,'r-',label='État estimé')
+plt.axhline(Xtrue,color='b',label='Valeur vraie')
+plt.legend()
+plt.title('Estimation de l\'état', fontweight='bold')
+plt.xlabel('Pas d\'observation')
+plt.ylabel('Tension')
+plt.savefig("simple_3DVAR2_state.png")
+#
+plt.figure()
+iobs = range(1,len(Observations))
+plt.plot(iobs,Variances[iobs],label='Variance d\'erreur a posteriori')
+plt.title('Estimation de la variance d\'erreur a posteriori', fontweight='bold')
+plt.xlabel('Pas d\'observation')
+plt.ylabel('$(Tension)^2$')
+plt.setp(plt.gca(),'ylim',[0,.01])
+plt.savefig("simple_3DVAR2_variance.png")
--- /dev/null
+Estimation variationnelle d'une trajectoire temporelle
+------------------------------------------------------
+ Observations bruitées acquises sur 50 pas de temps
+
+ État analysé à l'observation finale : [-0.37110687]
+
+ Variance a posteriori finale : [[0.009]]
+
--- /dev/null
+Deuxième exemple
+................
+
+Le 3DVAR peut aussi être utilisé pour une **analyse temporelle des observations
+d'un modèle dynamique donné**. Dans ce cas, l'analyse est conduite de manière
+itérative, lors de l'arrivée de chaque observation. On utilise pour cet exemple
+le même système dynamique simple [Welch06]_ que celui qui est analysé dans les
+:ref:`section_ref_algorithm_KalmanFilter_examples` du Filtre de Kalman.
+
+A chaque étape, l'analyse 3DVAR classique remet à jour uniquement l'état du
+système. Moyennant une modification des valeurs de covariances *a priori* par
+rapport aux hypothèses initiales du filtrage, cette réanalyse 3DVAR permet de
+converger vers la trajectoire vraie, comme l'illustre la figure associée, de
+manière ici un peu plus lente qu'avec un Filtre de Kalman.
+
+.. note::
+
+ Remarque concernant les covariances *a posteriori* : classiquement,
+ l'analyse itérative 3DVAR remet à jour uniquement l'état et non pas sa
+ covariance. Comme les hypothèses d'opérateurs et de covariance *a priori*
+ restent inchangées ici au court de l'évolution, la covariance *a
+ posteriori* est constante. Le tracé qui suit de cette covariance *a
+ posteriori* permet d'insister sur cette propriété tout à fait attendue de
+ l'analyse 3DVAR. Une hypothèse plus évoluée est proposée dans l'exemple qui
+ suit.
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+from numpy import array, random
+random.seed(1234567)
+Xtrue = -0.37727
+Yobs = []
+for i in range(51):
+ Yobs.append([random.normal(Xtrue, 0.1, size=(1,)),])
+#
+print("Estimation variationnelle d'une trajectoire temporelle")
+print("------------------------------------------------------")
+print(" Observations bruitées acquises sur %i pas de temps"%(len(Yobs)-1,))
+print("")
+from adao import adaoBuilder
+case = adaoBuilder.New()
+#
+case.setBackground (Vector = [0.])
+#
+case.setObservationOperator(Matrix = [1.])
+case.setObservationError (ScalarSparseMatrix = 0.3**2)
+#
+case.setEvolutionModel (Matrix = [1.])
+case.setEvolutionError (ScalarSparseMatrix = 1e-5)
+#
+case.setAlgorithmParameters(
+ Algorithm="3DVAR",
+ Parameters={
+ "EstimationOf":"State",
+ "StoreSupplementaryCalculations":[
+ "Analysis",
+ "APosterioriCovariance",
+ ],
+ },
+ )
+#
+# Boucle pour obtenir une analyse à l'arrivée de chaque observation
+#
+Ca = 1.
+for i in range(1,len(Yobs)):
+ case.setObservation(Vector = Yobs[i])
+ case.setBackgroundError(ScalarSparseMatrix = Ca) # Covariance remise à jour
+ case.execute( nextStep = True )
+ #
+ # Hypothèse de décroissance de la covariance a priori de l'erreur d'ébauche
+ #
+ if float(Ca) <= 0.1**2:
+ Ca = 0.1**2
+ else:
+ Ca = Ca * 0.9**2
+#
+Xa = case.get("Analysis")
+Pa = case.get("APosterioriCovariance")
+#
+print(" État analysé à l'observation finale :", Xa[-1])
+print("")
+print(" Variance a posteriori finale :", Pa[-1])
+print("")
+#
+#-------------------------------------------------------------------------------
+#
+Observations = array([yo[0] for yo in Yobs])
+Estimates = array([xa[0] for xa in case.get("Analysis")])
+Variances = array([pa[0,0] for pa in case.get("APosterioriCovariance")])
+#
+import matplotlib.pyplot as plt
+plt.rcParams['figure.figsize'] = (10, 4)
+#
+plt.figure()
+plt.plot(Observations,'kx',label='Mesures bruitées')
+plt.plot(Estimates,'r-',label='État estimé')
+plt.axhline(Xtrue,color='b',label='Valeur vraie')
+plt.legend()
+plt.title('Estimation de l\'état', fontweight='bold')
+plt.xlabel('Pas d\'observation')
+plt.ylabel('Tension')
+plt.savefig("simple_3DVAR3_state.png")
+#
+plt.figure()
+iobs = range(1,len(Observations))
+plt.plot(iobs,Variances[iobs],label='Variance d\'erreur a posteriori')
+plt.title('Estimation de la variance d\'erreur a posteriori', fontweight='bold')
+plt.xlabel('Pas d\'observation')
+plt.ylabel('$(Tension)^2$')
+plt.savefig("simple_3DVAR3_variance.png")
--- /dev/null
+Estimation variationnelle d'une trajectoire temporelle
+------------------------------------------------------
+ Observations bruitées acquises sur 50 pas de temps
+
+ État analysé à l'observation finale : [-0.37334336]
+
+ Variance a posteriori finale : [[0.009]]
+
--- /dev/null
+Troisième exemple
+.................
+
+A partir de l'exemple précédent, si l'on veut adapter la convergence temporelle
+du 3DVAR, on peut modifier par exemple les hypothèses de covariance *a priori*
+des erreurs d'ébauche au cours des itérations. Cette remise à jour est une
+**hypothèse** de l'utilisateur, et il y a de multiples possibilités qui vont
+dépendre de la physique du cas. On en illustre une ici.
+
+On choisit, arbitrairement, de faire décroître la covariance *a priori* des
+erreurs d'ébauche d'un facteur constant :math:`0.9^2=0.81` tant qu'elle reste
+supérieure à une valeur limite de :math:`0.1^2=0.01` (qui est la valeur fixe de
+covariance *a priori* des erreurs d'ébauche de l'exemple précédent), sachant
+qu'elle commence à la valeur `1` (qui est la valeur fixe de covariance *a
+priori* des erreurs d'ébauche utilisée pour le premier pas du filtrage de
+Kalman). Cette valeur est remise à jour à chaque pas, en la réinjectant comme
+covariance *a priori* de l'état qui est utilisé comme ébauche lors du pas
+suivant d'analyse, dans une boucle explicite.
+
+On constate dans ce cas que l'estimation d'état converge plus vite vers la
+valeur vraie, et que l'assimilation se comporte ensuite de manière similaire
+aux :ref:`section_ref_algorithm_KalmanFilter_examples` pour le Filtre de
+Kalman, ou à l'exemple précédent avec les covariances *a priori* manuellement
+adaptée. De plus, la covariance *a posteriori* décroît tant que l'on force la
+décroissance de la covariance *a priori*.
+
+.. note::
+
+ On insiste sur le fait que les variations de covariance *a priori*, qui
+ conditionnent les variations de covariance *a posteriori*, relèvent d'une
+ **hypothèse utilisateur** et non pas d'une obligation. Cette hypothèse doit
+ donc être **adaptée en fonction du cas physique**.
Exemples d'utilisation en Python (TUI)
++++++++++++++++++++++++++++++++++++++
-Voici un exemple très simple d'usage de l'algorithme proposé et de ses
-paramètres, écrit en :ref:`section_tui`, et dont les informations indiquées en
-entrée permettent de définir un cas équivalent en interface graphique.
+Voici un ou des exemples très simple d'usage de l'algorithme proposé et de ses
+paramètres, écrit en :ref:`section_tui`. De plus, les informations indiquées en
+entrée permettent aussi de définir un cas équivalent en interface graphique
+:ref:`section_gui_in_salome`.
Python, 3.6.5, 3.10.9
Numpy, 1.14.3, 1.24.2
- Scipy, 1.1.0, 1.10.0
+ Scipy, 0.19.1, 1.10.1
MatplotLib, 2.2.2, 3.7.0
GnuplotPy, 1.8, 1.8
NLopt, 2.4.2, 2.7.1
