doc/fr/scripts/simple_3DVAR.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 from numpy import array, ravel
   4 def QuadFunction( coefficients ):
   5     """
   6     Simulation quadratique aux points x : y = a x^2 + b x + c
   7     """
   8     a, b, c = list(ravel(coefficients))
   9     x_points = (-5, 0, 1, 3, 10)
  10     y_points = []
  11     for x in x_points:
  12         y_points.append( a*x*x + b*x + c )
  13     return array(y_points)
  14 #
  15 Xb   = array([1., 1., 1.])
  16 Yobs = array([57, 2, 3, 17, 192])
  17 #
  18 print("Résolution du problème de calage")
  19 print("--------------------------------")
  20 print("")
  21 from adao import adaoBuilder
  22 case = adaoBuilder.New()
  23 case.setBackground( Vector = Xb, Stored=True )
  24 case.setBackgroundError( ScalarSparseMatrix = 1.e6 )
  25 case.setObservation( Vector = Yobs, Stored=True )
  26 case.setObservationError( ScalarSparseMatrix = 1. )
  27 case.setObservationOperator( OneFunction = QuadFunction )
  28 case.setAlgorithmParameters(
  29     Algorithm='3DVAR',
  30     Parameters={
  31         'MaximumNumberOfIterations': 100,
  32         'StoreSupplementaryCalculations': [
  33             'CurrentState',
  34             ],
  35         },
  36     )
  37 case.setObserver(
  38     Info="  État intermédiaire en itération courante :",
  39     Template='ValuePrinter',
  40     Variable='CurrentState',
  41     )
  42 case.execute()
  43 print("")
  44 #
  45 #-------------------------------------------------------------------------------
  46 #
  47 print("Calage de %i coefficients pour une forme quadratique 1D sur %i mesures"%(
  48     len(case.get('Background')),
  49     len(case.get('Observation')),
  50     ))
  51 print("--------------------------------------------------------------------")
  52 print("")
  53 print("Vecteur d'observation.............:", ravel(case.get('Observation')))
  54 print("État d'ébauche a priori...........:", ravel(case.get('Background')))
  55 print("")
  56 print("Coefficients théoriques attendus..:", ravel((2,-1,2)))
  57 print("")
  58 print("Nombre d'itérations...............:", len(case.get('CurrentState')))
  59 print("Nombre de simulations.............:", len(case.get('CurrentState'))*4)
  60 print("Coefficients résultants du calage.:", ravel(case.get('Analysis')[-1]))
  61 #
  62 Xa = case.get('Analysis')[-1]
  63 import matplotlib.pyplot as plt
  64 plt.rcParams['figure.figsize'] = (10, 4)
  65 #
  66 plt.figure()
  67 plt.plot((-5,0,1,3,10),QuadFunction(Xb),'b-',label="Simulation à l'ébauche")
  68 plt.plot((-5,0,1,3,10),Yobs,            'kX',label='Observation',markersize=10)
  69 plt.plot((-5,0,1,3,10),QuadFunction(Xa),'r-',label="Simulation à l'optimum")
  70 plt.legend()
  71 plt.title('Calage de coefficients', fontweight='bold')
  72 plt.xlabel('Coordonnée arbitraire')
  73 plt.ylabel('Observations')
  74 plt.savefig("simple_3DVAR.png")