From: Jean-Philippe ARGAUD Date: Tue, 28 Feb 2023 12:51:43 +0000 (+0100) Subject: Documentation update et example improvement X-Git-Tag: V9_11_0a1~2 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=6f9fad9ddee0ebcd935ae1652ade66d34accdb65;p=modules%2Fadao.git Documentation update et example improvement --- diff --git a/doc/en/ref_algorithm_3DVAR.rst b/doc/en/ref_algorithm_3DVAR.rst index 08c4b77..b2023f1 100644 --- a/doc/en/ref_algorithm_3DVAR.rst +++ b/doc/en/ref_algorithm_3DVAR.rst @@ -255,18 +255,71 @@ StoreSupplementaryCalculations .. 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% @@ -275,6 +328,7 @@ StoreSupplementaryCalculations - :ref:`section_ref_algorithm_Blue` - :ref:`section_ref_algorithm_ExtendedBlue` +- :ref:`section_ref_algorithm_KalmanFilter` - :ref:`section_ref_algorithm_LinearityTest` .. ------------------------------------ .. diff --git a/doc/en/scripts/simple_3DVAR.png b/doc/en/scripts/simple_3DVAR.png deleted file mode 100644 index 6bb44ce..0000000 Binary files a/doc/en/scripts/simple_3DVAR.png and /dev/null differ diff --git a/doc/en/scripts/simple_3DVAR.py b/doc/en/scripts/simple_3DVAR.py deleted file mode 100644 index f207887..0000000 --- a/doc/en/scripts/simple_3DVAR.py +++ /dev/null @@ -1,74 +0,0 @@ -# -*- 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") diff --git a/doc/en/scripts/simple_3DVAR.res b/doc/en/scripts/simple_3DVAR.res deleted file mode 100644 index 8713620..0000000 --- a/doc/en/scripts/simple_3DVAR.res +++ /dev/null @@ -1,40 +0,0 @@ -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] diff --git a/doc/en/scripts/simple_3DVAR.rst b/doc/en/scripts/simple_3DVAR.rst deleted file mode 100644 index cb47c74..0000000 --- a/doc/en/scripts/simple_3DVAR.rst +++ /dev/null @@ -1,17 +0,0 @@ -.. 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. diff --git a/doc/en/scripts/simple_3DVAR1.png b/doc/en/scripts/simple_3DVAR1.png new file mode 100644 index 0000000..987dcaa Binary files /dev/null and b/doc/en/scripts/simple_3DVAR1.png differ diff --git a/doc/en/scripts/simple_3DVAR1.py b/doc/en/scripts/simple_3DVAR1.py new file mode 100644 index 0000000..7be90ad --- /dev/null +++ b/doc/en/scripts/simple_3DVAR1.py @@ -0,0 +1,74 @@ +# -*- 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") diff --git a/doc/en/scripts/simple_3DVAR1.res b/doc/en/scripts/simple_3DVAR1.res new file mode 100644 index 0000000..8713620 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR1.res @@ -0,0 +1,40 @@ +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] diff --git a/doc/en/scripts/simple_3DVAR1.rst b/doc/en/scripts/simple_3DVAR1.rst new file mode 100644 index 0000000..3a207f6 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR1.rst @@ -0,0 +1,25 @@ +.. 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. diff --git a/doc/en/scripts/simple_3DVAR1Plus.png b/doc/en/scripts/simple_3DVAR1Plus.png new file mode 100644 index 0000000..cb3ac04 Binary files /dev/null and b/doc/en/scripts/simple_3DVAR1Plus.png differ diff --git a/doc/en/scripts/simple_3DVAR1Plus.rst b/doc/en/scripts/simple_3DVAR1Plus.rst new file mode 100644 index 0000000..77f4191 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR1Plus.rst @@ -0,0 +1,14 @@ +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. diff --git a/doc/en/scripts/simple_3DVAR2.py b/doc/en/scripts/simple_3DVAR2.py new file mode 100644 index 0000000..a7e8db0 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR2.py @@ -0,0 +1,73 @@ +# -*- 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") diff --git a/doc/en/scripts/simple_3DVAR2.res b/doc/en/scripts/simple_3DVAR2.res new file mode 100644 index 0000000..ec9bd02 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR2.res @@ -0,0 +1,8 @@ +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]] + diff --git a/doc/en/scripts/simple_3DVAR2.rst b/doc/en/scripts/simple_3DVAR2.rst new file mode 100644 index 0000000..56ce659 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR2.rst @@ -0,0 +1,24 @@ +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. diff --git a/doc/en/scripts/simple_3DVAR2_state.png b/doc/en/scripts/simple_3DVAR2_state.png new file mode 100644 index 0000000..0f0f588 Binary files /dev/null and b/doc/en/scripts/simple_3DVAR2_state.png differ diff --git a/doc/en/scripts/simple_3DVAR2_variance.png b/doc/en/scripts/simple_3DVAR2_variance.png new file mode 100644 index 0000000..943e5e1 Binary files /dev/null and b/doc/en/scripts/simple_3DVAR2_variance.png differ diff --git a/doc/en/scripts/simple_3DVAR3.py b/doc/en/scripts/simple_3DVAR3.py new file mode 100644 index 0000000..2d9ad44 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR3.py @@ -0,0 +1,84 @@ +# -*- 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") diff --git a/doc/en/scripts/simple_3DVAR3.res b/doc/en/scripts/simple_3DVAR3.res new file mode 100644 index 0000000..21c9b65 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR3.res @@ -0,0 +1,8 @@ +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]] + diff --git a/doc/en/scripts/simple_3DVAR3.rst b/doc/en/scripts/simple_3DVAR3.rst new file mode 100644 index 0000000..0d6e9a5 --- /dev/null +++ b/doc/en/scripts/simple_3DVAR3.rst @@ -0,0 +1,31 @@ +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**. diff --git a/doc/en/scripts/simple_3DVAR3_state.png b/doc/en/scripts/simple_3DVAR3_state.png new file mode 100644 index 0000000..9f4ccc3 Binary files /dev/null and b/doc/en/scripts/simple_3DVAR3_state.png differ diff --git a/doc/en/scripts/simple_3DVAR3_variance.png b/doc/en/scripts/simple_3DVAR3_variance.png new file mode 100644 index 0000000..3c3bec3 Binary files /dev/null and b/doc/en/scripts/simple_3DVAR3_variance.png differ diff --git a/doc/en/snippets/Header2Algo09.rst b/doc/en/snippets/Header2Algo09.rst index 0024c9d..0119958 100644 --- a/doc/en/snippets/Header2Algo09.rst +++ b/doc/en/snippets/Header2Algo09.rst @@ -1,6 +1,6 @@ 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`. diff --git a/doc/en/snippets/ModuleCompatibility.rst b/doc/en/snippets/ModuleCompatibility.rst index 1138bf9..e5a5fcd 100644 --- a/doc/en/snippets/ModuleCompatibility.rst +++ b/doc/en/snippets/ModuleCompatibility.rst @@ -14,7 +14,7 @@ guarantee). 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 diff --git a/doc/fr/ref_algorithm_3DVAR.rst b/doc/fr/ref_algorithm_3DVAR.rst index d3c10b4..c8cf19d 100644 --- a/doc/fr/ref_algorithm_3DVAR.rst +++ b/doc/fr/ref_algorithm_3DVAR.rst @@ -260,19 +260,71 @@ StoreSupplementaryCalculations .. 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% @@ -281,6 +333,7 @@ StoreSupplementaryCalculations - :ref:`section_ref_algorithm_Blue` - :ref:`section_ref_algorithm_ExtendedBlue` +- :ref:`section_ref_algorithm_KalmanFilter` - :ref:`section_ref_algorithm_LinearityTest` .. ------------------------------------ .. diff --git a/doc/fr/scripts/simple_3DVAR.png b/doc/fr/scripts/simple_3DVAR.png deleted file mode 100644 index 00388ba..0000000 Binary files a/doc/fr/scripts/simple_3DVAR.png and /dev/null differ diff --git a/doc/fr/scripts/simple_3DVAR.py b/doc/fr/scripts/simple_3DVAR.py deleted file mode 100644 index 88c307c..0000000 --- a/doc/fr/scripts/simple_3DVAR.py +++ /dev/null @@ -1,74 +0,0 @@ -# -*- 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") diff --git a/doc/fr/scripts/simple_3DVAR.res b/doc/fr/scripts/simple_3DVAR.res deleted file mode 100644 index d28148a..0000000 --- a/doc/fr/scripts/simple_3DVAR.res +++ /dev/null @@ -1,40 +0,0 @@ -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] diff --git a/doc/fr/scripts/simple_3DVAR.rst b/doc/fr/scripts/simple_3DVAR.rst deleted file mode 100644 index c2c2c29..0000000 --- a/doc/fr/scripts/simple_3DVAR.rst +++ /dev/null @@ -1,18 +0,0 @@ -.. 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. diff --git a/doc/fr/scripts/simple_3DVAR1.png b/doc/fr/scripts/simple_3DVAR1.png new file mode 100644 index 0000000..9db9d39 Binary files /dev/null and b/doc/fr/scripts/simple_3DVAR1.png differ diff --git a/doc/fr/scripts/simple_3DVAR1.py b/doc/fr/scripts/simple_3DVAR1.py new file mode 100644 index 0000000..12da2a8 --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR1.py @@ -0,0 +1,74 @@ +# -*- 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") diff --git a/doc/fr/scripts/simple_3DVAR1.res b/doc/fr/scripts/simple_3DVAR1.res new file mode 100644 index 0000000..d28148a --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR1.res @@ -0,0 +1,40 @@ +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] diff --git a/doc/fr/scripts/simple_3DVAR1.rst b/doc/fr/scripts/simple_3DVAR1.rst new file mode 100644 index 0000000..c1295d7 --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR1.rst @@ -0,0 +1,25 @@ +.. 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`). diff --git a/doc/fr/scripts/simple_3DVAR1Plus.png b/doc/fr/scripts/simple_3DVAR1Plus.png new file mode 100644 index 0000000..78b62f0 Binary files /dev/null and b/doc/fr/scripts/simple_3DVAR1Plus.png differ diff --git a/doc/fr/scripts/simple_3DVAR1Plus.rst b/doc/fr/scripts/simple_3DVAR1Plus.rst new file mode 100644 index 0000000..c655224 --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR1Plus.rst @@ -0,0 +1,15 @@ +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. diff --git a/doc/fr/scripts/simple_3DVAR2.py b/doc/fr/scripts/simple_3DVAR2.py new file mode 100644 index 0000000..171d481 --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR2.py @@ -0,0 +1,73 @@ +# -*- 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") diff --git a/doc/fr/scripts/simple_3DVAR2.res b/doc/fr/scripts/simple_3DVAR2.res new file mode 100644 index 0000000..e0147aa --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR2.res @@ -0,0 +1,8 @@ +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]] + diff --git a/doc/fr/scripts/simple_3DVAR2.rst b/doc/fr/scripts/simple_3DVAR2.rst new file mode 100644 index 0000000..4e07966 --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR2.rst @@ -0,0 +1,25 @@ +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. diff --git a/doc/fr/scripts/simple_3DVAR2_state.png b/doc/fr/scripts/simple_3DVAR2_state.png new file mode 100644 index 0000000..512b8ab Binary files /dev/null and b/doc/fr/scripts/simple_3DVAR2_state.png differ diff --git a/doc/fr/scripts/simple_3DVAR2_variance.png b/doc/fr/scripts/simple_3DVAR2_variance.png new file mode 100644 index 0000000..91ff438 Binary files /dev/null and b/doc/fr/scripts/simple_3DVAR2_variance.png differ diff --git a/doc/fr/scripts/simple_3DVAR3.py b/doc/fr/scripts/simple_3DVAR3.py new file mode 100644 index 0000000..9d18eeb --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR3.py @@ -0,0 +1,84 @@ +# -*- 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") diff --git a/doc/fr/scripts/simple_3DVAR3.res b/doc/fr/scripts/simple_3DVAR3.res new file mode 100644 index 0000000..e7f8f2e --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR3.res @@ -0,0 +1,8 @@ +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]] + diff --git a/doc/fr/scripts/simple_3DVAR3.rst b/doc/fr/scripts/simple_3DVAR3.rst new file mode 100644 index 0000000..deabaea --- /dev/null +++ b/doc/fr/scripts/simple_3DVAR3.rst @@ -0,0 +1,32 @@ +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**. diff --git a/doc/fr/scripts/simple_3DVAR3_state.png b/doc/fr/scripts/simple_3DVAR3_state.png new file mode 100644 index 0000000..b8dd2a0 Binary files /dev/null and b/doc/fr/scripts/simple_3DVAR3_state.png differ diff --git a/doc/fr/scripts/simple_3DVAR3_variance.png b/doc/fr/scripts/simple_3DVAR3_variance.png new file mode 100644 index 0000000..d74815d Binary files /dev/null and b/doc/fr/scripts/simple_3DVAR3_variance.png differ diff --git a/doc/fr/snippets/Header2Algo09.rst b/doc/fr/snippets/Header2Algo09.rst index c2a017d..3e7ea72 100644 --- a/doc/fr/snippets/Header2Algo09.rst +++ b/doc/fr/snippets/Header2Algo09.rst @@ -1,6 +1,7 @@ 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`. diff --git a/doc/fr/snippets/ModuleCompatibility.rst b/doc/fr/snippets/ModuleCompatibility.rst index b9adf27..2474ca6 100644 --- a/doc/fr/snippets/ModuleCompatibility.rst +++ b/doc/fr/snippets/ModuleCompatibility.rst @@ -14,7 +14,7 @@ la version atteinte (mais cela reste sans garantie). 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