src/tests/daComposant/Plateforme/test002_Blue.py

   1 #-*-coding:iso-8859-1-*-
   2 __doc__ = """
   3     Cas-test vérifiant que l'application des coefficients de correction so et sb
   4     conduit à des matrices R et B pour lesquelles ces coefficients sont unitaires.
   5 """
   6 __author__ = "Jean-Philippe ARGAUD - Mars 2008"
   7
   8 import numpy
   9 from daCore.AssimilationStudy import AssimilationStudy
  10
  11 #===============================================================================
  12 def test(precision = 1.e-13, dimension = 3):
  13     """
  14     Cas-test vérifiant que l'application des coefficients de correction so et sb
  15     conduit à des matrices R et B pour lesquelles ces coefficients sont unitaires.
  16     """
  17     #
  18     # Définition des données "théoriques" vraies
  19     # ------------------------------------------
  20     xt = numpy.matrix(numpy.arange(dimension)).T
  21     Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
  22     Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
  23     #
  24     H  = numpy.matrix(numpy.core.identity(dimension))
  25     #
  26     xb = xt + Eb
  27     yo = H * xt + Eo
  28     #
  29     xb = xb.A1
  30     yo = yo.A1
  31     #
  32     # Définition des matrices d'erreurs
  33     # ---------------------------------
  34     R  = numpy.matrix(numpy.core.identity(dimension)).T
  35     B  = numpy.matrix(numpy.core.identity(dimension)).T
  36     #
  37     # Analyse BLUE
  38     # ------------
  39     ADD = AssimilationStudy()
  40     ADD.setBackground         (asVector     = xb )
  41     ADD.setBackgroundError    (asCovariance = B )
  42     ADD.setObservation        (asVector     = yo )
  43     ADD.setObservationError   (asCovariance = R )
  44     ADD.setObservationOperator(asMatrix     = H )
  45     #
  46     ADD.setControls()
  47     ADD.setAlgorithm(choice="Blue")
  48     #
  49     ADD.analyze()
  50     #
  51     xa        = numpy.array(ADD.get("Analysis").valueserie(0))
  52     d         = numpy.array(ADD.get("Innovation").valueserie(0))
  53     SigmaObs2 = float( numpy.dot(d,(yo-numpy.dot(H,xa)).A1) / R.trace() )
  54     SigmaBck2 = float( numpy.dot(d,numpy.dot(H,(xa - xb)).A1) /(H * B * H.T).trace() )
  55     #
  56     # Analyse BLUE avec correction des matrices R et B
  57     # Attention : ce second calcul de BLUE avec le meme objet ADD
  58     #             conduit à stocker les résultats dans le second step,
  59     #             donc il faut appeller "valueserie(1)"
  60     # ------------------------------------------------
  61     ADD.setBackgroundError (asCovariance = SigmaBck2*B )
  62     ADD.setObservationError(asCovariance = SigmaObs2*R )
  63     ADD.analyze()
  64     new_xa        = numpy.array(ADD.get("Analysis").valueserie(1))
  65     new_d         = numpy.array(ADD.get("Innovation").valueserie(1))
  66     new_SigmaObs2 = float( numpy.dot(new_d,(yo-numpy.dot(H,new_xa)).A1) / (SigmaObs2*R.trace()) )
  67     new_SigmaBck2 = float( numpy.dot(new_d,numpy.dot(H,(new_xa - xb)).A1) /(H * (SigmaBck2*B) * H.T).trace() )
  68     #
  69     # Vérification du résultat
  70     # ------------------------
  71     if max(abs(xa - new_xa)) > precision:
  72         raise ValueError("Résultat du test erroné (1)")
  73     elif max(abs(d - new_d)) > precision:
  74         raise ValueError("Résultat du test erroné (2)")
  75     elif abs(new_SigmaObs2-1.) > precision:
  76         print "new_SigmaObs2 =",new_SigmaObs2
  77         raise ValueError("Résultat du test erroné (3)")
  78     elif abs(new_SigmaBck2-1.) > precision :
  79         print "new_SigmaBck2 =",new_SigmaBck2
  80         raise ValueError("Résultat du test erroné (4)")
  81     else:
  82         print test.__doc__
  83         print "    Test correct, erreur maximale inférieure à %s"%precision
  84         print
  85
  86 #===============================================================================
  87 if __name__ == "__main__":
  88
  89     print
  90     print "AUTODIAGNOSTIC"
  91     print "=============="
  92     # numpy.random.seed(1000)
  93
  94     test(dimension = 100)