test/test6901/Verification_des_Assimilation_Algorithms.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2023 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, 2.e-1)
 112         print("  Les resultats obtenus sont corrects.")
 113         print("")
 114
 115     def test2(self):
 116         """Verification de la disponibilite de l'ensemble des algorithmes\n(Utilisation d'un operateur fonctionnel)"""
 117         print(self.test2.__doc__)
 118         Xa = {}
 119         M = numpy.diag([1.,2.,3.])
 120         def H(x): return M @ numpy.ravel( x )
 121         for algo in ("3DVAR", "Blue", "ExtendedBlue", "NonLinearLeastSquares", "DerivativeFreeOptimization"):
 122             print("")
 123             msg = "Algorithme en test : %s"%algo
 124             print(msg+"\n"+"-"*len(msg))
 125             #
 126             adaopy = adaoBuilder.New()
 127             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "Bounds":[[-1,10.],[-1,10.],[-1,10.]]})
 128             adaopy.setBackground         (Vector = [0,1,2])
 129             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 130             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 131             adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
 132             adaopy.setObservationOperator(OneFunction = H)
 133             adaopy.setObserver("Analysis",Template="ValuePrinter")
 134             adaopy.execute()
 135             Xa[algo] = adaopy.get("Analysis")[-1]
 136             del adaopy
 137         #
 138         M = numpy.diag([1.,2.,3.])
 139         def H(x): return M @ numpy.ravel( x )
 140         for algo in ("ExtendedKalmanFilter", "KalmanFilter", "EnsembleKalmanFilter", "UnscentedKalmanFilter", "4DVAR"):
 141             print("")
 142             msg = "Algorithme en test : %s"%algo
 143             print(msg+"\n"+"-"*len(msg))
 144             #
 145             adaopy = adaoBuilder.New()
 146             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "SetSeed":1000})
 147             adaopy.setBackground         (Vector = [0,1,2])
 148             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 149             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 150             adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
 151             adaopy.setObservationOperator(OneFunction = H)
 152             adaopy.setEvolutionError     (ScalarSparseMatrix = 1.)
 153             adaopy.setEvolutionModel     (Matrix = "1 0 0;0 1 0;0 0 1")
 154             adaopy.setObserver("Analysis",Template="ValuePrinter")
 155             adaopy.execute()
 156             Xa[algo] = adaopy.get("Analysis")[-1]
 157             del adaopy
 158         #
 159         M = numpy.identity(3)
 160         def H(x): return M @ numpy.ravel( x )
 161         for algo in ("ParticleSwarmOptimization", "QuantileRegression", ):
 162             print("")
 163             msg = "Algorithme en test : %s"%algo
 164             print(msg+"\n"+"-"*len(msg))
 165             #
 166             adaopy = adaoBuilder.New()
 167             adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"BoxBounds":3*[[-1,3]], "SetSeed":1000})
 168             adaopy.setBackground         (Vector = [0,1,2])
 169             adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 170             adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 171             adaopy.setObservationError   (DiagonalSparseMatrix = "1 2 3")
 172             adaopy.setObservationOperator(OneFunction = H)
 173             adaopy.setObserver("Analysis",Template="ValuePrinter")
 174             adaopy.execute()
 175             Xa[algo] = adaopy.get("Analysis")[-1]
 176             del adaopy
 177         #
 178         print("")
 179         msg = "Tests des ecarts attendus :"
 180         print(msg+"\n"+"="*len(msg))
 181         verify_similarity_of_algo_results(("3DVAR", "Blue", "ExtendedBlue", "4DVAR", "DerivativeFreeOptimization"), Xa, 5.e-5)
 182         verify_similarity_of_algo_results(("KalmanFilter", "ExtendedKalmanFilter", "UnscentedKalmanFilter"), Xa, 2.e-14)
 183         verify_similarity_of_algo_results(("KalmanFilter", "EnsembleKalmanFilter"), Xa, 2e-1)
 184         print("  Les resultats obtenus sont corrects.")
 185         print("")
 186
 187 def almost_equal_vectors(v1, v2, precision = 1.e-15, msg = ""):
 188     """Comparaison de deux vecteurs"""
 189     print("    Difference maximale %s: %.2e"%(msg, max(abs(v2 - v1))))
 190     return max(abs(v2 - v1)) < precision
 191
 192 def verify_similarity_of_algo_results(serie = [], Xa = {}, precision = 1.e-15):
 193     print("  Comparaisons :")
 194     for algo1 in serie:
 195         for algo2 in serie:
 196             if algo1 is algo2: break
 197             assert almost_equal_vectors( Xa[algo1], Xa[algo2], precision, "entre %s et %s "%(algo1, algo2) )
 198     print("  Algorithmes dont les resultats sont similaires a %.0e : %s\n"%(precision, serie,))
 199     sys.stdout.flush()
 200
 201 #===============================================================================
 202 if __name__ == "__main__":
 203     print("\nAUTODIAGNOSTIC\n==============")
 204     sys.stderr = sys.stdout
 205     unittest.main(verbosity=2)