+++ /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
-
- print "xb", xb
- print "B", B
- print "yo", yo
- print "R", R
-
- #
- # 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
