salomeresdir = $(prefix)/share/salome/resources/${MODULE_NAME}
# Directory for installing tests files
-testsdir = $(prefix)/tests
+testsdir = $(prefix)/tests
+testsdasalomedir = $(testsdir)/daSalome
+testsplateformedir = $(testsdir)/daComposant/Plateforme
# Directories for installing admin files
admlocaldir = $(prefix)/adm_local
src/tests/Makefile
src/tests/daSalome/Makefile
src/tests/daSalome/test000_Blue_AnalysisFile.py
+ src/tests/daComposant/Makefile
+ src/tests/daComposant/Plateforme/Makefile
+ src/tests/daComposant/Plateforme/context.py
bin/Makefile
bin/qtEficas_adao_study.py
doc/Makefile
+SUBDIRS = daComposant
+
if SALOME_MODULE
-SUBDIRS = daSalome
+SUBDIRS += daSalome
endif
--- /dev/null
+SUBDIRS = Plateforme
+
--- /dev/null
+include $(top_srcdir)/adm_local/make_common_starter.am
+
+DATA_INST = \
+ test000_Etude_ADD.py \
+ test001_Blue.py \
+ test002_Blue.py \
+ test003_Blue.py \
+ test004_Blue.py \
+ test005_Blue.py \
+ test006_Blue_ReduceVariance.py \
+ test007_Blue.py \
+ test008_Kalman.py \
+ test009_Blue_EspaceVeta.py \
+ test010_Kalman_sur_trajectoire_1D.py \
+ test012_LinearLeastSquares.py \
+ test013_EnsembleBlue.py \
+ test014_Persistence_et_Blue_en_boucle.py \
+ test015_3DVAR.py \
+ test016_3DVAR_par_fonction.py \
+ test017_3DVAR_par_fonction.py \
+ test018_3DVAR_par_fonction_avec_bornes.py \
+ test101_RMS.py \
+ test102_PlotVector.py \
+ test103_PlotVectors.py \
+ test104_HomogeneiteKhi2.py \
+ test105_VarianceOrder.py \
+ test106_CompareVarianceFisher.py \
+ test107_GaussianAdequation.py \
+ test108_StopReductionVariance.py \
+ test109_HLinearity.py \
+ test300_Optimize_CG.py \
+ test301_Optimize_BFGS.py \
+ test302_Optimize_LBFGS.py \
+ test303_Optimize_TNC.py \
+ test901_Informations_dynamiques.py \
+ test902_Demonstrations_de_logging.py
+
+GEN_DATA_INST = context.py
+
+testsplateforme_DATA = ${DATA_INST} ${GEN_DATA_INST}
+
+EXTRA_DIST = ${DATA_INST} context.py.in
+
--- /dev/null
+#-*-coding:utf-8*-
+
+# Permet de définir ou se trouve daCore
+import sys ; sys.path.insert(0, "@PYTHON_SITE@/salome/daCore")
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que si les covariances d'erreur B et R sont identiques et
+ unitaires, l'analyse est située au milieu de l'ébauche [0,1,2] et de
+ l'observation [0.5,1.5,2.5].
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+import logging
+# logging.getLogger().setLevel(logging.DEBUG)
+
+#===============================================================================
+def test(precision = 1.e-13):
+ """
+ Cas-test vérifiant que si les covariances d'erreur B et R sont identiques et
+ unitaires, l'analyse est située au milieu de l'ébauche [0,1,2] et de
+ l'observation [0.5,1.5,2.5].
+ """
+ #
+ # Définition de l'étude d'assimilation
+ # ------------------------------------
+ ADD = AssimilationStudy("Ma premiere etude")
+ #
+ ADD.setBackground (asVector = [0,1,2])
+ ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
+ ADD.setObservation (asVector = [0.5,1.5,2.5])
+ ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
+ ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ Xa = ADD.get("Analysis")
+ print
+ print " Nombre d'analyses :",Xa.stepnumber()
+ print " Analyse résultante :",Xa.valueserie(0)
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(numpy.array(Xa.valueserie(0))-numpy.array([0.25, 1.25, 2.25])) > precision:
+ raise ValueError("Résultat du test erroné")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant sur le Blue que si l'erreur est nulle, l'incrément
+ d'analyse est nul.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(precision = 1.e-13, dimension = 3):
+ """
+ Cas-test vérifiant sur le Blue que si l'erreur est nulle, l'incrément
+ d'analyse est nul.
+ """
+ #
+ # Définition des données
+ # ----------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.zeros((dimension,))).T
+ Eb = numpy.matrix(numpy.zeros((dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ xb = xt + Eb
+ yo = H * xt + Eo
+ #
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Analyse
+ # -------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(xa - xb)) > precision:
+ raise ValueError("Résultat du test erroné (1)")
+ elif max(abs(d)) > precision:
+ raise ValueError("Résultat du test erroné (2)")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ # numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que l'application des coefficients de correction so et sb
+ conduit à des matrices R et B pour lesquelles ces coefficients sont unitaires.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(precision = 1.e-13, dimension = 3):
+ """
+ Cas-test vérifiant que l'application des coefficients de correction so et sb
+ conduit à des matrices R et B pour lesquelles ces coefficients sont unitaires.
+ """
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ xb = xt + Eb
+ yo = H * xt + Eo
+ #
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices d'erreurs
+ # ---------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ SigmaObs2 = float( numpy.dot(d,(yo-numpy.dot(H,xa)).A1) / R.trace() )
+ SigmaBck2 = float( numpy.dot(d,numpy.dot(H,(xa - xb)).A1) /(H * B * H.T).trace() )
+ #
+ # Analyse BLUE avec correction des matrices R et B
+ # Attention : ce second calcul de BLUE avec le meme objet ADD
+ # conduit à stocker les résultats dans le second step,
+ # donc il faut appeller "valueserie(1)"
+ # ------------------------------------------------
+ ADD.setBackgroundError (asCovariance = SigmaBck2*B )
+ ADD.setObservationError(asCovariance = SigmaObs2*R )
+ ADD.analyze()
+ new_xa = numpy.array(ADD.get("Analysis").valueserie(1))
+ new_d = numpy.array(ADD.get("Innovation").valueserie(1))
+ new_SigmaObs2 = float( numpy.dot(new_d,(yo-numpy.dot(H,new_xa)).A1) / (SigmaObs2*R.trace()) )
+ new_SigmaBck2 = float( numpy.dot(new_d,numpy.dot(H,(new_xa - xb)).A1) /(H * (SigmaBck2*B) * H.T).trace() )
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(xa - new_xa)) > precision:
+ raise ValueError("Résultat du test erroné (1)")
+ elif max(abs(d - new_d)) > precision:
+ raise ValueError("Résultat du test erroné (2)")
+ elif abs(new_SigmaObs2-1.) > precision:
+ print "new_SigmaObs2 =",new_SigmaObs2
+ raise ValueError("Résultat du test erroné (3)")
+ elif abs(new_SigmaBck2-1.) > precision :
+ print "new_SigmaBck2 =",new_SigmaBck2
+ raise ValueError("Résultat du test erroné (4)")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ # numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que si l'erreur sur le background est nulle et que
+ l'erreur sur les observations est connue, alors l'analyse donne le "milieu"
+ du background et des observations.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(precision = 1.e-13, dimension = 3):
+ """
+ Cas-test vérifiant que si l'erreur sur le background est nulle et que
+ l'erreur sur les observations est connue, alors l'analyse donne le "milieu"
+ du background et des observations.
+ """
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eb = numpy.matrix(numpy.zeros((dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ xb = xt + Eb
+ yo = H * xt + Eo
+ #
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ Xa = ADD.get("Analysis")
+ xa = numpy.matrix(Xa.valueserie(0)).T
+ SigmaObs2 = ADD.get("SigmaObs2")
+ SigmaBck2 = ADD.get("SigmaBck2")
+ d = ADD.get("Innovation")
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(xa.A1 - xb - Eo.A1/2.)) > precision:
+ raise ValueError("Résultat du test erroné (1)")
+ elif max(abs(yo - (H * xa).A1 - Eo.A1/2.)) > precision:
+ raise ValueError("Résultat du test erroné (2)")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ # numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que si l'erreur sur le background est nulle et que
+ l'erreur sur les observations est connue, alors l'analyse donne le "milieu"
+ du background et des observations.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(precision = 1.e-13, dimension = 3):
+ """
+ Cas-test vérifiant que si l'on rajoute l'évaluation de l'opérateur
+ d'observation au background, on obtient la même valeur que pour le BLUE
+ normal.
+ """
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eb = numpy.matrix(numpy.zeros((dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ xb = xt + Eb
+ yo = H * xt + Eo
+ Hxb = H*xb
+ #
+ xb = xb.A1
+ yo = yo.A1
+ HXb = Hxb.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension))
+ B = numpy.matrix(numpy.core.identity(dimension))
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ SigmaObs2 = float( numpy.dot(d,(yo-numpy.dot(H,xa)).A1) / R.trace() )
+ SigmaBck2 = float( numpy.dot(d,numpy.dot(H,(xa - xb)).A1) /(H * B * H.T).trace() )
+ #
+ # Analyse BLUE avec une évaluation au point Xb
+ # Attention : ce second calcul de BLUE avec le meme objet ADD
+ # conduit à stocker les résultats dans le second step,
+ # donc il faut appeller "valueserie(1)"
+ # ------------------------------------------------
+ ADD.setObservationOperator(asMatrix = H, appliedToX = {"HXb":HXb} )
+ ADD.analyze()
+ #
+ new_xa = numpy.array(ADD.get("Analysis").valueserie(1))
+ new_d = numpy.array(ADD.get("Innovation").valueserie(1))
+ new_SigmaObs2 = float( numpy.dot(new_d,(yo-numpy.dot(H,new_xa)).A1) / R.trace() )
+ new_SigmaBck2 = float( numpy.dot(new_d,numpy.dot(H,(new_xa - xb)).A1) /(H * B * H.T).trace() )
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(xa - new_xa)) > precision:
+ raise ValueError("Résultat du test erroné (1)")
+ elif max(abs(d - new_d)) > precision:
+ raise ValueError("Résultat du test erroné (2)")
+ elif abs((new_SigmaObs2-SigmaObs2)/SigmaObs2) > precision:
+ raise ValueError("Résultat du test erroné (3)")
+ elif abs((new_SigmaBck2-SigmaBck2)/SigmaBck2) > precision :
+ raise ValueError("Résultat du test erroné (4)")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ # numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que si H est l'identité et que les matrices de covariance
+ d'erreurs sont liées par R = alpha * B, alors l'ecart type de OMA est
+ proportionnel a l'ecart type de l'innovation d selon la relation :
+ rms(OMA) = alpha/(1. + alpha) rms(d)
+"""
+__author__ = "Sophie RICCI - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(precision = 1.e-13, dimension = 3, alpha = 2.):
+ """
+ Cas-test vérifiant que si H est l'identité et que les matrices de covariance
+ d'erreurs sont liées par R = alpha * B, alors l'ecart type de OMA est
+ proportionnel a l'ecart type de l'innovation d selon la relation :
+ rms(OMA) = alpha/(1. + alpha) rms(d)
+ """
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eb = numpy.matrix(numpy.zeros((dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ xb = xt + Eb
+ yo = H * xt + Eo
+ #
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ R = alpha * B
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ xa = ADD.get("Analysis").valueserie(0)
+ d = ADD.get("Innovation").valueserie(0)
+ #
+ # Calcul RMS pour d et OMA
+ # ------------------------
+ ADD.setDiagnostic("RMS",
+ name = "Calcul de la RMS sur l'innovation et OMA",
+ )
+ RMS = ADD.get("Calcul de la RMS sur l'innovation et OMA")
+ #
+ # La RMS de l'innovation d
+ # ------------------------
+ RMS.calculate(d,numpy.zeros(len(d)))
+ # Le calcul ci-dessus doit être identique à : RMS.calculate(xb,yo)
+ #
+ # La RMS de l'écart OMA
+ # ---------------------
+ RMS.calculate(xa,yo)
+ #
+ # Vérification du résultat
+ # ------------------------
+ if (RMS.valueserie(1) - (alpha/(1. + alpha)) * RMS.valueserie(0)) > precision:
+ raise ValueError("Résultat du test erroné")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Vérification de la réduction de variance opérée par un BLUE lors de
+ l'analyse
+"""
+__author__ = "Sophie RICCI - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(dimension = 10):
+ """
+ Cas-test vérifiant que l'analyse BLUE permet de réduire la variance entre
+ les écarts OMB et les écarts OMA
+ """
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ xb = xt + Eb
+ yo = H * xt + Eo
+ #
+ xb = xb
+ yo = yo
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ OMA = yo.A1 - xa
+ #
+ # Application du test
+ # -------------------
+ ADD.setDiagnostic("ReduceVariance",
+ name = "Reduction de la variance entre OMB et OMA")
+ #
+ D = ADD.get("Reduction de la variance entre OMB et OMA")
+ #
+ D.calculate( vectorOMB = d, vectorOMA = OMA )
+ #
+ # Vérification du résultat
+ # ------------------------
+ if not D.valueserie(0) :
+ raise ValueError("Résultat du test erroné (1)")
+ else :
+ print test.__doc__
+ print " Test correct"
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant les relations d'ordre attendues sur les écarts RMS entre
+ les valeurs analysees et la valeur vraie, pour 3 analyses BLUE réalisées
+ avec des poids extrêmes dans R et B
+"""
+__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(dimension = 3):
+ """
+ Cas-test vérifiant les relations d'ordre attendues sur les écarts RMS entre
+ les valeurs analysees et la valeur vraie, pour 3 analyses BLUE réalisées
+ avec des poids extrêmes dans R et B
+ """
+ print test.__doc__
+ #
+ # Définition des données
+ # ----------------------
+ xt = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ H = numpy.matrix(numpy.core.identity(dimension))
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ xb = xt + Eb
+ yo = H * xt
+ xt = xt.A1
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setObservation (asVector = yo )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservationOperator(asMatrix = H )
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ # Définition des matrices de covariances d'erreur : ébauche parfaite
+ # ------------------------------------------------------------------
+ alpha1 = 10000.0
+ R = alpha1 * B
+ ADD.setObservationError (asCovariance = R )
+ ADD.analyze()
+ x1 = ADD.get("Analysis").valueserie(0)
+ #
+ # Définition des matrices de covariances d'erreurs : poids identiques
+ # -------------------------------------------------------------------
+ alpha2 = 1.0
+ R = alpha2 * B
+ ADD.setObservationError (asCovariance = R )
+ ADD.analyze()
+ x2 = ADD.get("Analysis").valueserie(1)
+ #
+ # Définition des matrices de covariances d'erreurs : observations parfaites
+ # -------------------------------------------------------------------------
+ alpha3 = 0.0001
+ R = alpha3 * B
+ ADD.setObservationError (asCovariance = R )
+ ADD.analyze()
+ x3 = ADD.get("Analysis").valueserie(2)
+ #
+ # Calcul des écarts RMS
+ # ---------------------
+ ADD.setDiagnostic("RMS", "Calcul de la RMS entre analyse et yo")
+ RMS = ADD.get("Calcul de la RMS entre analyse et yo")
+ #
+ RMS.calculate(x1,yo)
+ RMS.calculate(x2,yo)
+ RMS.calculate(x3,yo)
+ RMS_yo_x1 = RMS.valueserie(0)
+ RMS_yo_x2 = RMS.valueserie(1)
+ RMS_yo_x3 = RMS.valueserie(2)
+ #
+ print " Cas ébauche parfaite : R/B = %.1e"%alpha1,"RMS = %.7f"%RMS_yo_x1
+ print " Cas poids identiques : R/B = %.1e"%alpha2,"RMS = %.7f"%RMS_yo_x2
+ print " Cas observations parfaites : R/B = %.1e"%alpha3,"RMS = %.7f"%RMS_yo_x3
+ if ( (RMS_yo_x3 <= RMS_yo_x2) and (RMS_yo_x2 <= RMS_yo_x1) ) :
+ print " La reponse de l'assimilation est cohérente avec la modification du rapport B/R."
+ print
+ print " Test correct"
+ print
+ else :
+ raise ValueError("Résultat du test erroné")
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant le fonctionnement du filtre de Kalman sur un système
+ dynamique de trajectoire 1D constante
+"""
+__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+import Persistence
+
+#===============================================================================
+def test(dimension = 3):
+ """
+ Cas-test vérifiant le fonctionnement du filtre de Kalman sur un système
+ dynamique de trajectoire 1D constante
+ """
+ print test.__doc__
+ #
+ # Définition des données
+ # ----------------------
+ a_size = (dimension,)
+ #
+ # Valeur vraie
+ xt = -0.4
+ Xt = Persistence.OneScalar("Valeur vraie", basetype=float)
+ Xt.store(xt)
+ for i in range(dimension):
+ Xt.store(xt)
+ #
+ # Observations bruitées
+ yo = numpy.random.normal(xt, 0.1, size=a_size)
+ Yo = Persistence.OneScalar("Observations", basetype=float)
+ Yo.store(0.)
+ for v in yo:
+ Yo.store(v)
+ #
+ # Création de l'étude et résolution
+ # ---------------------------------
+ ADD = AssimilationStudy("Assimilation temporelle de Kalman")
+ #
+ ADD.setBackground (asVector = "0.")
+ ADD.setBackgroundError (asCovariance = "1.")
+ #
+ ADD.setObservationOperator(asMatrix = "1.")
+ ADD.setObservation (asPersistentVector = Yo)
+ ADD.setObservationError (asCovariance = "1.e-2")
+ #
+ ADD.setEvolutionModel (asMatrix = "1")
+ ADD.setEvolutionError (asCovariance = "1.e-5")
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Kalman")
+ #
+ ADD.analyze()
+ #
+ Xa = ADD.get("Analysis")
+ print " Valeur vraie visée........................:",xt
+ print " Ebauche, i.e. valeur initiale d'analyse...:",Xa.valueserie(0)[0]
+ print " Nombre d'analyses (sans l'ébauche)........:",Xa.stepnumber()-1
+ print " Moyenne des analyses......................:",Xa.stepmean()
+ #
+ # Biais des erreurs
+ EpsY = []
+ for i in range(Yo.stepnumber()):
+ EpsY.append(Yo.valueserie(i) - Xt.valueserie(i))
+ print " Biais des erreurs <Obs-Vraie>.............:",numpy.array(EpsY).mean()
+ print " Variance des erreurs <Obs-Vraie>..........:",numpy.array(EpsY).var()
+ EpsY = []
+ for i in range(Xa.stepnumber()):
+ EpsY.append(Xa.valueserie(i)[0] - Xt.valueserie(i))
+ print " Biais des erreurs <Ana-Vraie>.............:",numpy.array(EpsY).mean()
+ print " Variance des erreurs <Ana-Vraie>..........:",numpy.array(EpsY).var()
+ print
+ #
+ ADD.setDiagnostic("PlotVectors", "Affichage de Xa et Xt")
+ MonPlot = ADD.get("Affichage de Xa et Xt")
+ MonPlot.calculate(
+ ( [ x[0] for x in Xa.valueserie()], Xt.valueserie(), Yo.valueserie() ),
+ title = "Analyse de Kalman sur trajectoire constante",
+ ltitle = ["Analyse", "Valeur vraie", "Observations"],
+ filename = "kalman_sur_trajectoire_constante.ps",
+ pause = False,
+ )
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Vérification du calcul de BLUE dans l'espace des états plutôt que dans
+ l'espace des observations.
+"""
+__author__ = "Sophie RICCI - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(precision = 1.e-13):
+ """
+ Vérification du calcul de BLUE dans l'espace des états plutôt que dans
+ l'espace des observations.
+ """
+ #
+ # Définition des données
+ # ------------------------------------------
+ H = numpy.matrix(([1., 1.])).T
+ #
+ xb = 6.
+ xt = 3.
+ yo = H * xt
+ #
+ dimension = yo.size
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ B = numpy.matrix(([1.]))
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(xa - 4.)) > precision:
+ raise ValueError("Résultat du test erroné")
+ else:
+ print test.__doc__
+ print " L'analyse Blue dans l'espace de contrôle est correcte."
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant le fonctionnement du filtre de Kalman sur un système
+ dynamique de trajectoire 1D multiplicative : X(n+1) = G * X(n)
+"""
+__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+import Persistence
+
+#===============================================================================
+def test(dimension = 3):
+ """
+ Cas-test vérifiant le fonctionnement du filtre de Kalman sur un système
+ dynamique de trajectoire 1D multiplicative : X(n+1) = G * X(n)
+ """
+ print test.__doc__
+ #
+ # Définition des données
+ # ----------------------
+ Xt = Persistence.OneScalar("Valeur vraie", basetype=float)
+ gain = 1.01
+ Xt.store(2.5)
+ for i in range(dimension):
+ Xt.store( Xt.valueserie(-1) * gain )
+ Yo = Persistence.OneScalar("Observations", basetype=float)
+ Yo.store(0.)
+ for i in range(dimension):
+ Yo.store(numpy.random.normal(Xt.valueserie(i+1), 0.8, size=(1,)))
+ #
+ # Création de l'étude et résolution
+ # ---------------------------------
+ ADD = AssimilationStudy("Assimilation temporelle de Kalman")
+ #
+ ADD.setBackground (asVector = "0")
+ ADD.setBackgroundError (asCovariance = "1")
+ #
+ ADD.setObservationOperator(asMatrix = "1")
+ ADD.setObservation (asPersistentVector = Yo)
+ ADD.setObservationError (asCovariance = "100")
+ #
+ ADD.setEvolutionModel (asMatrix = [gain,])
+ ADD.setEvolutionError (asCovariance = "1.")
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Kalman")
+ #
+ ADD.analyze()
+ Xa = ADD.get("Analysis")
+ print " Ebauche, i.e. valeur initiale d'analyse...:",Xa.valueserie(0)[0]
+ print " Nombre d'analyses (sans l'ébauche)........:",Xa.stepnumber()-1
+ print " Moyenne des analyses......................:",Xa.stepmean()
+ #
+ # Biais des erreurs
+ EpsY = []
+ for i in range(Yo.stepnumber()):
+ EpsY.append(Yo.valueserie(i) - Xt.valueserie(i))
+ print " Biais des erreurs <Obs-Vraie>.............:",numpy.array(EpsY).mean()
+ print " Variance des erreurs <Obs-Vraie>..........:",numpy.array(EpsY).var()
+ EpsY = []
+ for i in range(Xa.stepnumber()):
+ EpsY.append(Xa.valueserie(i)[0] - Xt.valueserie(i))
+ print " Biais des erreurs <Ana-Vraie>.............:",numpy.array(EpsY).mean()
+ print " Variance des erreurs <Ana-Vraie>..........:",numpy.array(EpsY).var()
+ print
+ #
+ ADD.setDiagnostic("PlotVectors", "Affichage de Xa et Xt")
+ MonPlot = ADD.get("Affichage de Xa et Xt")
+ MonPlot.calculate(
+ ( [ x[0] for x in Xa.valueserie()], Xt.valueserie(), Yo.valueserie() ),
+ title = "Analyse de Kalman sur trajectoire constante",
+ ltitle = ["Analyse", "Valeur vraie", "Observations"],
+ filename = "kalman_sur_trajectoire_multiplicative.ps",
+ pause = False,
+ )
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Analyse moindre carres sans ebauche
+"""
+__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+#===============================================================================
+def test(dimension = 100, precision = 1.e-13):
+ """
+ Analyse moindre carres sans ebauche
+ """
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ H = numpy.identity(dimension)
+ yo = H * xt
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.identity(dimension)
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ # Les valeurs de xb et B ne sont pas utilisées dans l'algorithme
+ # pour lequel on ne considere pas d'ébauche
+ ADD.setBackground (asVector = numpy.zeros((dimension,)) )
+ ADD.setBackgroundError (asCovariance = numpy.zeros((dimension,dimension)) )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ #
+ ADD.setAlgorithm(choice="LinearLeastSquares")
+ #
+ ADD.analyze()
+ #
+ xa = ADD.get("Analysis").valueserie(0)
+ if max(abs(xa - xt.A1)) > precision :
+ raise ValueError("Resultat du test errone")
+ else :
+ print test.__doc__
+ print " Test correct"
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test(3)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que si les covariances d'erreur B et R sont identiques et
+ unitaires, l'analyse est située au milieu de l'ébauche [0,1,2] et de
+ l'observation [0.5,1.5,2.5], avec une erreur d'un ordre inférieur à celle
+ introduite dans R (si l'erreur est de 1 dans R, la précision de vérification
+ est de 0.1*0.1).
+"""
+__author__ = "Jean-Philippe ARGAUD - Novembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+import Persistence
+
+import logging
+# logging.getLogger().setLevel(logging.DEBUG)
+
+#===============================================================================
+def test(precision = 1.e-2):
+ """
+ Cas-test vérifiant que si les covariances d'erreur B et R sont identiques et
+ unitaires, l'analyse est située au milieu de l'ébauche [0,1,2] et de
+ l'observation [0.5,1.5,2.5], avec une erreur d'un ordre inférieur à celle
+ introduite dans R (si l'erreur est de 1 dans R, la précision de vérification
+ est de 0.1*0.1).
+ """
+ #
+ # Définition de l'étude d'assimilation
+ # ------------------------------------
+ ADD = AssimilationStudy("Ma premiere etude")
+ #
+ Xb = Persistence.OneVector("Ebauche", basetype=numpy.matrix)
+ for i in range(100):
+ Xb.store( numpy.matrix( [0,10,20], numpy.float ).T )
+ #
+ ADD.setBackground (asPersistentVector = Xb )
+ ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
+ ADD.setObservation (asVector = [0.5,10.5,20.5])
+ ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
+ ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="EnsembleBlue")
+ #
+ ADD.analyze()
+ #
+ Xa = ADD.get("Analysis")
+ Analyse_moyenne = numpy.matrix( Xa.valueserie() ).mean(axis=0).A1
+ print
+ print " Ebauche :",[0,1,2]
+ print " Analyse moyenne :",Analyse_moyenne
+ print " Nombre d'analyses :",Xa.stepnumber()
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(Analyse_moyenne-numpy.array([0.25, 10.25, 20.25]))/10 > precision:
+ raise ValueError("Résultat du test erroné")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que la Persistence ou le BLUE en boucle donnent bien
+ les résultats attendus.
+"""
+__author__ = "Jean-Philippe ARGAUD - Janvier 2009"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+from Persistence import Persistence
+
+#===============================================================================
+def test(precision = 1.e-13, dimension = 3, nbsteps = 4):
+ """
+ Cas-test vérifiant que la Persistence ou le BLUE en boucle donnent bien
+ les résultats attendus.
+ """
+ vect1 = [1, 2, 1, 2, 1]
+ vect2 = [-3, -3, 0, -3, -3]
+ vect3 = [-1, 1, -5, 1, -1]
+ vect4 = 2*[0.29, 0.97, 0.73, 0.01, 0.20]
+
+ print
+ print " TEST DE LA PERSISTENCE"
+ print " ----------------------"
+ OBJET_DE_TEST = Persistence("My object", unit="", basetype=numpy.array)
+ print " Stockage de 3 vecteurs de longueur identique"
+ OBJET_DE_TEST.store(vect1)
+ OBJET_DE_TEST.store(vect2)
+ OBJET_DE_TEST.store(vect3)
+ print " Stockage d'un quatrième vecteur de longueur différente"
+ OBJET_DE_TEST.store(vect4)
+ print " Taille \"shape\" du dernier objet stocké",OBJET_DE_TEST.shape()
+ print " Taille \"len\" du dernier objet stocké",len(OBJET_DE_TEST)
+
+ print " Affichage des objets avec leur type"
+ for k in range(4):
+ xa = OBJET_DE_TEST.valueserie(k)
+ print " %2i ==> %s, taille %2i, 3ème valeur : %s, objet : %s"%(k,type(xa),len(xa),xa[2],xa)
+
+ del OBJET_DE_TEST
+
+ print
+ print " TEST DE BOUCLE AUTOUR D'UN BLUE"
+ print " -------------------------------"
+ yo = 0.5 + numpy.arange(dimension)
+ B = numpy.matrix(numpy.core.identity(dimension))
+ R = numpy.matrix(numpy.core.identity(dimension))
+ H = numpy.matrix(numpy.core.identity(dimension))
+
+ ADD = AssimilationStudy("Ma premiere etude BLUE")
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ ADD.setAlgorithm(choice="Blue")
+
+ calculs1 = []
+ for i in range(nbsteps):
+ xb = numpy.arange(dimension)
+ xb[min(dimension-1,2)] = i
+ #
+ ADD.setBackground(asVector = xb)
+ ADD.analyze()
+
+ print
+ print " Nombre d'analyses :", ADD.get("Analysis").stepnumber()
+ print " Observation :", yo
+ print " Ebauche :", xb
+ xa = ADD.get("Analysis").valueserie(i)
+ d = ADD.get("Innovation").valueserie(i)
+ print " Analyse résultante :", xa
+ print " so :", float( numpy.dot(d,(yo-numpy.dot(H,xa)).A1) / R.trace() )
+ print " sb :", float( numpy.dot(d,numpy.dot(H,(xa - xb)).A1) /(H * B * H.T).trace() )
+ print " Innovation :", d
+ print " Détails de xa :", type(xa), len(xa), xa[2]
+ calculs1.append(xa[2])
+ del ADD, yo, B, R, H, xb
+
+ print
+ print " TEST DE BOUCLE AUTOUR D'UN BLUE AVEC appliedToX"
+ print " -----------------------------------------------"
+ yo = 0.5 + numpy.arange(dimension)
+ B = numpy.matrix(numpy.core.identity(dimension))
+ R = numpy.matrix(numpy.core.identity(dimension))
+ H = numpy.matrix(numpy.core.identity(dimension))
+
+ ADD = AssimilationStudy("Ma premiere etude BLUE")
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ ADD.setAlgorithm(choice="Blue")
+
+ calculs2 = []
+ for i in range(nbsteps):
+ xb = numpy.arange(dimension)
+ xb[min(dimension-1,2)] = i
+ HXb = numpy.dot(H,xb)
+ #
+ ADD.setObservationOperator(asMatrix = H,
+ appliedToX = {"HXb":HXb})
+ ADD.setBackground(asVector = xb)
+ ADD.analyze()
+
+ print
+ print " Nombre d'analyses :", ADD.get("Analysis").stepnumber()
+ print " Observation :", yo
+ print " Ebauche :", xb
+ print " HXb :", HXb
+ xa = ADD.get("Analysis").valueserie(i)
+ d = ADD.get("Innovation").valueserie(i)
+ print " Analyse résultante :", xa
+ print " so :", float( numpy.dot(d,(yo-numpy.dot(H,xa)).A1) / R.trace() )
+ print " sb :", float( numpy.dot(d,numpy.dot(H,(xa - xb)).A1) /(H * B * H.T).trace() )
+ print " Innovation :", d
+ print " Détails de xa :", type(xa), len(xa), xa[2]
+ calculs2.append(xa[2])
+ del ADD, yo, B, R, H, xb
+
+ #
+ # Vérification du résultat
+ # ------------------------
+ resultats = ( 2.5 + numpy.arange(nbsteps) )/2.
+ calculs1 = numpy.array(calculs1)
+ calculs2 = numpy.array(calculs2)
+ if max(abs(calculs1 - resultats)) > precision:
+ raise ValueError("Résultat du test erroné (1)")
+ elif max(abs(calculs2 - resultats)) > precision:
+ raise ValueError("Résultat du test erroné (2)")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ # numpy.random.seed(1000)
+
+ test(dimension=10)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant sur le 3D-VAR que si l'erreur est nulle, l'incrément
+ d'analyse est nul.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2009"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+import Persistence
+
+import logging
+# Si on désire plus d'information sur le déroulement du calcul, on peut
+# décommenter l'une des lignes qui suit :
+# logging.getLogger().setLevel(logging.INFO)
+# logging.getLogger().setLevel(logging.DEBUG)
+
+#===============================================================================
+def test(precision = 1.e-13, dimension = 3):
+ """
+ Cas-test vérifiant sur le 3D-VAR que si l'erreur est nulle, l'incrément
+ d'analyse est nul.
+ """
+ #
+ # Définition des données
+ # ----------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.zeros((dimension,))).T
+ Eb = numpy.matrix(numpy.zeros((dimension,))).T
+ # Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ # Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ xb = xt + Eb
+ yo = H * xt + Eo
+ #
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Analyse
+ # -------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="3DVAR")
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(xa - xb)) > precision:
+ raise ValueError("Résultat du test erroné (1)")
+ elif max(abs(d)) > precision:
+ raise ValueError("Résultat du test erroné (2)")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 3)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant sur le 3D-VAR que si l'erreur est nulle, l'incrément
+ d'analyse est nul.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2009"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+import Persistence
+
+import logging
+# Si on désire plus d'information sur le déroulement du calcul, on peut
+# décommenter l'une des lignes qui suit :
+# logging.getLogger().setLevel(logging.INFO)
+# logging.getLogger().setLevel(logging.DEBUG)
+
+#===============================================================================
+def test(precision = 1.e-13, dimension = 3):
+ """
+ Cas-test vérifiant sur le 3D-VAR que si l'erreur est nulle, l'incrément
+ d'analyse est nul.
+ """
+ #
+ # Définition des données
+ # ----------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.zeros((dimension,))).T
+ Eb = numpy.matrix(numpy.zeros((dimension,))).T
+ # Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ # Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ # Définition de l'effet de l'opérateur H comme une fonction
+ # ---------------------------------------------------------
+ def FunctionH( X ):
+ return H * X
+ def AdjointH( (X, Y) ):
+ return H.T * Y
+ #
+ xb = xt + Eb
+ yo = FunctionH( xt ) + Eo
+ #
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Analyse
+ # -------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asFunction = {"Tangent":FunctionH,
+ "Adjoint":AdjointH} )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="3DVAR")
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(xa - xb)) > precision:
+ raise ValueError("Résultat du test erroné (1)")
+ elif max(abs(d)) > precision:
+ raise ValueError("Résultat du test erroné (2)")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 3)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant sur le 3D-VAR que si l'erreur d'observation est nulle, on
+ trouve comme analyse la demi-somme de l'ébauche et de la valeur vraie.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2009"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+import Persistence
+
+import logging
+# Si on désire plus d'information sur le déroulement du calcul, on peut
+# décommenter l'une des lignes qui suit :
+# logging.getLogger().setLevel(logging.INFO)
+# logging.getLogger().setLevel(logging.DEBUG)
+
+#===============================================================================
+def test(precision = 1.e-10, dimension = 3):
+ """
+ Cas-test vérifiant sur le 3D-VAR que si l'erreur d'observation est nulle, on
+ trouve comme analyse la demi-somme de l'ébauche et de la valeur vraie.
+ """
+ #
+ # Définition des données
+ # ----------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.zeros((dimension,))).T
+ # Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ # Eb = numpy.matrix(numpy.zeros((dimension,))).T
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ # Définition de l'effet de l'opérateur H comme une fonction
+ # ---------------------------------------------------------
+ def FunctionH( X ):
+ return H * X
+ def AdjointH( (X, Y) ):
+ return H.T * Y
+ #
+ xb = xt + Eb
+ yo = FunctionH( xt ) + Eo
+ #
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Analyse
+ # -------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asFunction = {"Direct":FunctionH,
+ "Tangent":FunctionH,
+ "Adjoint":AdjointH} )
+ #
+ ADD.setAlgorithm(choice="3DVAR")
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(xa - (xb+xt.A1)/2)) > precision:
+ raise ValueError("Résultat du test erroné (1)")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 300)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant sur le 3D-VAR que si l'erreur d'observation est nulle, on
+ trouve comme analyse la demi-somme de l'ébauche et de la valeur vraie.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2009"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+import Persistence
+
+import logging
+# Si on désire plus d'information sur le déroulement du calcul, on peut
+# décommenter l'une des lignes qui suit :
+# logging.getLogger().setLevel(logging.INFO)
+# logging.getLogger().setLevel(logging.DEBUG)
+
+#===============================================================================
+def test(precision = 1.e-10, dimension = 3, minimum = 1.):
+ """
+ Cas-test vérifiant sur le 3D-VAR que si l'erreur d'observation est nulle, on
+ trouve comme analyse la demi-somme de l'ébauche et de la valeur vraie.
+ """
+ #
+ # Définition des données
+ # ----------------------
+ xt = numpy.matrix(numpy.arange(dimension)).T
+ Eo = numpy.matrix(numpy.zeros((dimension,))).T
+ # Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ # Eb = numpy.matrix(numpy.zeros((dimension,))).T
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ # Définition de l'effet de l'opérateur H comme une fonction
+ # ---------------------------------------------------------
+ def FunctionH( X ):
+ return H * X
+ def AdjointH( (X, Y) ):
+ return H.T * Y
+ #
+ xb = xt + Eb
+ yo = FunctionH( xt ) + Eo
+ #
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = numpy.matrix(numpy.core.identity(dimension)).T
+ B = numpy.matrix(numpy.core.identity(dimension)).T
+ #
+ # Definition des bornes
+ # ---------------------
+ Bounds = dimension*[[minimum,None]]
+ #
+ # Analyse
+ # -------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asFunction = {"Direct":FunctionH,
+ "Tangent":FunctionH,
+ "Adjoint":AdjointH} )
+ #
+ ADD.setAlgorithm(choice="3DVAR")
+ ADD.setAlgorithmParameters(asDico={
+ "Minimizer":"LBFGSB",
+ "Bounds" :Bounds,
+ })
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ #
+ # Vérification du résultat
+ # ------------------------
+ if xa.min() < minimum:
+ raise ValueError("Résultat du test erroné (1)")
+ else:
+ print test.__doc__
+ print " Test correct, valeur minimale de %s respectée"%minimum
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 300, minimum = 2.)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant le calcul de RMS.
+"""
+__author__ = "Jean-Philippe ARGAUD - Juillet 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(precision = 1.e-13):
+ """
+ Cas-test vérifiant des calculs de RMS.
+ """
+ #
+ ADD = AssimilationStudy()
+ #
+ ADD.setDiagnostic("RMS", "Calcul de RMS multi-pas")
+ #
+ # La ligne suivante permet de simplifier les écritures ultérieures pour
+ # les "calculate", mais n'est pas indispensable : on aurait pu conserver à
+ # chaque appel la commande "ADD.get("...")"
+ #
+ RMS = ADD.get("Calcul de RMS multi-pas")
+ #
+ vect1 = [1, 2, 1, 2, 1]
+ vect2 = [2, 1, 2, 1, 2]
+ RMS.calculate(vect1,vect2)
+ vect1 = [1, 3, 1, 3, 1]
+ vect2 = [2, 2, 2, 2, 2]
+ RMS.calculate(vect1,vect2)
+ vect1 = [1, 1, 1, 1, 1]
+ vect2 = [2, 2, 2, 2, 2]
+ RMS.calculate(vect1,vect2)
+ vect1 = [1, 1, 1, 1, 1]
+ vect2 = [4, -2, 4, -2, -2]
+ RMS.calculate(vect1,vect2)
+ vect1 = [0.29, 0.97, 0.73, 0.01, 0.20]
+ vect2 = [0.92, 0.86, 0.11, 0.72, 0.54]
+ RMS.calculate(vect1,vect2)
+ vect1 = [-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]
+ vect2 = [0,0,0,0,0,0,0,0,0,0]
+ RMS.calculate(vect1,vect2)
+ #
+ Valeurs_attendues = [1.0, 1.0, 1.0, 3.0, 0.53162016515553656, 0.73784217096601323]
+ #
+ # Vérification du résultat
+ # ------------------------
+ ecart = abs( max( numpy.array(RMS.valueserie()) - numpy.array(Valeurs_attendues) ) )
+ if ecart > precision:
+ raise "Résultat du test erroné"
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ # numpy.random.seed(1000)
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant l'affichage multi-pas Gnuplot d'un vecteur.
+"""
+__author__ = "Jean-Philippe ARGAUD - Juillet 2008"
+
+execfile("context.py")
+
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(dimension = 100):
+ """
+ Cas-test vérifiant l'affichage multi-pas Gnuplot d'un vecteur.
+ """
+ #
+ ADD = AssimilationStudy()
+ #
+ ADD.setDiagnostic("PlotVector", "Affichage multi-pas Gnuplot d'un vecteur")
+ #
+ MonPlot = ADD.get("Affichage multi-pas Gnuplot d'un vecteur")
+ #
+ vect = [1, 2, 1, 2, 1]
+ MonPlot.calculate(vect, title = "Vecteur 1", xlabel = "Axe X", ylabel = "Axe Y", pause = False )
+ vect = [1, 3, 1, 3, 1]
+ MonPlot.calculate(vect, title = "Vecteur 2", filename = "vecteur.ps", pause = False)
+ vect = [-1, 1, 1, 1, -1]
+ MonPlot.calculate(vect, title = "Vecteur 3", pause = False)
+ vect = [0.29, 0.97, 0.73, 0.01, 0.20]
+ MonPlot.calculate(vect, title = "Vecteur 4", pause = False)
+ vect = [-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]
+ MonPlot.calculate(vect, title = "Vecteur 5", pause = False)
+ vect = dimension*[0.29, 0.97, 0.73, 0.01, 0.20]
+ MonPlot.calculate(vect, title = "Vecteur 6 : long construit par repetition", pause = False)
+ vect = [0.29, 0.97, 0.73, 0.01, 0.20]
+ MonPlot.calculate(vect, title = "Vecteur 7", pause = False)
+ temps = [0.1,0.2,0.3,0.4,0.5]
+ MonPlot.calculate(vect, temps, title = "Vecteur 8 avec axe du temps modifie", pause = False)
+ #
+ # Vérification du résultat
+ # ------------------------
+ print test.__doc__
+ print " Test correct"
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant l'affichage multi-pas Gnuplot d'une liste de vecteurs.
+"""
+__author__ = "Jean-Philippe ARGAUD - Juillet 2008"
+
+execfile("context.py")
+
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(dimension = 100):
+ """
+ Cas-test vérifiant l'affichage multi-pas Gnuplot d'une liste de vecteurs.
+ """
+ #
+ ADD = AssimilationStudy()
+ #
+ ADD.setDiagnostic("PlotVectors", "Affichage multi-pas Gnuplot d'une liste de vecteurs")
+ #
+ MonPlot = ADD.get("Affichage multi-pas Gnuplot d'une liste de vecteurs")
+ #
+ vect1 = [1, 2, 1, 2, 1]
+ MonPlot.calculate([vect1,], title = "Vecteur 1", xlabel = "Axe X", ylabel = "Axe Y", pause = False )
+ vect2 = [1, 3, 1, 3, 1]
+ MonPlot.calculate([vect1,vect2], title = "Vecteurs 1 et 2", filename = "liste_de_vecteurs.ps", pause = False )
+ vect3 = [-1, 1, -1, 1, -1]
+ MonPlot.calculate((vect1,vect2,vect3), title = "Vecteurs 1 a 3", pause = False )
+ vect4 = 100*[0.29, 0.97, 0.73, 0.01, 0.20]
+ MonPlot.calculate([vect4,], title = "Vecteur 4 : long construit par repetition", pause = False )
+ MonPlot.calculate(
+ (vect1,vect2,vect3),
+ [0.1,0.2,0.3,0.4,0.5],
+ title = "Vecteurs 1 a 3, temps modifie", pause = False)
+ print
+ #
+ # Vérification du résultat
+ # ------------------------
+ print test.__doc__
+ print " Test correct"
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-\r
+__doc__ = """\r
+ Test d'homogenéité des distributions de OMB et OMA lors d'une analyse Blue\r
+"""\r
+__author__ = "Sophie RICCI - Septembre 2008"\r
+\r
+execfile("context.py")\r
+\r
+import numpy\r
+from AssimilationStudy import AssimilationStudy\r
+\r
+#===============================================================================\r
+def test(dimension = 75):\r
+ """\r
+ Test d'homogenéité des distributions de OMB et OMA lors d'une analyse Blue\r
+ """\r
+ numpy.random.seed(1000)\r
+ #\r
+ # Définition des données "théoriques" vraies\r
+ # ------------------------------------------\r
+ xt = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T\r
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T\r
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T\r
+ #\r
+ H = numpy.matrix(numpy.core.identity(dimension))\r
+ #\r
+ xb = xt + Eb\r
+ yo = H * xt + Eo\r
+ #\r
+ xb = xb.A1\r
+ yo = yo.A1\r
+ #\r
+ # Définition des matrices de covariances d'erreurs\r
+ # ------------------------------------------------\r
+ R = numpy.matrix(numpy.core.identity(dimension)).T\r
+ B = numpy.matrix(numpy.core.identity(dimension)).T\r
+ #\r
+ # Analyse BLUE\r
+ # ------------\r
+ ADD = AssimilationStudy()\r
+ ADD.setBackground (asVector = xb )\r
+ ADD.setBackgroundError (asCovariance = B )\r
+ ADD.setObservation (asVector = yo )\r
+ ADD.setObservationError (asCovariance = R )\r
+ ADD.setObservationOperator(asMatrix = H )\r
+ #\r
+ ADD.setControls()\r
+ ADD.setAlgorithm(choice="Blue")\r
+ #\r
+ ADD.analyze()\r
+ #\r
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))\r
+ d = numpy.array(ADD.get("Innovation").valueserie(0))\r
+ OMA = yo - xa\r
+ #\r
+ # Application du test d'adéquation du Khi-2 \r
+ # -----------------------------------------\r
+ ADD.setDiagnostic("HomogeneiteKhi2",\r
+ name = "Test d'homogeneite entre OMB et OMA par calcul du khi2",\r
+ parameters = { "tolerance":0.05, "nbclasses":8 , "dxclasse":None})\r
+ #\r
+ # Instanciation de l'objet testkhi2\r
+ # ---------------------------------\r
+ D = ADD.get("Test d'homogeneite entre OMB et OMA par calcul du khi2")\r
+ #\r
+ # Calcul du test et résultat\r
+ # --------------------------\r
+ D.calculate(d, OMA)\r
+ print " Reponse du test", D.valueserie()\r
+ print\r
+\r
+#==============================================================================\r
+if __name__ == "__main__":\r
+\r
+ print\r
+ print "AUTODIAGNOSTIC"\r
+ print "=============="\r
+ \r
+ test()\r
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Diagnostic sur les variances dans B et R par rapport à l'ébauche Xb et aux
+ observations Y. On teste si on a les conditions :
+ 1%*xb < sigma_b < 10%*xb
+ et
+ 1%*yo < sigma_o < 10%*yo
+ lors d une anlayse BLUE.
+"""
+__author__ = "Sophie RICCI - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(dimension = 3):
+ """
+ Diagnostic sur les variances dans B et R par rapport à l'ébauche Xb et aux
+ observations Y. On teste si on a les conditions :
+ 1%*xb < sigma_b < 10%*xb
+ et
+ 1%*yo < sigma_o < 10%*yo
+ lors d une anlayse BLUE.
+ """
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ #
+ H = numpy.matrix(numpy.core.identity(dimension))
+ #
+ xb = xt + Eb
+ yo = H * xt + Eo
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = 1.e-3 * yo.mean() * yo.mean() * numpy.matrix(numpy.core.identity(dimension))
+ B = 1.e-3 * xb.mean() * xb.mean() * numpy.matrix(numpy.core.identity(dimension))
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ #
+ # Application du test
+ # -------------------
+ ADD.setDiagnostic("VarianceOrder", name = "Ordre des matrices de covariance")
+ #
+ D = ADD.get("Ordre des matrices de covariance")
+ #
+ D.calculate( Xb = xb, B = B, Y = yo, R = R )
+ #
+ # Verification du resultat
+ # ------------------------
+ if not D.valueserie(0) :
+ raise ValueError("Resultat du test errone ")
+ else :
+ print test.__doc__
+ print " Test correct"
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Test d'égalité des variances de OMB et OMA lors d'une analyse Blue
+ au sens du test de Fisher.
+"""
+__author__ = "Sophie RICCI - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(dimension = 500):
+ """
+ Test d'égalité des variances de OMB et OMA lors d'une analyse Blue
+ au sens du test de Fisher.
+ """
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension,))).T
+ H = numpy.matrix(numpy.identity(dimension))
+ xb = xt + Eb
+ yo = H * xt + Eo
+ xb = xb.A1
+ yo = yo.A1
+ #
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ R = 1000. * numpy.matrix(numpy.identity(dimension)).T
+ B = numpy.matrix(numpy.identity(dimension)).T
+ #
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ ADD.analyze()
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))
+ d = numpy.array(ADD.get("Innovation").valueserie(0))
+ OMA = yo - xa
+ #
+ # Application du test d adequation du Khi-2
+ # -----------------------------------------
+ ADD.setDiagnostic("CompareVarianceFisher",
+ name = "Test de comparaison des variances de OMB et OMA par test de Fisher",
+ parameters = { "tolerance":0.05 })
+ #
+ # Instanciation du diagnostic
+ # ---------------------------
+ D = ADD.get("Test de comparaison des variances de OMB et OMA par test de Fisher")
+ #
+ # Calcul
+ # ------
+ D.calculate(d, OMA)
+ if not D.valueserie(0) :
+ raise ValueError("L'analyse ne change pas de manière significative la variance. Le test est erroné.")
+ else :
+ print test.__doc__
+ print " L'analyse effectuée change de manière significative la variance."
+ print " Test correct"
+ print
+
+#==============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-\r
+__doc__ = """\r
+ Test d adequation des distributions de OMB et OMA avec une distribution \r
+ gaussienne dont la moyenne et la std sont calculees sur l echantillon.\r
+ L analyse est un Blue.\r
+"""\r
+__author__ = "Sophie RICCI - Septembre 2008"\r
+\r
+execfile("context.py")\r
+\r
+import numpy\r
+from AssimilationStudy import AssimilationStudy\r
+\r
+#===============================================================================\r
+def test(dimension = 10):\r
+ """\r
+ Test d adequation des distributions de OMB et OMA avec une distribution \r
+ gaussienne dont la moyenne et la std sont calculees sur l echantillon.\r
+ L analyse est un Blue.\r
+ """\r
+ #\r
+ # Définition des données "théoriques" vraies\r
+ # ------------------------------------------\r
+ xt = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension))).T\r
+ Eo = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension))).T\r
+ Eb = numpy.matrix(numpy.random.normal(0.,1.,size=(dimension))).T\r
+ H = numpy.matrix(numpy.identity(dimension))\r
+ xb = xt + Eb\r
+ yo = H * xt + Eo\r
+ xb = xb.A1\r
+ yo = yo.A1\r
+ # Définition des matrices de covariances d'erreurs\r
+ # ------------------------------------------------\r
+ R = numpy.matrix(numpy.identity(dimension)).T\r
+ B = numpy.matrix(numpy.identity(dimension)).T\r
+ #\r
+ # Analyse BLUE\r
+ # ------------\r
+ ADD = AssimilationStudy()\r
+ ADD.setBackground (asVector = xb )\r
+ ADD.setBackgroundError (asCovariance = B )\r
+ ADD.setObservation (asVector = yo )\r
+ ADD.setObservationError (asCovariance = R )\r
+ ADD.setObservationOperator(asMatrix = H )\r
+ #\r
+ ADD.setControls()\r
+ ADD.setAlgorithm(choice="Blue")\r
+ #\r
+ ADD.analyze()\r
+ #\r
+ xa = numpy.array(ADD.get("Analysis").valueserie(0))\r
+ d = numpy.array(ADD.get("Innovation").valueserie(0))\r
+ OMA = yo - xa\r
+ \r
+ #\r
+ # Application du test d adequation du Khi-2 \r
+ # -------------------------------------------------\r
+ ADD.setDiagnostic("GaussianAdequation",\r
+ name = "Test d adequation a une gaussienne par calcul du khi2",\r
+ parameters = { "tolerance":0.05, "nbclasses":8., "dxclasse":None })\r
+ #\r
+ # Instanciation de l'objet testkhi2\r
+ # --------------------------------------------------------------------\r
+ D = ADD.get("Test d adequation a une gaussienne par calcul du khi2")\r
+ #\r
+ # Calcul\r
+ # --------------------------------------------------------------------\r
+ print test.__doc__\r
+ D.calculate(d)\r
+ if not D.valueserie(0) :\r
+ raise ValueError("L'adéquation à une gaussienne pour la variable OMB n'est pasvalide.")\r
+ else :\r
+ print " L'adéquation à une gaussienne pour la variable OMB est valide."\r
+ D.calculate(OMA)\r
+ if not D.valueserie(1) :\r
+ raise ValueError("L'adéquation à une gaussienne pour la variable OMA n'est pasvalide.")\r
+ else :\r
+ print " L'adéquation à une gaussienne pour la variable OMA est valide."\r
+ print\r
+\r
+#==============================================================================\r
+if __name__ == "__main__":\r
+\r
+ print\r
+ print "AUTODIAGNOSTIC"\r
+ print "=============="\r
+ \r
+ test()\r
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Recherche de l'arrêt de la réduction de la variance (VAR(OMB)-VAR(OMA))
+ lors d'itérations sur une analyse Blue avec R = alpha*B et H = Id.
+ - avec remise à jour de l'ébauche xb = xa (updatexb = True)
+ - avec correction de R et B par so et sb (sosb = True)
+"""
+__author__ = "Sophie RICCI - Septembre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+from scipy import asarray
+
+#===============================================================================
+def test(dimension = 3, alpha = 1., N = 10, updatexb = True, sosb = False) :
+ #
+ # Définition des données "théoriques" vraies
+ # ------------------------------------------
+ xt = numpy.arange(dimension)
+ Eo = numpy.random.normal(0.,1.,size=(dimension,))
+ Eb = numpy.zeros((dimension,))
+ H = numpy.identity(dimension)
+ xb = xt + Eb
+ yo = numpy.dot(H,xt) + Eo
+ # Définition des matrices de covariances d'erreurs
+ # ------------------------------------------------
+ B = numpy.identity(dimension)
+ R = alpha * B
+ # Analyse BLUE
+ # ------------
+ ADD = AssimilationStudy()
+ ADD.setObservation (asVector = yo )
+ ADD.setObservationOperator(asMatrix = H )
+ ADD.setControls()
+ ADD.setAlgorithm(choice="Blue")
+ SigmaBck2 = 1.
+ SigmaObs2 = 1.
+ VectSigmaObs2, VectSigmaBck2 = [],[]
+ vectd , vectOMA = [],[]
+
+ # Iterations du Blue
+ for i in range (0, N) :
+ ADD.setBackground (asVector = xb )
+ # Mise a jour de R et B par so et sb
+ if sosb :
+ newB = SigmaBck2*B
+ newR = SigmaObs2*R
+ ADD.setBackgroundError (asCovariance = newB )
+ ADD.setObservationError(asCovariance = newR )
+ else :
+ newB = B
+ newR = R
+ ADD.setBackgroundError (asCovariance = newB )
+ ADD.setObservationError(asCovariance = newR )
+ ADD.analyze()
+ xa = ADD.get("Analysis").valueserie(i)
+ d = ADD.get("Innovation").valueserie(i)
+ # Construit le vecteur des so et sb
+ SigmaObs2 = ADD.get("SigmaObs2").valueserie(i)
+ SigmaBck2 = ADD.get("SigmaBck2").valueserie(i)
+ VectSigmaObs2.append(SigmaObs2)
+ VectSigmaBck2.append(SigmaBck2)
+
+ # Calcul de la variance de OMB et OMA
+ OMB = yo -xb
+ var_OMB = OMB.var()
+ vectd.append(var_OMB)
+
+ OMA = yo-xa
+ var_OMA = OMA.var()
+ vectOMA.append(var_OMA)
+
+ # Update de l ebauche par l analyse
+ if updatexb :
+ xb = xa
+
+ # Construction du vecteur de difference VAR(OMB)-VAR(0MA)
+ vectd = asarray(vectd)
+ vectOMA = asarray(vectOMA)
+ vector = asarray(vectd) - asarray(vectOMA)
+
+ # Plot de VAR(d) - VAR(OMA) au cours des iterations
+ # --------------------------------------------------
+ ADD.setDiagnostic("PlotVector", "Affichage multi-pas Gnuplot d'un vecteur")
+ MonPlot = ADD.get("Affichage multi-pas Gnuplot d'un vecteur")
+# MonPlot.calculate(vector, title = " VAR(d) - VAR(OMA) ", xlabel = "Axe X", ylabel = "Axe Y", filename = "Plot_StopReductionVariance_VAROMB-VAROMA.ps", pause = True)
+# MonPlot.calculate(VectSigmaObs2, title = " SigmaObs2 ", xlabel = "Axe X", ylabel = "Axe Y", filename = "testiter_so.ps", pause = True)
+# MonPlot.calculate(VectSigmaBck2, title = " SigmaBck2 ", xlabel = "Axe X", ylabel = "Axe Y", filename = "testiter_sb.ps", pause = True)
+
+
+ # Application du diagnostic sur l arret de la reduction de la variance
+ # ----------------------------------------------------------------------
+ ADD.setDiagnostic("StopReductionVariance",
+ name = "Arret de la reduction de la variance entre OMB et OMA")
+ #
+ D = ADD.get("Arret de la reduction de la variance entre OMB et OMA")
+ D.calculate( vector = vector, CutOffSlope = 0.005, MultiSlope0 = None)
+
+ # Verification du resultat
+ # ------------------------
+ print __doc__
+ print " La variance n'est plus significativement réduite après l'itération", D.valueserie(0)
+ print " Test correct"
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test()
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Diagnotic de test sur la validité de l'hypothèse de linéarité de l'opérateur
+ H entre xp et xm
+"""
+__author__ = "Jean-Philippe ARGAUD - Octobre 2008"
+
+execfile("context.py")
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test(tolerance = 0.1, dimension = 3):
+ """
+ Diagnotic de test sur la validité de l'hypothèse de linéarité de l'opérateur
+ H entre xp et xm
+ """
+ #
+ # Définition des données
+ # ----------------------
+ dxparam = 1.
+ Hxm = numpy.random.normal(0.,1.,size=(dimension,))
+ Hxp = Hxm + 2*dxparam
+ Hx = (Hxp + Hxm)/2.
+ H = (Hxp - Hxm)/2.
+ #
+ # Instanciation de l'objet diagnostic
+ # -----------------------------------
+ ADD = AssimilationStudy()
+ ADD.setDiagnostic("HLinearity",
+ name = "Test le linearite de Hlin",
+ parameters = { "tolerance":tolerance })
+ D = ADD.get("Test le linearite de Hlin")
+ #
+ # Calcul
+ # ------
+ D.calculate(
+ Hlin = H,
+ deltaparam = dxparam,
+ Hxp = Hxp,
+ Hxm = Hxm,
+ Hx = Hx)
+ #
+ # Vérification du résultat
+ # ------------------------
+ if not D.valueserie(0) :
+ raise ValueError("Résultat du test erroné")
+ else:
+ print test.__doc__
+ print " Test correct, tolerance du test fixée à %s"%tolerance
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 1.e4) # Fonctionne bien jusqu'à 1.e7
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que la minimisation CG de Scipy fonctionne
+ sur un cas simple.
+"""
+__author__ = "Jean-Philippe ARGAUD - Avril 2009"
+
+import numpy
+import scipy.optimize
+
+import logging
+# Si on désire plus d'information sur le déroulement du calcul, on peut
+# décommenter l'une des lignes qui suit :
+# logging.getLogger().setLevel(logging.INFO)
+# logging.getLogger().setLevel(logging.DEBUG)
+
+#===============================================================================
+class CostFunction:
+ """
+ Classe permettant de rassembler toutes les informations disponibles sur la
+ fonction-coût nécessaires à la minimisation.
+
+ Il y a 2 méthodes :
+ - value : renvoie la valeur de la fonction-coût, i.e. J
+ - gradient : renvoie son gradient, i.e. grad(J)
+
+ La fonction-coût choisie est une simple norme quadratique sur l'écart entre
+ la variable de minimisation X et une valeur constante X0. Le gradient est
+ donc de deux fois cet écart, et le minimum est atteint quand X=X0.
+ """
+ def __init__(self, X0 = None ):
+ self.X0 = X0
+ self.i = 0
+ logging.debug("")
+ logging.debug("Initialisations pour J :")
+ logging.debug(" X0 = %s"%self.X0)
+
+ def value(self, X = None ):
+ #
+ self.i += 1
+ #
+ J = numpy.dot( ( X - self.X0 ), ( X - self.X0 ) )
+ J = float( J )
+ #
+ logging.debug("")
+ logging.debug("Etape de minimisation numéro %i"%self.i)
+ logging.debug("------------------------------")
+ logging.debug("Calcul de la valeur de J :")
+ logging.debug(" X = %s"%X)
+ logging.debug(" J = %s"%J)
+ return J
+
+ def gradient(self, X = None ):
+ #
+ gradJ = 2. * ( X - self.X0 )
+ #
+ logging.debug("Calcul du gradient de J :")
+ logging.debug(" grad(J) = %s"%gradJ)
+ return gradJ
+
+ def iterations(self):
+ return self.i
+
+#===============================================================================
+def test(precision = 1.e-07, dimension = 3):
+ """
+ Cas-test vérifiant que la minimisation CG de Scipy fonctionne
+ sur un cas simple.
+ """
+ #
+ # Définition de l'objet contenant la fonction-coût
+ # ------------------------------------------------
+ X0 = numpy.random.normal(0.,1.,size=(dimension,))
+ J = CostFunction( X0 )
+ #
+ X_initial = 3. * X0
+ #
+ # X_optimal, J_optimal, Informations
+ X_optimal, J_optimal, func_calls, grad_calls, warnflag = scipy.optimize.fmin_cg(
+ f = J.value,
+ x0 = X_initial,
+ fprime = J.gradient,
+ args = (),
+ gtol = precision,
+ full_output = True,
+ )
+ #
+ GradJ_opt = J.gradient( X_optimal )
+ #
+ logging.info("")
+ logging.info("Résultats finaux :")
+ logging.info(" X0 = %s"%X0)
+ logging.info(" X_optimal = %s"%X_optimal)
+ logging.info(" J_optimal = %s"%J_optimal)
+ logging.info(" GradJ_opt = %s"%GradJ_opt)
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(GradJ_opt)) > precision:
+ raise ValueError("Résultat du test erroné sur le gradient de J")
+ elif J_optimal > precision:
+ raise ValueError("Résultat du test erroné sur J")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print " Nombre de calculs de la fonctionnelle J : %i"%J.iterations()
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que la minimisation BFGS de Scipy fonctionne
+ sur un cas simple.
+"""
+__author__ = "Jean-Philippe ARGAUD - Avril 2009"
+
+import numpy
+import scipy.optimize
+
+import logging
+# Si on désire plus d'information sur le déroulement du calcul, on peut
+# décommenter l'une des lignes qui suit :
+# logging.getLogger().setLevel(logging.INFO)
+# logging.getLogger().setLevel(logging.DEBUG)
+
+from test300_Optimize_CG import CostFunction
+
+#===============================================================================
+def test(precision = 1.e-07, dimension = 3):
+ """
+ Cas-test vérifiant que la minimisation BFGS de Scipy fonctionne
+ sur un cas simple.
+ """
+ #
+ # Définition de l'objet contenant la fonction-coût
+ # ------------------------------------------------
+ X0 = numpy.random.normal(0.,1.,size=(dimension,))
+ J = CostFunction( X0 )
+ #
+ X_initial = 3. * X0
+ #
+ X_optimal, J_optimal, GradJ_opt, Hopt, func_calls, grad_calls, warnflag = scipy.optimize.fmin_bfgs(
+ f = J.value,
+ x0 = X_initial,
+ fprime = J.gradient,
+ args = (),
+ gtol = precision,
+ full_output = True,
+ )
+ #
+ logging.info("")
+ logging.info("Résultats finaux :")
+ logging.info(" X0 = %s"%X0)
+ logging.info(" X_optimal = %s"%X_optimal)
+ logging.info(" J_optimal = %s"%J_optimal)
+ logging.info(" GradJ_opt = %s"%GradJ_opt)
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(GradJ_opt)) > precision:
+ raise ValueError("Résultat du test erroné sur le gradient de J")
+ elif J_optimal > precision:
+ raise ValueError("Résultat du test erroné sur J")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print " Nombre de calculs de la fonctionnelle J : %i"%J.iterations()
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que la minimisation L-BFGS-B de Scipy fonctionne
+ sur un cas simple.
+"""
+__author__ = "Jean-Philippe ARGAUD - Avril 2009"
+
+import numpy
+import scipy.optimize
+
+import logging
+# Si on désire plus d'information sur le déroulement du calcul, on peut
+# décommenter l'une des lignes qui suit :
+# logging.getLogger().setLevel(logging.INFO)
+# logging.getLogger().setLevel(logging.DEBUG)
+
+from test300_Optimize_CG import CostFunction
+
+#===============================================================================
+def test(precision = 1.e-07, dimension = 3):
+ """
+ Cas-test vérifiant que la minimisation L-BFGS-B de Scipy fonctionne
+ sur un cas simple.
+ """
+ #
+ # Définition de l'objet contenant la fonction-coût
+ # ------------------------------------------------
+ X0 = numpy.random.normal(0.,1.,size=(dimension,))
+ J = CostFunction( X0 )
+ #
+ X_initial = 3. * X0
+ #
+ X_optimal, J_optimal, Informations = scipy.optimize.fmin_l_bfgs_b(
+ func = J.value,
+ x0 = X_initial,
+ fprime = J.gradient,
+ args = (),
+ pgtol = precision,
+ )
+ #
+ logging.info("")
+ logging.info("Résultats finaux :")
+ logging.info(" X0 = %s"%X0)
+ logging.info(" X_optimal = %s"%X_optimal)
+ logging.info(" J_optimal = %s"%J_optimal)
+ logging.info(" GradJ_opt = %s"%Informations['grad'])
+ #
+ # Vérification du résultat
+ # ------------------------
+ if max(abs(Informations['grad'])) > precision:
+ raise ValueError("Résultat du test erroné sur le gradient de J")
+ elif J_optimal > precision:
+ raise ValueError("Résultat du test erroné sur J")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print " Nombre de calculs de la fonctionnelle J : %i"%J.iterations()
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant que la minimisation TNC de Scipy fonctionne
+ sur un cas simple.
+"""
+__author__ = "Jean-Philippe ARGAUD - Avril 2009"
+
+import numpy
+import scipy.optimize
+
+import logging
+# Si on désire plus d'information sur le déroulement du calcul, on peut
+# décommenter l'une des lignes qui suit :
+# logging.getLogger().setLevel(logging.INFO)
+# logging.getLogger().setLevel(logging.DEBUG)
+
+if logging.getLogger().level < 30:
+ message = scipy.optimize.tnc.MSG_ALL
+else:
+ message = scipy.optimize.tnc.MSG_NONE
+
+from test300_Optimize_CG import CostFunction
+
+#===============================================================================
+def test(precision = 1.e-07, dimension = 3):
+ """
+ Cas-test vérifiant que la minimisation TNC de Scipy fonctionne
+ sur un cas simple.
+ """
+ #
+ # Définition de l'objet contenant la fonction-coût
+ # ------------------------------------------------
+ X0 = numpy.random.normal(0.,1.,size=(dimension,))
+ J = CostFunction( X0 )
+ #
+ X_initial = 3. * X0
+ #
+ X_optimal, neval, rc = scipy.optimize.fmin_tnc(
+ func = J.value,
+ x0 = X_initial,
+ fprime = J.gradient,
+ args = (),
+ approx_grad = False,
+ pgtol = precision,
+ messages = message,
+ )
+ #
+ J_optimal = J.value( X_optimal )
+ GradJ_opt = J.gradient( X_optimal )
+ #
+ logging.info("")
+ logging.info("Résultats finaux :")
+ logging.info(" X0 = %s"%X0)
+ logging.info(" X_optimal = %s"%X_optimal)
+ logging.info(" J_optimal = %s"%J_optimal)
+ logging.info(" GradJ_opt = %s"%GradJ_opt)
+ #
+ # Vérification du résultat
+ # ------------------------
+ if J_optimal > precision:
+ raise ValueError("Résultat du test erroné sur J")
+ else:
+ print test.__doc__
+ print " Test correct, erreur maximale inférieure à %s"%precision
+ print " Nombre de calculs de la fonctionnelle J : %i"%J.iterations()
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+ numpy.random.seed(1000)
+
+ test(dimension = 100)
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test vérifiant les messages de sortie pour toutes les informations
+ dynamiques
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+execfile("context.py")
+
+from AssimilationStudy import AssimilationStudy
+
+#===============================================================================
+def test():
+ """
+ Cas-test vérifiant les messages de sortie pour toutes les informations
+ dynamiques
+ """
+ #
+ # Définition de l'étude d'assimilation
+ # ------------------------------------
+ ADD = AssimilationStudy("Verifications des informations dynamiques")
+ #
+ print test.__doc__
+ # print "Chemin des algorithmes :",ADD.get_algorithms_main_path()
+ print "Algorithmes disponibles :",ADD.get_available_algorithms()
+ print
+ # print "Chemin des diagnostics :",ADD.get_diagnostics_main_path()
+ print "Diagnostics disponibles :",ADD.get_available_diagnostics()
+ print
+ chemin = "../../Sources"
+ print "Ajout du chemin :",chemin
+ ADD.add_algorithms_path(chemin)
+ print "Algorithmes disponibles :",ADD.get_available_algorithms()
+ print
+ print "Ajout du chemin :",chemin
+ ADD.add_diagnostics_path(chemin)
+ print "Diagnostics disponibles :",ADD.get_available_diagnostics()
+ print
+ ADD.setAlgorithm(choice=ADD.get_available_algorithms()[0])
+ liste = ADD.get().keys()
+ liste.sort()
+ print "Variables et diagnostics disponibles :",liste
+ print
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ test()
--- /dev/null
+#-*-coding:iso-8859-1-*-
+__doc__ = """
+ Cas-test démontrant les possibilités de logging, sous la forme de fonctions
+ successives illustrant le fonctionnement par défaut.
+"""
+__author__ = "Jean-Philippe ARGAUD - Octotbre 2008"
+
+import os
+execfile("context.py")
+
+#===============================================================================
+def demo1():
+ from AssimilationStudy import AssimilationStudy
+ print """
+ DEMO 1
+ ------
+ L'initialisation de l'environnement de logging a été automatiquement faite à
+ l'import de AssimilationStudy.
+
+ Seuls les messages d'un niveau au moins égal à warning sont disponibles par
+ défaut. Cela permet de truffer le code de messages de DEBUG ou d'INFO sans
+ qu'ils apparaissent à l'affichage standard.
+"""
+ import logging
+ logging.debug("Un message de debug")
+ logging.info("Un message d'info")
+ logging.warning("Un message d'avertissement")
+ logging.error("Un message d'erreur")
+ logging.critical("Un message critique")
+
+def demo2():
+ from AssimilationStudy import AssimilationStudy
+ print """
+ DEMO 2
+ ------
+ On recommence le cas précédent, mais en affichant tous les messages. Cela
+ permet de deboguer le code en ayant les messages non affichés précédemment.
+
+ La commande de disponibilité de tous les niveaux est atteignable à travers le
+ module standard "logging" (avec une minuscule) :
+ logging.getLogger().setLevel(...)
+"""
+ import logging
+
+ logging.getLogger().setLevel(logging.DEBUG)
+
+ logging.debug("Un message de debug")
+ logging.info("Un message d'info")
+ logging.warning("Un message d'avertissement")
+ logging.error("Un message d'erreur")
+ logging.critical("Un message critique")
+
+def demo3():
+ print """
+ DEMO 3
+ ------
+ On peut disposer des messages conjointement à la console et dans un fichier.
+
+ Pour cela, il faut importer le module Logging n'importe où (après le module
+ AssimilationStudy ou en son absence, mais pas avant). On en profite aussi pour
+ initier le logging général avec le niveau INFO, donc le message de debug
+ précédemment affiché ne l'est plus.
+
+"""
+ import Logging
+ Logging.Logging().setLogfile()
+
+ if os.path.isfile("AssimilationStudy.log"):
+ print " ---> Le fichier attendu a bien été créé."
+ else:
+ raise ValueError("Le fichier attendu n'a pas été créé.")
+
+ import logging
+ logging.getLogger().setLevel(logging.INFO)
+
+ logging.debug("Un message de debug")
+ logging.info("Un message d'info")
+ logging.warning("Un message d'avertissement")
+ logging.error("Un message d'erreur")
+ logging.critical("Un message critique")
+
+ print
+ print " On peut vérifier le contenu du fichier \"AssimilationStudy.log\"."
+
+#===============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ demo1()
+ demo2()
+ demo3()
+ print
+ print " Pour les autres modes avancés de contôle du fichier et des niveaux"
+ print " on se reportera à la documentation du module \"Logging\"."
+ print
DATA_INST = \
test000_Blue_AnalysisCode.py test000_Blue_AnalysisFile.py test000_Blue.py
-tests_DATA = ${DATA_INST}
+testsdasalome_DATA = ${DATA_INST}
+
EXTRA_DIST = test000_Blue_AnalysisCode.py test000_Blue_AnalysisFile.py.in test000_Blue.py
