src/tests/daComposant/Plateforme/test003_Blue.py

   1 #-*-coding:iso-8859-1-*-
   2 __doc__ = """
   3     Cas-test vérifiant que si l'erreur sur le background est nulle et que
   4     l'erreur sur les observations est connue, alors l'analyse donne le "milieu"
   5     du background et des observations.
   6 """
   7 __author__ = "Jean-Philippe ARGAUD - Mars 2008"
   8
   9 import numpy
  10 from daCore.AssimilationStudy import AssimilationStudy
  11
  12 #===============================================================================
  13 def test(precision = 1.e-13, dimension = 3):
  14     """
  15     Cas-test vérifiant que si l'erreur sur le background est nulle et que
  16     l'erreur sur les observations est connue, alors l'analyse donne le "milieu"
  17     du background et des observations.
  18     """
  19     #
  20     # Définition des données "théoriques" vraies
  21     # ------------------------------------------
  22     xt = numpy.matrix(numpy.arange(dimension)).T
  23     Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
  24     Eb = numpy.matrix(numpy.zeros((dimension,))).T
  25     #
  26     H  = numpy.matrix(numpy.core.identity(dimension))
  27     #
  28     xb = xt + Eb
  29     yo = H * xt + Eo
  30     #
  31     xb = xb.A1
  32     yo = yo.A1
  33     #
  34     # Définition des matrices de covariances d'erreurs
  35     # ------------------------------------------------
  36     R  = numpy.matrix(numpy.core.identity(dimension)).T
  37     B  = numpy.matrix(numpy.core.identity(dimension)).T
  38     #
  39     # Analyse BLUE
  40     # ------------
  41     ADD = AssimilationStudy()
  42     ADD.setBackground         (asVector     = xb )
  43     ADD.setBackgroundError    (asCovariance = B )
  44     ADD.setObservation        (asVector     = yo )
  45     ADD.setObservationError   (asCovariance = R )
  46     ADD.setObservationOperator(asMatrix     = H )
  47     #
  48     ADD.setControls()
  49     ADD.setAlgorithm(choice="Blue")
  50     #
  51     ADD.analyze()
  52     #
  53     Xa = ADD.get("Analysis")
  54     xa        = numpy.matrix(Xa.valueserie(0)).T
  55     SigmaObs2 = ADD.get("SigmaObs2")
  56     SigmaBck2 = ADD.get("SigmaBck2")
  57     d         = ADD.get("Innovation")
  58     #
  59     # Vérification du résultat
  60     # ------------------------
  61     if max(abs(xa.A1 - xb - Eo.A1/2.)) > precision:
  62         raise ValueError("Résultat du test erroné (1)")
  63     elif max(abs(yo - (H * xa).A1 - Eo.A1/2.)) > precision:
  64         raise ValueError("Résultat du test erroné (2)")
  65     else:
  66         print test.__doc__
  67         print "    Test correct, erreur maximale inférieure à %s"%precision
  68         print
  69
  70 #===============================================================================
  71 if __name__ == "__main__":
  72
  73     print
  74     print "AUTODIAGNOSTIC"
  75     print "=============="
  76     # numpy.random.seed(1000)
  77
  78     test(dimension = 100)