src/tests/daComposant/Plateforme/test004_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'on rajoute l'évaluation de l'opérateur
  16     d'observation au background, on obtient la même valeur que pour le BLUE
  17     normal.
  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     Hxb = H*xb
  31     #
  32     xb = xb.A1
  33     yo = yo.A1
  34     HXb = Hxb.A1
  35     #
  36     # Définition des matrices de covariances d'erreurs
  37     # ------------------------------------------------
  38     R  = numpy.matrix(numpy.core.identity(dimension))
  39     B  = numpy.matrix(numpy.core.identity(dimension))
  40     #
  41     # Analyse BLUE
  42     # ------------
  43     ADD = AssimilationStudy()
  44     ADD.setBackground         (asVector     = xb )
  45     ADD.setBackgroundError    (asCovariance = B )
  46     ADD.setObservation        (asVector     = yo )
  47     ADD.setObservationError   (asCovariance = R )
  48     ADD.setObservationOperator(asMatrix     = H )
  49     #
  50     ADD.setControls()
  51     ADD.setAlgorithm(choice="Blue")
  52     #
  53     ADD.analyze()
  54     #
  55     xa        = numpy.array(ADD.get("Analysis").valueserie(0))
  56     d         = numpy.array(ADD.get("Innovation").valueserie(0))
  57     SigmaObs2 = float( numpy.dot(d,(yo-numpy.dot(H,xa)).A1) / R.trace() )
  58     SigmaBck2 = float( numpy.dot(d,numpy.dot(H,(xa - xb)).A1) /(H * B * H.T).trace() )
  59     #
  60     # Analyse BLUE avec une évaluation au point Xb
  61     # Attention : ce second calcul de BLUE avec le meme objet ADD
  62     #             conduit à stocker les résultats dans le second step,
  63     #             donc il faut appeller "valueserie(1)"
  64     # ------------------------------------------------
  65     ADD.setObservationOperator(asMatrix     = H, appliedToX = {"HXb":HXb} )
  66     ADD.analyze()
  67     #
  68     new_xa        = numpy.array(ADD.get("Analysis").valueserie(1))
  69     new_d         = numpy.array(ADD.get("Innovation").valueserie(1))
  70     new_SigmaObs2 = float( numpy.dot(new_d,(yo-numpy.dot(H,new_xa)).A1) / R.trace() )
  71     new_SigmaBck2 = float( numpy.dot(new_d,numpy.dot(H,(new_xa - xb)).A1) /(H * B * H.T).trace() )
  72     #
  73     # Vérification du résultat
  74     # ------------------------
  75     if max(abs(xa - new_xa)) > precision:
  76         raise ValueError("Résultat du test erroné (1)")
  77     elif max(abs(d - new_d)) > precision:
  78         raise ValueError("Résultat du test erroné (2)")
  79     elif abs((new_SigmaObs2-SigmaObs2)/SigmaObs2) > precision:
  80         raise ValueError("Résultat du test erroné (3)")
  81     elif abs((new_SigmaBck2-SigmaBck2)/SigmaBck2) > precision :
  82         raise ValueError("Résultat du test erroné (4)")
  83     else:
  84         print test.__doc__
  85         print "    Test correct, erreur maximale inférieure à %s"%precision
  86         print
  87
  88 #===============================================================================
  89 if __name__ == "__main__":
  90
  91     print
  92     print "AUTODIAGNOSTIC"
  93     print "=============="
  94     # numpy.random.seed(1000)
  95
  96     test(dimension = 100)