src/tests/daComposant/Plateforme/test008_Kalman.py

   1 #-*-coding:iso-8859-1-*-
   2 __doc__ = """
   3     Cas-test vérifiant le fonctionnement du filtre de Kalman sur un système
   4     dynamique de trajectoire 1D constante
   5 """
   6 __author__ = "Jean-Philippe ARGAUD - Septembre 2008"
   7
   8 import numpy
   9 from daCore.AssimilationStudy import AssimilationStudy
  10 from daCore.Persistence       import OneScalar
  11
  12 #===============================================================================
  13 def test(dimension = 3):
  14     """
  15     Cas-test vérifiant le fonctionnement du filtre de Kalman sur un système
  16     dynamique de trajectoire 1D constante
  17     """
  18     print test.__doc__
  19     #
  20     # Définition des données
  21     # ----------------------
  22     a_size = (dimension,)
  23     #
  24     # Valeur vraie
  25     xt = -0.4
  26     Xt = OneScalar("Valeur vraie", basetype=float)
  27     Xt.store(xt)
  28     for i in range(dimension):
  29         Xt.store(xt)
  30     #
  31     # Observations bruitées
  32     yo = numpy.random.normal(xt, 0.1, size=a_size)
  33     Yo = OneScalar("Observations", basetype=float)
  34     Yo.store(0.)
  35     for v in yo:
  36         Yo.store(v)
  37     #
  38     # Création de l'étude et résolution
  39     # ---------------------------------
  40     ADD = AssimilationStudy("Assimilation temporelle de Kalman")
  41     #
  42     ADD.setBackground         (asVector     = "0.")
  43     ADD.setBackgroundError    (asCovariance = "1.")
  44     #
  45     ADD.setObservationOperator(asMatrix     = "1.")
  46     ADD.setObservation        (asPersistentVector = Yo)
  47     ADD.setObservationError   (asCovariance = "1.e-2")
  48     #
  49     ADD.setEvolutionModel     (asMatrix     = "1")
  50     ADD.setEvolutionError     (asCovariance = "1.e-5")
  51     #
  52     ADD.setControls()
  53     ADD.setAlgorithm(choice="Kalman")
  54     #
  55     ADD.analyze()
  56     #
  57     Xa = ADD.get("Analysis")
  58     print "    Valeur vraie visée........................:",xt
  59     print "    Ebauche, i.e. valeur initiale d'analyse...:",Xa.valueserie(0)[0]
  60     print "    Nombre d'analyses (sans l'ébauche)........:",Xa.stepnumber()-1
  61     print "    Moyenne des analyses......................:",Xa.stepmean()
  62     #
  63     # Biais des erreurs
  64     EpsY = []
  65     for i in range(Yo.stepnumber()):
  66         EpsY.append(Yo.valueserie(i) - Xt.valueserie(i))
  67     print "    Biais des erreurs <Obs-Vraie>.............:",numpy.array(EpsY).mean()
  68     print "    Variance des erreurs <Obs-Vraie>..........:",numpy.array(EpsY).var()
  69     EpsY = []
  70     for i in range(Xa.stepnumber()):
  71         EpsY.append(Xa.valueserie(i)[0] - Xt.valueserie(i))
  72     print "    Biais des erreurs <Ana-Vraie>.............:",numpy.array(EpsY).mean()
  73     print "    Variance des erreurs <Ana-Vraie>..........:",numpy.array(EpsY).var()
  74     print
  75     #
  76     ADD.setDiagnostic("PlotVectors", "Affichage de Xa et Xt")
  77     MonPlot = ADD.get("Affichage de Xa et Xt")
  78     MonPlot.calculate(
  79         ( [ x[0] for x in Xa.valueserie()], Xt.valueserie(), Yo.valueserie() ),
  80         title = "Analyse de Kalman sur trajectoire constante",
  81         ltitle = ["Analyse", "Valeur vraie", "Observations"],
  82         filename = "kalman_sur_trajectoire_constante.ps",
  83         pause = False,
  84         )
  85
  86 #===============================================================================
  87 if __name__ == "__main__":
  88
  89     print
  90     print "AUTODIAGNOSTIC"
  91     print "=============="
  92     numpy.random.seed(1000)
  93
  94     test(100)