test/test6901/Verification_des_Assimilation_Algorithms.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2018 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 # ==============================================================================
  25 import adaoBuilder, numpy, sys
  26 def test1():
  27     """Verification de la disponibilite de l'ensemble des algorithmes\n(Utilisation d'un operateur matriciel)"""
  28     print(test1.__doc__)
  29     Xa = {}
  30     for algo in ("3DVAR", "Blue", "ExtendedBlue", "LinearLeastSquares", "NonLinearLeastSquares", "DerivativeFreeOptimization"):
  31         print("")
  32         msg = "Algorithme en test : %s"%algo
  33         print(msg+"\n"+"-"*len(msg))
  34         #
  35         adaopy = adaoBuilder.New()
  36         adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "Bounds":[[-1,10.],[-1,10.],[-1,10.]]})
  37         adaopy.setBackground         (Vector = [0,1,2])
  38         adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
  39         adaopy.setObservation        (Vector = [0.5,1.5,2.5])
  40         adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
  41         adaopy.setObservationOperator(Matrix = "1 0 0;0 2 0;0 0 3")
  42         adaopy.setObserver("Analysis",Template="ValuePrinter")
  43         adaopy.execute()
  44         Xa[algo] = adaopy.get("Analysis")[-1]
  45         del adaopy
  46     #
  47     for algo in ("ExtendedKalmanFilter", "KalmanFilter", "UnscentedKalmanFilter", "EnsembleKalmanFilter", "4DVAR"):
  48         print("")
  49         msg = "Algorithme en test : %s"%algo
  50         print(msg+"\n"+"-"*len(msg))
  51         #
  52         adaopy = adaoBuilder.New()
  53         adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "SetSeed":1000})
  54         adaopy.setBackground         (Vector = [0,1,2])
  55         adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
  56         adaopy.setObservation        (Vector = [0.5,1.5,2.5])
  57         adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
  58         adaopy.setObservationOperator(Matrix = "1 0 0;0 2 0;0 0 3")
  59         adaopy.setEvolutionError     (ScalarSparseMatrix = 1.)
  60         adaopy.setEvolutionModel     (Matrix = "1 0 0;0 1 0;0 0 1")
  61         adaopy.setObserver("Analysis",Template="ValuePrinter")
  62         adaopy.execute()
  63         Xa[algo] = adaopy.get("Analysis")[-1]
  64         del adaopy
  65     #
  66     for algo in ("ParticleSwarmOptimization", "QuantileRegression", ):
  67         print("")
  68         msg = "Algorithme en test : %s"%algo
  69         print(msg+"\n"+"-"*len(msg))
  70         #
  71         adaopy = adaoBuilder.New()
  72         adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"BoxBounds":3*[[-1,3]], "SetSeed":1000})
  73         adaopy.setBackground         (Vector = [0,1,2])
  74         adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
  75         adaopy.setObservation        (Vector = [0.5,1.5,2.5])
  76         adaopy.setObservationError   (DiagonalSparseMatrix = "1 2 3")
  77         adaopy.setObservationOperator(Matrix = "1 0 0;0 1 0;0 0 1")
  78         adaopy.setObserver("Analysis",Template="ValuePrinter")
  79         adaopy.execute()
  80         Xa[algo] = adaopy.get("Analysis")[-1]
  81         del adaopy
  82     #
  83     for algo in ("EnsembleBlue", ):
  84         print("")
  85         msg = "Algorithme en test : %s"%algo
  86         print(msg+"\n"+"-"*len(msg))
  87         #
  88         adaopy = adaoBuilder.New()
  89         adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"SetSeed":1000, })
  90         adaopy.setBackground         (VectorSerie = 100*[[0,1,2]])
  91         adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
  92         adaopy.setObservation        (Vector = [0.5,1.5,2.5])
  93         adaopy.setObservationError   (DiagonalSparseMatrix = "1 2 3")
  94         adaopy.setObservationOperator(Matrix = "1 0 0;0 1 0;0 0 1")
  95         adaopy.setObserver("Analysis",Template="ValuePrinter")
  96         adaopy.execute()
  97         Xa[algo] = adaopy.get("Analysis")[-1]
  98         del adaopy
  99     #
 100     print("")
 101     msg = "Tests des ecarts attendus :"
 102     print(msg+"\n"+"="*len(msg))
 103     verify_similarity_of_algo_results(("3DVAR", "Blue", "ExtendedBlue", "4DVAR", "DerivativeFreeOptimization"), Xa, 5.e-5)
 104     verify_similarity_of_algo_results(("LinearLeastSquares", "NonLinearLeastSquares"), Xa, 5.e-7)
 105     verify_similarity_of_algo_results(("KalmanFilter", "ExtendedKalmanFilter", "UnscentedKalmanFilter"), Xa, 1.e-14)
 106     verify_similarity_of_algo_results(("KalmanFilter", "EnsembleKalmanFilter"), Xa, 5.e-2)
 107     print("  Les resultats obtenus sont corrects.")
 108     print("")
 109     #
 110     return 0
 111
 112 def test2():
 113     """Verification de la disponibilite de l'ensemble des algorithmes\n(Utilisation d'un operateur fonctionnel)"""
 114     print(test2.__doc__)
 115     Xa = {}
 116     M = numpy.matrix("1 0 0;0 2 0;0 0 3")
 117     def H(x): return M * numpy.asmatrix(numpy.ravel( x )).T
 118     for algo in ("3DVAR", "Blue", "ExtendedBlue", "NonLinearLeastSquares", "DerivativeFreeOptimization"):
 119         print("")
 120         msg = "Algorithme en test : %s"%algo
 121         print(msg+"\n"+"-"*len(msg))
 122         #
 123         adaopy = adaoBuilder.New()
 124         adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "Bounds":[[-1,10.],[-1,10.],[-1,10.]]})
 125         adaopy.setBackground         (Vector = [0,1,2])
 126         adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 127         adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 128         adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
 129         adaopy.setObservationOperator(OneFunction = H)
 130         adaopy.setObserver("Analysis",Template="ValuePrinter")
 131         adaopy.execute()
 132         Xa[algo] = adaopy.get("Analysis")[-1]
 133         del adaopy
 134     #
 135     M = numpy.matrix("1 0 0;0 2 0;0 0 3")
 136     def H(x): return M * numpy.asmatrix(numpy.ravel( x )).T
 137     for algo in ("ExtendedKalmanFilter", "KalmanFilter", "EnsembleKalmanFilter", "UnscentedKalmanFilter", "4DVAR"):
 138         print("")
 139         msg = "Algorithme en test : %s"%algo
 140         print(msg+"\n"+"-"*len(msg))
 141         #
 142         adaopy = adaoBuilder.New()
 143         adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"EpsilonMinimumExponent":-10, "SetSeed":1000})
 144         adaopy.setBackground         (Vector = [0,1,2])
 145         adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 146         adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 147         adaopy.setObservationError   (DiagonalSparseMatrix = "1 1 1")
 148         adaopy.setObservationOperator(OneFunction = H)
 149         adaopy.setEvolutionError     (ScalarSparseMatrix = 1.)
 150         adaopy.setEvolutionModel     (Matrix = "1 0 0;0 1 0;0 0 1")
 151         adaopy.setObserver("Analysis",Template="ValuePrinter")
 152         adaopy.execute()
 153         Xa[algo] = adaopy.get("Analysis")[-1]
 154         del adaopy
 155     #
 156     M = numpy.matrix("1 0 0;0 1 0;0 0 1")
 157     def H(x): return M * numpy.asmatrix(numpy.ravel( x )).T
 158     for algo in ("ParticleSwarmOptimization", "QuantileRegression", ):
 159         print("")
 160         msg = "Algorithme en test : %s"%algo
 161         print(msg+"\n"+"-"*len(msg))
 162         #
 163         adaopy = adaoBuilder.New()
 164         adaopy.setAlgorithmParameters(Algorithm=algo, Parameters={"BoxBounds":3*[[-1,3]], "SetSeed":1000})
 165         adaopy.setBackground         (Vector = [0,1,2])
 166         adaopy.setBackgroundError    (ScalarSparseMatrix = 1.)
 167         adaopy.setObservation        (Vector = [0.5,1.5,2.5])
 168         adaopy.setObservationError   (DiagonalSparseMatrix = "1 2 3")
 169         adaopy.setObservationOperator(OneFunction = H)
 170         adaopy.setObserver("Analysis",Template="ValuePrinter")
 171         adaopy.execute()
 172         Xa[algo] = adaopy.get("Analysis")[-1]
 173         del adaopy
 174     #
 175     print("")
 176     msg = "Tests des ecarts attendus :"
 177     print(msg+"\n"+"="*len(msg))
 178     verify_similarity_of_algo_results(("3DVAR", "Blue", "ExtendedBlue", "4DVAR", "DerivativeFreeOptimization"), Xa, 5.e-5)
 179     verify_similarity_of_algo_results(("KalmanFilter", "ExtendedKalmanFilter", "UnscentedKalmanFilter"), Xa, 1.e14)
 180     verify_similarity_of_algo_results(("KalmanFilter", "EnsembleKalmanFilter"), Xa, 5.e-2)
 181     print("  Les resultats obtenus sont corrects.")
 182     print("")
 183     #
 184     return 0
 185
 186 def almost_equal_vectors(v1, v2, precision = 1.e-15, msg = ""):
 187     """Comparaison de deux vecteurs"""
 188     print("    Difference maximale %s: %.2e"%(msg, max(abs(v2 - v1))))
 189     return max(abs(v2 - v1)) < precision
 190
 191 def verify_similarity_of_algo_results(serie = [], Xa = {}, precision = 1.e-15):
 192     print("  Comparaisons :")
 193     for algo1 in serie:
 194         for algo2 in serie:
 195             if algo1 is algo2: break
 196             assert almost_equal_vectors( Xa[algo1], Xa[algo2], precision, "entre %s et %s "%(algo1, algo2) )
 197     print("  Algorithmes dont les resultats sont similaires a %.0e : %s\n"%(precision, serie,))
 198     sys.stdout.flush()
 199
 200 #===============================================================================
 201 if __name__ == "__main__":
 202     print('\nAUTODIAGNOSTIC\n')
 203     test1()
 204     test2()