From e904b7adf04062bc1896ce1030e649f4f06d1e73 Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Mon, 4 Feb 2019 20:40:38 +0100 Subject: [PATCH] Minor updates for module behavior and tests --- src/daComposant/daCore/BasicObjects.py | 2 + test/test1002/Performances.py | 2 + test/test6904/CTestTestfile.cmake | 1 + .../Definition_complete_de_cas_3DVAR.py | 33 ++-- .../Definition_complete_de_cas_Blue.py | 35 ++-- .../Definition_complete_de_cas_NLLS.py | 164 ++++++++++++++++++ 6 files changed, 202 insertions(+), 35 deletions(-) create mode 100644 test/test6904/Definition_complete_de_cas_NLLS.py diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index fe10247..6b37224 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -581,6 +581,7 @@ class Algorithm(object): - IndexOfOptimum : index de l'état optimal courant lors d'itérations - Innovation : l'innovation : d = Y - H(X) - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn) + - JacobianMatrixAtBackground : matrice jacobienne à l'état d'ébauche - JacobianMatrixAtCurrentState : matrice jacobienne à l'état courant - JacobianMatrixAtOptimum : matrice jacobienne à l'optimum - KalmanGainAtOptimum : gain de Kalman à l'optimum @@ -628,6 +629,7 @@ class Algorithm(object): self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum") self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation") self.StoredVariables["InnovationAtCurrentState"] = Persistence.OneVector(name = "InnovationAtCurrentState") + self.StoredVariables["JacobianMatrixAtBackground"] = Persistence.OneMatrix(name = "JacobianMatrixAtBackground") self.StoredVariables["JacobianMatrixAtCurrentState"] = Persistence.OneMatrix(name = "JacobianMatrixAtCurrentState") self.StoredVariables["JacobianMatrixAtOptimum"] = Persistence.OneMatrix(name = "JacobianMatrixAtOptimum") self.StoredVariables["KalmanGainAtOptimum"] = Persistence.OneMatrix(name = "KalmanGainAtOptimum") diff --git a/test/test1002/Performances.py b/test/test1002/Performances.py index db4845b..4b3cde6 100644 --- a/test/test1002/Performances.py +++ b/test/test1002/Performances.py @@ -85,3 +85,5 @@ if __name__ == "__main__": testNumpy01(dimension = 1.e4) testNumpy02(dimension = 3.e3) print("") + print(" Les résultats obtenus sont corrects.") + print("") diff --git a/test/test6904/CTestTestfile.cmake b/test/test6904/CTestTestfile.cmake index 552e54d..87c77ad 100644 --- a/test/test6904/CTestTestfile.cmake +++ b/test/test6904/CTestTestfile.cmake @@ -22,6 +22,7 @@ SET(TEST_NAMES Definition_complete_de_cas_3DVAR Definition_complete_de_cas_Blue + Definition_complete_de_cas_NLLS ) FOREACH(tfile ${TEST_NAMES}) diff --git a/test/test6904/Definition_complete_de_cas_3DVAR.py b/test/test6904/Definition_complete_de_cas_3DVAR.py index 5b8e0a7..35d3135 100644 --- a/test/test6904/Definition_complete_de_cas_3DVAR.py +++ b/test/test6904/Definition_complete_de_cas_3DVAR.py @@ -19,7 +19,7 @@ # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com # # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D -"Verification d'un exemple de la documentation" +"Full definition of a use case for the standard user" import sys import unittest @@ -27,8 +27,8 @@ import numpy # ============================================================================== # -# Construction artificielle d'un exemple de donnees utilisateur -# ------------------------------------------------------------- +# Artificial user data example +# ---------------------------- alpha = 5. beta = 7 gamma = 9.0 @@ -38,11 +38,11 @@ betamin, betamax = 3, 13 gammamin, gammamax = 1.5, 15.5 # def simulation(x): - "Fonction de simulation H pour effectuer Y=H(X)" + "Observation operator H for Y=H(X)" import numpy - __x = numpy.matrix(numpy.ravel(numpy.matrix(x))).T - __H = numpy.matrix("1 0 0;0 2 0;0 0 3; 1 2 3") - return __H * __x + __x = numpy.ravel(x) + __H = numpy.array([[1, 0, 0], [0, 2, 0], [0, 0, 3], [1, 2, 3]]) + return numpy.dot(__H,__x) # def multisimulation( xserie ): yserie = [] @@ -50,17 +50,16 @@ def multisimulation( xserie ): yserie.append( simulation( x ) ) return yserie # -# Observations obtenues par simulation -# ------------------------------------ +# Twin experiment observations +# ---------------------------- observations = simulation((2, 3, 4)) # ============================================================================== class InTest(unittest.TestCase): def test1(self): - print("""Exemple de la doc : - - Exploitation independante des resultats d'un cas de calcul - ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + print(""" + Full definition of a use case for the standard user + +++++++++++++++++++++++++++++++++++++++++++++++++++ """) # import numpy @@ -93,11 +92,11 @@ class InTest(unittest.TestCase): } ) case.set( 'Background', - Vector = numpy.array(Xb), # array, list, tuple, matrix + Vector = Xb, # array, list, tuple, matrix Stored = True, # Bool ) case.set( 'Observation', - Vector = numpy.array(observations), # array, list, tuple, matrix + Vector = observations, # array, list, tuple, matrix Stored = False, # Bool ) case.set( 'BackgroundError', @@ -144,8 +143,8 @@ class InTest(unittest.TestCase): # ---------- ecart = assertAlmostEqualArrays(Xoptimum, [ 2., 3., 4.]) # - print(" L'écart absolu maximal obtenu lors du test est de %.2e."%ecart) - print(" Les résultats obtenus sont corrects.") + print(" The maximal absolute error in the test is of %.2e."%ecart) + print(" The results are correct.") print("") # return Xoptimum diff --git a/test/test6904/Definition_complete_de_cas_Blue.py b/test/test6904/Definition_complete_de_cas_Blue.py index 3daa993..43c96da 100644 --- a/test/test6904/Definition_complete_de_cas_Blue.py +++ b/test/test6904/Definition_complete_de_cas_Blue.py @@ -19,7 +19,7 @@ # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com # # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D -"Verification d'un exemple de la documentation" +"Full definition of a use case for the standard user" import sys import unittest @@ -27,8 +27,8 @@ import numpy # ============================================================================== # -# Construction artificielle d'un exemple de donnees utilisateur -# ------------------------------------------------------------- +# Artificial user data example +# ---------------------------- alpha = 5. beta = 7 gamma = 9.0 @@ -38,11 +38,11 @@ betamin, betamax = 3, 13 gammamin, gammamax = 1.5, 15.5 # def simulation(x): - "Fonction de simulation H pour effectuer Y=H(X)" + "Observation operator H for Y=H(X)" import numpy - __x = numpy.matrix(numpy.ravel(numpy.matrix(x))).T - __H = numpy.matrix("1 0 0;0 2 0;0 0 3; 1 2 3") - return __H * __x + __x = numpy.ravel(x) + __H = numpy.array([[1, 0, 0], [0, 2, 0], [0, 0, 3], [1, 2, 3]]) + return numpy.dot(__H,__x) # def multisimulation( xserie ): yserie = [] @@ -50,17 +50,16 @@ def multisimulation( xserie ): yserie.append( simulation( x ) ) return yserie # -# Observations obtenues par simulation -# ------------------------------------ +# Twin experiment observations +# ---------------------------- observations = simulation((2, 3, 4)) # ============================================================================== class InTest(unittest.TestCase): def test1(self): - print("""Exemple de la doc : - - Exploitation independante des resultats d'un cas de calcul - ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + print(""" + Full definition of a use case for the standard user + +++++++++++++++++++++++++++++++++++++++++++++++++++ """) # import numpy @@ -78,7 +77,7 @@ class InTest(unittest.TestCase): # -------- case = adaoBuilder.New() case.set( 'AlgorithmParameters', - Algorithm = 'Blue', # Mots-clé réservé + Algorithm = 'Blue', # Mots-clé réservé Parameters = { # Dictionnaire "StoreSupplementaryCalculations":[# Liste de mots-clés réservés "CostFunctionJAtCurrentOptimum", @@ -90,11 +89,11 @@ class InTest(unittest.TestCase): } ) case.set( 'Background', - Vector = numpy.array(Xb), # array, list, tuple, matrix + Vector = Xb, # array, list, tuple, matrix Stored = True, # Bool ) case.set( 'Observation', - Vector = numpy.array(observations), # array, list, tuple, matrix + Vector = observations, # array, list, tuple, matrix Stored = False, # Bool ) case.set( 'BackgroundError', @@ -141,8 +140,8 @@ class InTest(unittest.TestCase): # ---------- ecart = assertAlmostEqualArrays(Xoptimum, [ 2., 3., 4.]) # - print(" L'écart absolu maximal obtenu lors du test est de %.2e."%ecart) - print(" Les résultats obtenus sont corrects.") + print("The maximal absolute error in the test is of %.2e."%ecart) + print("The results are correct.") print("") # return Xoptimum diff --git a/test/test6904/Definition_complete_de_cas_NLLS.py b/test/test6904/Definition_complete_de_cas_NLLS.py new file mode 100644 index 0000000..033e355 --- /dev/null +++ b/test/test6904/Definition_complete_de_cas_NLLS.py @@ -0,0 +1,164 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2008-2019 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D +"Full definition of a use case for the standard user" + +import sys +import unittest +import numpy + +# ============================================================================== +# +# Artificial user data example +# ---------------------------- +alpha = 5. +beta = 7 +gamma = 9.0 +# +alphamin, alphamax = 0., 10. +betamin, betamax = 3, 13 +gammamin, gammamax = 1.5, 15.5 +# +def simulation(x): + "Observation operator H for Y=H(X)" + import numpy + __x = numpy.ravel(x) + __H = numpy.array([[1, 0, 0], [0, 2, 0], [0, 0, 3], [1, 2, 3]]) + return numpy.dot(__H,__x) +# +def multisimulation( xserie ): + yserie = [] + for x in xserie: + yserie.append( simulation( x ) ) + return yserie +# +# Twin experiment observations +# ---------------------------- +observations = simulation((2, 3, 4)) + +# ============================================================================== +class InTest(unittest.TestCase): + def test1(self): + print(""" + Full definition of a use case for the standard user + +++++++++++++++++++++++++++++++++++++++++++++++++++ + """) + # + import numpy + from adao import adaoBuilder + # + # Mise en forme des entrees + # ------------------------- + Xb = (alpha, beta, gamma) + Bounds = ( + (alphamin, alphamax), + (betamin, betamax ), + (gammamin, gammamax)) + # + # TUI ADAO + # -------- + case = adaoBuilder.New() + case.set( 'AlgorithmParameters', + Algorithm = 'NonLinearLeastSquares', # Mots-clé réservé + Parameters = { # Dictionnaire + "Bounds":Bounds, # Liste de paires de Real ou de None + "MaximumNumberOfSteps":100, # Int >= 0 + "CostDecrementTolerance":1.e-7, # Real > 0 + "StoreSupplementaryCalculations":[# Liste de mots-clés réservés + "CostFunctionJAtCurrentOptimum", + "CostFunctionJoAtCurrentOptimum", + "CurrentOptimum", + "SimulatedObservationAtCurrentOptimum", + "SimulatedObservationAtOptimum", + ], + } + ) + case.set( 'Background', + Vector = Xb, # array, list, tuple, matrix + Stored = True, # Bool + ) + case.set( 'Observation', + Vector = observations, # array, list, tuple, matrix + Stored = False, # Bool + ) + case.set( 'ObservationError', + Matrix = None, # None ou matrice carrée + ScalarSparseMatrix = 1.0, # None ou Real > 0 + DiagonalSparseMatrix = None, # None ou vecteur + ) + case.set( 'ObservationOperator', + OneFunction = multisimulation, # MultiFonction [Y] = F([X]) + Parameters = { # Dictionnaire + "DifferentialIncrement":0.0001, # Real > 0 + "CenteredFiniteDifference":False, # Bool + }, + InputFunctionAsMulti = True, # Bool + ) + case.set( 'Observer', + Variable = "CurrentState", # Mot-clé + Template = "ValuePrinter", # Mot-clé + String = None, # None ou code Python + Info = None, # None ou string + + ) + case.execute() + # + # Exploitation independante + # ------------------------- + Xbackground = case.get("Background") + Xoptimum = case.get("Analysis")[-1] + FX_at_optimum = case.get("SimulatedObservationAtOptimum")[-1] + J_values = case.get("CostFunctionJAtCurrentOptimum")[:] + print("") + print("Number of internal iterations...: %i"%len(J_values)) + print("Initial state...................: %s"%(numpy.ravel(Xbackground),)) + print("Optimal state...................: %s"%(numpy.ravel(Xoptimum),)) + print("Simulation at optimal state.....: %s"%(numpy.ravel(FX_at_optimum),)) + print("") + # + # Fin du cas + # ---------- + ecart = assertAlmostEqualArrays(Xoptimum, [ 2., 3., 4.]) + # + print("The maximal absolute error in the test is of %.2e."%ecart) + print("The results are correct.") + print("") + # + return Xoptimum + +# ============================================================================== +def assertAlmostEqualArrays(first, second, places=7, msg=None, delta=None): + "Compare two vectors, like unittest.assertAlmostEqual" + import numpy + if msg is not None: + print(msg) + if delta is not None: + if ( (numpy.asarray(first) - numpy.asarray(second)) > float(delta) ).any(): + raise AssertionError("%s != %s within %s places"%(first,second,delta)) + else: + if ( (numpy.asarray(first) - numpy.asarray(second)) > 10**(-int(places)) ).any(): + raise AssertionError("%s != %s within %i places"%(first,second,places)) + return max(abs(numpy.asarray(first) - numpy.asarray(second))) + +# ============================================================================== +if __name__ == '__main__': + print("\nAUTODIAGNOSTIC\n==============") + unittest.main() -- 2.39.2