test/test6901/Verification_des_Assimilation_Algorithms.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2019 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22 "Verification de la disponibilite de l'ensemble des algorithmes"
  23
  24 import sys
  25 import unittest
  26 import numpy
  27 from adao import adaoBuilder
  28
  29 # ==============================================================================
  30 class Test_Adao(unittest.TestCase):
  31     def test1(self):
  32         """Verification de la disponibilite de l'ensemble des algorithmes\n(Utilisation d'un operateur matriciel)"""
  33         print(self.test1.__doc__.strip()+"\n")
  34         Xa = {}
  35         for algo in ("3DVAR", "Blue", "ExtendedBlue", "LinearLeastSquares", "NonLinearLeastSquares", "DerivativeFreeOptimization"):
  36             print("")
  37             msg = "Algorithme en test : %s"%algo
  38             print(msg+"\n"+"-"*len(msg))
  39             #
  40             adaopy = adaoBuilder.New()
  41             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "Bounds":[[-1,10.],[-1,10.],[-1,10.]]})
  42             adaopy.setBackground         (Vector = [0,1,2])
  43             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
  44             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
  45             adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
  46             adaopy.setObservationOperator(Matrix = "1 0 0;0 2 0;0 0 3")
  47             adaopy.setObserver("Analysis",Template="ValuePrinter")
  48             adaopy.execute()
  49             Xa[algo] = adaopy.get("Analysis")[-1]
  50             del adaopy
  51         #
  52         for algo in ("ExtendedKalmanFilter", "KalmanFilter", "UnscentedKalmanFilter", "EnsembleKalmanFilter", "4DVAR"):
  53             print("")
  54             msg = "Algorithme en test : %s"%algo
  55             print(msg+"\n"+"-"*len(msg))
  56             #
  57             adaopy = adaoBuilder.New()
  58             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "SetSeed":1000})
  59             adaopy.setBackground         (Vector = [0,1,2])
  60             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
  61             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
  62             adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
  63             adaopy.setObservationOperator(Matrix = "1 0 0;0 2 0;0 0 3")
  64             adaopy.setEvolutionError     (ScalarSparseMatrix = 1.)
  65             adaopy.setEvolutionModel     (Matrix = "1 0 0;0 1 0;0 0 1")
  66             adaopy.setObserver("Analysis",Template="ValuePrinter")
  67             adaopy.execute()
  68             Xa[algo] = adaopy.get("Analysis")[-1]
  69             del adaopy
  70         #
  71         for algo in ("ParticleSwarmOptimization", "QuantileRegression", ):
  72             print("")
  73             msg = "Algorithme en test : %s"%algo
  74             print(msg+"\n"+"-"*len(msg))
  75             #
  76             adaopy = adaoBuilder.New()
  77             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"BoxBounds":3*[[-1,3]], "SetSeed":1000})
  78             adaopy.setBackground         (Vector = [0,1,2])
  79             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
  80             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
  81             adaopy.setObservationError   (DiagonalSparseMatrix = "1 2 3")
  82             adaopy.setObservationOperator(Matrix = "1 0 0;0 1 0;0 0 1")
  83             adaopy.setObserver("Analysis",Template="ValuePrinter")
  84             adaopy.execute()
  85             Xa[algo] = adaopy.get("Analysis")[-1]
  86             del adaopy
  87         #
  88         for algo in ("EnsembleBlue", ):
  89             print("")
  90             msg = "Algorithme en test : %s"%algo
  91             print(msg+"\n"+"-"*len(msg))
  92             #
  93             adaopy = adaoBuilder.New()
  94             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"SetSeed":1000, })
  95             adaopy.setBackground         (VectorSerie = 100*[[0,1,2]])
  96             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
  97             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
  98             adaopy.setObservationError   (DiagonalSparseMatrix = "1 2 3")
  99             adaopy.setObservationOperator(Matrix = "1 0 0;0 1 0;0 0 1")
 100             adaopy.setObserver("Analysis",Template="ValuePrinter")
 101             adaopy.execute()
 102             Xa[algo] = adaopy.get("Analysis")[-1]
 103             del adaopy
 104         #
 105         print("")
 106         msg = "Tests des ecarts attendus :"
 107         print(msg+"\n"+"="*len(msg))
 108         verify_similarity_of_algo_results(("3DVAR", "Blue", "ExtendedBlue", "4DVAR", "DerivativeFreeOptimization"), Xa, 5.e-5)
 109         verify_similarity_of_algo_results(("LinearLeastSquares", "NonLinearLeastSquares"), Xa, 5.e-7)
 110         verify_similarity_of_algo_results(("KalmanFilter", "ExtendedKalmanFilter", "UnscentedKalmanFilter"), Xa, 1.e-14)
 111         verify_similarity_of_algo_results(("KalmanFilter", "EnsembleKalmanFilter"), Xa, 5.e-2)
 112         print("  Les resultats obtenus sont corrects.")
 113         print("")
 114         #
 115         return 0
 116
 117     def test2(self):
 118         """Verification de la disponibilite de l'ensemble des algorithmes\n(Utilisation d'un operateur fonctionnel)"""
 119         print(self.test2.__doc__)
 120         Xa = {}
 121         M = numpy.matrix("1 0 0;0 2 0;0 0 3")
 122         def H(x): return M * numpy.asmatrix(numpy.ravel( x )).T
 123         for algo in ("3DVAR", "Blue", "ExtendedBlue", "NonLinearLeastSquares", "DerivativeFreeOptimization"):
 124             print("")
 125             msg = "Algorithme en test : %s"%algo
 126             print(msg+"\n"+"-"*len(msg))
 127             #
 128             adaopy = adaoBuilder.New()
 129             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "Bounds":[[-1,10.],[-1,10.],[-1,10.]]})
 130             adaopy.setBackground         (Vector = [0,1,2])
 131             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 132             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 133             adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
 134             adaopy.setObservationOperator(OneFunction = H)
 135             adaopy.setObserver("Analysis",Template="ValuePrinter")
 136             adaopy.execute()
 137             Xa[algo] = adaopy.get("Analysis")[-1]
 138             del adaopy
 139         #
 140         M = numpy.matrix("1 0 0;0 2 0;0 0 3")
 141         def H(x): return M * numpy.asmatrix(numpy.ravel( x )).T
 142         for algo in ("ExtendedKalmanFilter", "KalmanFilter", "EnsembleKalmanFilter", "UnscentedKalmanFilter", "4DVAR"):
 143             print("")
 144             msg = "Algorithme en test : %s"%algo
 145             print(msg+"\n"+"-"*len(msg))
 146             #
 147             adaopy = adaoBuilder.New()
 148             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "SetSeed":1000})
 149             adaopy.setBackground         (Vector = [0,1,2])
 150             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 151             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 152             adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
 153             adaopy.setObservationOperator(OneFunction = H)
 154             adaopy.setEvolutionError     (ScalarSparseMatrix = 1.)
 155             adaopy.setEvolutionModel     (Matrix = "1 0 0;0 1 0;0 0 1")
 156             adaopy.setObserver("Analysis",Template="ValuePrinter")
 157             adaopy.execute()
 158             Xa[algo] = adaopy.get("Analysis")[-1]
 159             del adaopy
 160         #
 161         M = numpy.matrix("1 0 0;0 1 0;0 0 1")
 162         def H(x): return M * numpy.asmatrix(numpy.ravel( x )).T
 163         for algo in ("ParticleSwarmOptimization", "QuantileRegression", ):
 164             print("")
 165             msg = "Algorithme en test : %s"%algo
 166             print(msg+"\n"+"-"*len(msg))
 167             #
 168             adaopy = adaoBuilder.New()
 169             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"BoxBounds":3*[[-1,3]], "SetSeed":1000})
 170             adaopy.setBackground         (Vector = [0,1,2])
 171             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 172             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 173             adaopy.setObservationError   (DiagonalSparseMatrix = "1 2 3")
 174             adaopy.setObservationOperator(OneFunction = H)
 175             adaopy.setObserver("Analysis",Template="ValuePrinter")
 176             adaopy.execute()
 177             Xa[algo] = adaopy.get("Analysis")[-1]
 178             del adaopy
 179         #
 180         print("")
 181         msg = "Tests des ecarts attendus :"
 182         print(msg+"\n"+"="*len(msg))
 183         verify_similarity_of_algo_results(("3DVAR", "Blue", "ExtendedBlue", "4DVAR", "DerivativeFreeOptimization"), Xa, 5.e-5)
 184         verify_similarity_of_algo_results(("KalmanFilter", "ExtendedKalmanFilter", "UnscentedKalmanFilter"), Xa, 1.e14)
 185         verify_similarity_of_algo_results(("KalmanFilter", "EnsembleKalmanFilter"), Xa, 5.e-2)
 186         print("  Les resultats obtenus sont corrects.")
 187         print("")
 188         #
 189         return 0
 190
 191 def almost_equal_vectors(v1, v2, precision = 1.e-15, msg = ""):
 192     """Comparaison de deux vecteurs"""
 193     print("    Difference maximale %s: %.2e"%(msg, max(abs(v2 - v1))))
 194     return max(abs(v2 - v1)) < precision
 195
 196 def verify_similarity_of_algo_results(serie = [], Xa = {}, precision = 1.e-15):
 197     print("  Comparaisons :")
 198     for algo1 in serie:
 199         for algo2 in serie:
 200             if algo1 is algo2: break
 201             assert almost_equal_vectors( Xa[algo1], Xa[algo2], precision, "entre %s et %s "%(algo1, algo2) )
 202     print("  Algorithmes dont les resultats sont similaires a %.0e : %s\n"%(precision, serie,))
 203     sys.stdout.flush()
 204
 205 #===============================================================================
 206 if __name__ == "__main__":
 207     print("\nAUTODIAGNOSTIC\n==============")
 208     sys.stderr = sys.stdout
 209     unittest.main(verbosity=2)