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