src/tests/daSalome/test017_3DVAR_par_fonction.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
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD - Mars 2009"
-import sys ; sys.path.insert(0, "../daCore")
import logging
-from daCore import Persistence
-from daCore.BasicObjects import Algorithm
-from daCore import PlatformInfo ; m = PlatformInfo.SystemUsage()
+from daCore import BasicObjects, PlatformInfo
+m = PlatformInfo.SystemUsage()
import numpy
import scipy.optimize
disp = 0
# ==============================================================================
-class ElementaryAlgorithm(Algorithm):
+class ElementaryAlgorithm(BasicObjects.Algorithm):
def __init__(self):
- Algorithm.__init__(self)
+ BasicObjects.Algorithm.__init__(self)
self._name = "3DVAR"
logging.debug("%s Initialisation"%self._name)
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD - Mars 2008"
-import sys
import logging
-from daCore import Persistence
-from daCore.BasicObjects import Algorithm
-from daCore import PlatformInfo ; m = PlatformInfo.SystemUsage()
+from daCore import BasicObjects, PlatformInfo
+m = PlatformInfo.SystemUsage()
# ==============================================================================
-class ElementaryAlgorithm(Algorithm):
+class ElementaryAlgorithm(BasicObjects.Algorithm):
def __init__(self):
- Algorithm.__init__(self)
+ BasicObjects.Algorithm.__init__(self)
self._name = "BLUE"
logging.debug("%s Initialisation"%self._name)
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sebastien MASSART, Jean-Philippe ARGAUD - Novembre 2008"
-import sys ; sys.path.insert(0, "../daCore")
import logging
+from daCore import BasicObjects, PlatformInfo
+m = PlatformInfo.SystemUsage()
+
import numpy
-import Persistence
-from BasicObjects import Algorithm
-import PlatformInfo ; m = PlatformInfo.SystemUsage()
# ==============================================================================
-class ElementaryAlgorithm(Algorithm):
+class ElementaryAlgorithm(BasicObjects.Algorithm):
def __init__(self):
- Algorithm.__init__(self)
+ BasicObjects.Algorithm.__init__(self)
self._name = "ENSEMBLEBLUE"
logging.debug("%s Initialisation"%self._name)
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
-import sys ; sys.path.insert(0, "../daCore")
import logging
-import Persistence
-from BasicObjects import Algorithm
-import PlatformInfo ; m = PlatformInfo.SystemUsage()
+from daCore import BasicObjects, PlatformInfo
+m = PlatformInfo.SystemUsage()
# ==============================================================================
-class ElementaryAlgorithm(Algorithm):
+class ElementaryAlgorithm(BasicObjects.Algorithm):
def __init__(self):
- Algorithm.__init__(self)
+ BasicObjects.Algorithm.__init__(self)
self._name = "KALMAN"
logging.debug("%s Initialisation"%self._name)
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Septembre 2008"
-import sys ; sys.path.insert(0, "../daCore")
import logging
-import Persistence
-from BasicObjects import Algorithm
-import PlatformInfo ; m = PlatformInfo.SystemUsage()
+from daCore import BasicObjects, PlatformInfo
+m = PlatformInfo.SystemUsage()
# ==============================================================================
-class ElementaryAlgorithm(Algorithm):
+class ElementaryAlgorithm(BasicObjects.Algorithm):
def __init__(self):
- Algorithm.__init__(self)
+ BasicObjects.Algorithm.__init__(self)
self._name = "LINEARLEASTSQUARES"
def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None):
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
Renvoie les résultats disponibles après l'exécution de la méthode
d'assimilation, ou les diagnostics disponibles. Attention, quand un
- diagnostic porte le même nom qu'un variable stockée, c'est la variable
+ diagnostic porte le même nom qu'une variable stockée, c'est la variable
stockée qui est renvoyée, et le diagnostic est inatteignable.
"""
if key is not None:
allvariables.update( self.__StoredDiagnostics )
return allvariables
+ def get_available_variables(self, key=None):
+ """
+ Renvoie les variables potentiellement utilisables pour l'étude,
+ identifiés par les chaînes de caractères. L'algorithme doit avoir été
+ préalablement choisi sinon la méthode renvoie "None".
+ """
+
+ if type( self.__algorithm ) is type( {} ):
+ return None
+ if key is not None:
+ if self.__algorithm.has_key(key):
+ return self.__algorithm.get( key )
+ else:
+ raise ValueError("The requested key \"%s\" does not exists as a stored variable."%key)
+ else:
+ variables = self.__algorithm.get().keys()
+ variables.sort()
+ return variables
+
def get_available_algorithms(self):
"""
- Renvoie la liste des algorithmes identifiés par les chaînes de
- caractères
+ Renvoie la liste des algorithmes potentiellement utilisables, identifiés
+ par les chaînes de caractères.
"""
files = []
for directory in sys.path:
def get_available_diagnostics(self):
"""
- Renvoie la liste des diagnostics identifiés par les chaînes de
- caractères
+ Renvoie la liste des diagnostics potentiellement utilisables, identifiés
+ par les chaînes de caractères.
"""
files = []
for directory in sys.path:
sys.path = list(set(sys.path)) # Conserve en unique exemplaire chaque chemin
return 1
+ # -----------------------------------------------------------
+ def setDataObserver(self,
+ VariableName = None,
+ HookFunction = None,
+ HookParameters = None,
+ Scheduler = None,
+ ):
+ """
+ Permet d'associer un observer à une variable nommée gérée en interne,
+ activable selon des règles définies dans le Scheduler.
+ """
+ #
+ if type( self.__algorithm ) is dict:
+ raise ValueError("No observer can be build before choosing an algorithm.")
+ #
+ # Vérification du nom de variable et typage
+ # -----------------------------------------
+ if type( VariableName ) is str:
+ VariableNames = [VariableName,]
+ elif type( VariableName ) is list:
+ VariableNames = map( str, VariableName )
+ else:
+ raise ValueError("The observer requires a name or a list of names of variables.")
+ #
+ # Association interne de l'observer à la variable
+ # -----------------------------------------------
+ for n in VariableNames:
+ if not self.__algorithm.has_key( n ):
+ raise ValueError("An observer requires to be set on a variable named %s which does not exist."%n)
+ else:
+ self.__algorithm.StoredVariables[ n ].setDataObserver(
+ Scheduler = Scheduler,
+ HookFunction = HookFunction,
+ HookParameters = HookParameters,
+ )
def prepare_to_pickle(self):
self.__algorithmFile = None
self.__diagnosticFile = None
ADD.analyze()
print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
+ print "Ebauche :", [0, 1, 2]
+ print "Observation :", [0.5, 1.5, 2.5]
+ print "Demi-somme :", list((numpy.array([0, 1, 2])+numpy.array([0.5, 1.5, 2.5]))/2)
+ print " qui doit être identique à :"
print "Analyse résultante :", ADD.get("Analysis").valueserie(0)
print "Innovation :", ADD.get("Innovation").valueserie(0)
print
print "Variables et diagnostics disponibles :", liste
print
+ def obs(var=None,info=None):
+ print " ---> Mise en oeuvre de l'observer"
+ print " var =",var.valueserie(-1)
+ print " info =",info
+ ADD.setDataObserver( 'Analysis', HookFunction=obs, Scheduler = [2, 4], HookParameters = "Second observer")
+ # Attention, il faut décaler le stockage de 1 pour suivre le pas interne
+ # car le pas 0 correspond à l'analyse ci-dessus.
+ for i in range(1,6):
+ print
+ print "Action sur la variable observée, étape :",i
+ ADD.get('Analysis').store( [i, i, i] )
+ print
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
#
self.__steps = []
self.__values = []
+ #
+ self.__dynamic = False
+ #
+ self.__dataobservers = []
def basetype(self, basetype=None):
"""
self.__steps.append(self.__step)
#
self.__values.append(self.__basetype(value))
+ #
+ if self.__dynamic: self.__replot()
+ for hook, parameters, scheduler in self.__dataobservers:
+ if self.__step in scheduler:
+ hook( self, parameters )
def shape(self):
"""
except:
raise TypeError("Base type is incompatible with numpy")
+ def __preplot(self,
+ title = "",
+ xlabel = "",
+ ylabel = "",
+ ltitle = None,
+ geometry = "600x400",
+ persist = False,
+ pause = True,
+ ):
+ import os
+ #
+ # Vérification de la disponibilité du module Gnuplot
+ try:
+ import Gnuplot
+ self.__gnuplot = Gnuplot
+ except:
+ raise ImportError("The Gnuplot module is required to plot the object.")
+ #
+ # Vérification et compléments sur les paramètres d'entrée
+ if persist:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -persist -geometry '+geometry
+ else:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -geometry '+geometry
+ if ltitle is None:
+ ltitle = ""
+ self.__g = self.__gnuplot.Gnuplot() # persist=1
+ self.__g('set terminal '+self.__gnuplot.GnuplotOpts.default_term)
+ self.__g('set style data lines')
+ self.__g('set grid')
+ self.__g('set autoscale')
+ self.__g('set xlabel "'+str(xlabel).encode('ascii','replace')+'"')
+ self.__g('set ylabel "'+str(ylabel).encode('ascii','replace')+'"')
+ self.__title = title
+ self.__ltitle = ltitle
+ self.__pause = pause
+
def plot(self, item=None, step=None,
steps = None,
title = "",
ltitle = None,
geometry = "600x400",
filename = "",
+ dynamic = False,
persist = False,
pause = True,
):
- filename : base de nom de fichier Postscript pour une sauvegarde,
qui est automatiquement complétée par le numéro du
fichier calculé par incrément simple de compteur
+ - dynamic : effectue un affichage des valeurs à chaque stockage
+ (au-delà du second) La méthode "plot" permet de déclarer
+ l'affichage dynamique, et c'est la méthode "__replot"
+ qui est utilisée pour l'effectuer
- persist : booléen indiquant que la fenêtre affichée sera
conservée lors du passage au dessin suivant
Par défaut, persist = False
Par défaut, pause = True
"""
import os
- #
- # Vérification de la disponibilité du module Gnuplot
- try:
- import Gnuplot
- self.__gnuplot = Gnuplot
- except:
- raise ImportError("The Gnuplot module is required to plot the object.")
- #
- # Vérification et compléments sur les paramètres d'entrée
- if persist:
- self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -persist -geometry '+geometry
- else:
- self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -geometry '+geometry
- if ltitle is None:
- ltitle = ""
- self.__g = self.__gnuplot.Gnuplot() # persist=1
- self.__g('set terminal '+self.__gnuplot.GnuplotOpts.default_term)
- self.__g('set style data lines')
- self.__g('set grid')
- self.__g('set autoscale')
- self.__g('set xlabel "'+str(xlabel).encode('ascii','replace')+'"')
- self.__g('set ylabel "'+str(ylabel).encode('ascii','replace')+'"')
+ if not self.__dynamic:
+ self.__preplot(title, xlabel, ylabel, ltitle, geometry, persist, pause )
+ if dynamic:
+ self.__dynamic = True
+ if len(self.__values) == 0: return 0
#
# Tracé du ou des vecteurs demandés
indexes = []
if os.path.isfile(stepfilename):
raise ValueError("Error: a file with this name \"%s\" already exists."%stepfilename)
self.__g.hardcopy(filename=stepfilename, color=1)
- if pause:
+ if self.__pause:
raw_input('Please press return to continue...\n')
+ def __replot(self):
+ """
+ Affichage dans le cas du suivi dynamique de la variable
+ """
+ if self.__dynamic and len(self.__values) < 2: return 0
+ #
+ import os
+ self.__g('set title "'+str(self.__title).encode('ascii','replace'))
+ Steps = range(len(self.__values))
+ self.__g.plot( self.__gnuplot.Data( Steps, self.__values, title=self.__ltitle ) )
+ #
+ if self.__pause:
+ raw_input('Please press return to continue...\n')
+
# ---------------------------------------------------------
def stepmean(self):
"""
- ltitle : titre associé au vecteur tracé
- geometry : taille en pixels de la fenêtre et position du coin haut
gauche, au format X11 : LxH+X+Y (défaut : 600x400)
- - filename : nom de fichier Postscript pour une sauvegarde,
+ - filename : nom de fichier Postscript pour une sauvegarde
- persist : booléen indiquant que la fenêtre affichée sera
conservée lors du passage au dessin suivant
Par défaut, persist = False
if pause:
raw_input('Please press return to continue...\n')
+ # ---------------------------------------------------------
+ def setDataObserver(self,
+ HookFunction = None,
+ HookParameters = None,
+ Scheduler = None,
+ ):
+ """
+ Méthode d'association à la variable d'un triplet définissant un observer
+
+ Le Scheduler attendu est une fréquence, une simple liste d'index ou un
+ xrange des index.
+ """
+ #
+ # Vérification du Scheduler
+ # -------------------------
+ maxiter = int( 1e9 )
+ if type(Scheduler) is int: # Considéré comme une fréquence à partir de 0
+ Schedulers = xrange( 0, maxiter, int(Scheduler) )
+ elif type(Scheduler) is xrange: # Considéré comme un itérateur
+ Schedulers = Scheduler
+ elif type(Scheduler) is list: # Considéré comme des index explicites
+ Schedulers = map( long, Scheduler )
+ else: # Dans tous les autres cas, activé par défaut
+ Schedulers = xrange( 0, maxiter )
+ #
+ # Stockage interne de l'observer dans la variable
+ # -----------------------------------------------
+ self.__dataobservers.append( [HookFunction, HookParameters, Schedulers] )
+
# ==============================================================================
class OneScalar(Persistence):
"""
print "Taille \"len\" du dernier objet stocké",len(OBJET_DE_TEST)
del OBJET_DE_TEST
print
+
+ print "======> Affichage dynamique d'objets"
+ OBJET_DE_TEST = Persistence("My object", unit="", basetype=float)
+ D = OBJET_DE_TEST
+ D.plot(
+ dynamic = True,
+ title = "Valeur suivie",
+ xlabel = "Pas",
+ ylabel = "Valeur",
+ pause = False,
+ )
+ for i in range(1,11):
+ D.store( i*i )
+ print "Taille \"shape\" du dernier objet stocké",OBJET_DE_TEST.shape()
+ print "Taille \"len\" du dernier objet stocké",len(OBJET_DE_TEST)
+ print "Nombre d'objets stockés",OBJET_DE_TEST.stepnumber()
+ del OBJET_DE_TEST
+ print
+
+ print "======> Affectation simple d'observateurs dynamiques"
+ def obs(var=None,info=None):
+ print " ---> Mise en oeuvre de l'observer"
+ print " var =",var.valueserie(-1)
+ print " info =",info
+ OBJET_DE_TEST = Persistence("My object", unit="", basetype=list)
+ D = OBJET_DE_TEST
+ D.setDataObserver( HookFunction = obs )
+ for i in range(5):
+ # print
+ print "Action de 1 observer sur la variable observée, étape :",i
+ D.store( [i, i, i] )
+ del OBJET_DE_TEST
+ print
+
+ print "======> Affectation multiple d'observateurs dynamiques"
+ def obs(var=None,info=None):
+ print " ---> Mise en oeuvre de l'observer"
+ print " var =",var.valueserie(-1)
+ print " info =",info
+ OBJET_DE_TEST = Persistence("My object", unit="", basetype=list)
+ D = OBJET_DE_TEST
+ D.setDataObserver(
+ HookFunction = obs,
+ Scheduler = [2, 4],
+ HookParameters = "Second observer",
+ )
+ D.setDataObserver(
+ HookFunction = obs,
+ Scheduler = xrange(1,3),
+ HookParameters = "Troisième observer",
+ )
+ for i in range(5):
+ # print
+ print "Action de 2 observers sur la variable observée, étape :",i
+ D.store( [i, i, i] )
+ del OBJET_DE_TEST
+ print
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
+#
+# Copyright (C) 2008-2010 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2009 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs
- dependants au sens du test de Student.
- Ce diagnostic utilise le calcul de la p-value pour le test de Student
- pour 2 vecteurs dependants
- En input : la tolerance
- En output : le resultat du diagnostic est une reponse booleenne au test :
- True si les moyennes sont egales au sens du Test de Student
- False dans le cas contraire.
-"""
-__author__ = "Sophie RICCI - Octobre 2008"
-
-import sys ; sys.path.insert(0, "../daCore")
-
-import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from ComputeStudent import DependantVectors
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
- """
- Diagnostic qui effectueIndependantVectorsEqualVariance le test d egalite des moyennes de 2 vecteurs
- dependants au sens du test de Student.
- Ce diagnostic utilise le calcul de la p-value pour le test de Student
- pour 2 vecteurs dependants
- En input : la tolerance
- En output : le resultat du diagnostic est une reponse booleenne au test :
- True si les moyennes sont egales au sens du Test de Student
- False dans le cas contraire.
- """
- def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
- if not self.parameters.has_key("tolerance"):
- raise ValueError("A parameter named \"tolerance\" is required.")
-
- def formula(self, V1, V2):
- """
- Effectue le calcul de la p-value de Student pour deux vecteurs.
- """
- [aire, Q, reponse, message] = DependantVectors(
- vector1 = V1,
- vector2 = V2,
- tolerance = self.parameters["tolerance"] )
- logging.info( message )
- answerStudentTest = False
- if (aire < (100.*self.parameters["tolerance"])) :
- answerStudentTest = False
- else:
- answerStudentTest = True
- return answerStudentTest
-
- def calculate(self, vector1 = None, vector2 = None, step = None):
- """
- Active la formule de calcul
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Student value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- if V1.size != V2.size:
- raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
- value = self.formula( V1, V2 )
- self.store( value = value, step = step)
-
-# ==============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
-
- print " Test d'égalite des moyennes au sens de Student pour deux vecteurs"
- print " dépendants."
- print
- #
- # Initialisation des inputs et appel du diagnostic
- # --------------------------------------------------------------------
- tolerance = 0.05
- D = ElementaryDiagnostic("ComputeMeanStudent_DependVect", parameters = {
- "tolerance":tolerance,
- })
- #
- # Tirage de l'echantillon aleatoire
- # --------------------------------------------------------------------
- x1 = numpy.array(([-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]))
- x2 = numpy.array(([-0.23, 1.36, 0.32, 0.24, -0.66, -0.19, -0.31, 0.56, 1.21, 0.99]))
- #
- # Calcul
- # --------------------------------------------------------------------
- D.calculate(x1, x2)
- #
- if D.valueserie(0) :
- print " L'hypothèse d'égalité des moyennes est valide."
- print
- else :
- raise ValueError("The egality of the means is NOT valid")
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2009 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs
- independants supposes de variances differentes au sens du test de Student.
- En input : la tolerance
- En output : le resultat du diagnostic est une reponse booleenne au test :
- True si les moyennes sont egales au sens du Test de Student
- False dans le cas contraire.
-"""
-__author__ = "Sophie RICCI - Octobre 2008"
-
-import sys ; sys.path.insert(0, "../daCore")
-
-import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from ComputeStudent import IndependantVectorsDifferentVariance
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
- """
- Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs
- independants supposes de variances differentes au sens du test de Student.
- En input : la tolerance
- En output : le resultat du diagnostic est une reponse booleenne au test :
- True si les moyennes sont egales au sens du Test de Student
- False dans le cas contraire.
- """
- def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
- if not self.parameters.has_key("tolerance"):
- raise ValueError("A parameter named \"tolerance\" is required.")
-
- def formula(self, V1, V2):
- """
- Effectue le calcul de la p-value de Student pour deux vecteurs
- independants supposes de variances differentes.
- """
- [aire, Q, reponse, message] = IndependantVectorsDifferentVariance(
- vector1 = V1,
- vector2 = V2,
- tolerance = self.parameters["tolerance"],
- )
- logging.info( message )
- answerStudentTest = False
- if (aire < (100.*self.parameters["tolerance"])) :
- answerStudentTest = False
- else:
- answerStudentTest = True
- return answerStudentTest
-
- def calculate(self, vector1 = None, vector2 = None, step = None):
- """
- Active la formule de calcul
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Student value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- if V1.size != V2.size:
- raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
- value = self.formula( V1, V2 )
- self.store( value = value, step = step)
-
-# ==============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
-
- print " Test d'égalite des moyennes au sens de Student pour deux vecteurs"
- print " indépendants supposés de variances différentes."
- print
- #
- # Initialisation des inputs et appel du diagnostic
- # --------------------------------------------------------------------
- tolerance = 0.05
- D = ElementaryDiagnostic("IndependantVectorsDifferentVariance", parameters = {
- "tolerance":tolerance,
- })
- #
- # Tirage de l'echantillon aleatoire
- # --------------------------------------------------------------------
- x1 = numpy.array(([-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]))
- x2 = numpy.array(([-0.23, 1.36, 0.32, 0.24, -0.66, -0.19, -0.31, 0.56, 1.21, 0.99]))
- #
- # Calcul
- # --------------------------------------------------------------------
- D.calculate(x1, x2)
-#
- if D.valueserie(0) :
- print " L'hypothèse d'égalité des moyennes est valide."
- print
- else :
- raise ValueError("The egality of the means is NOT valid")
-
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2009 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs
- independants supposes de variances egales au sens du test de Student.
- En input : la tolerance
- En output : le resultat du diagnostic est une reponse booleenne au test :
- True si les moyennes sont egales au sens du Test de Student
- False dans le cas contraire.
-"""
-__author__ = "Sophie RICCI - Octobre 2008"
-
-import sys ; sys.path.insert(0, "../daCore")
-
-import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from ComputeStudent import IndependantVectorsEqualVariance
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
- """
- Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs independants supposes de variances egales au sens du test de Student.
- En input : la tolerance
- En output : le resultat du diagnostic est une reponse booleenne au test :
- True si les moyennes sont egales au sens du Test de Student
- False dans le cas contraire.
- """
- def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
- if not self.parameters.has_key("tolerance"):
- raise ValueError("A parameter named \"tolerance\" is required.")
-
- def formula(self, V1, V2):
- """
- Effectue le calcul de la p-value de Student pour deux vecteurs
- independants supposes de variances egales.
- """
- [aire, Q, reponse, message] = IndependantVectorsEqualVariance(
- vector1 = V1,
- vector2 = V2,
- tolerance = self.parameters["tolerance"],
- )
- logging.info( message )
- answerStudentTest = False
- if (aire < (100.*self.parameters["tolerance"])) :
- answerStudentTest = False
- else:
- answerStudentTest = True
- return answerStudentTest
-
- def calculate(self, vector1 = None, vector2 = None, step = None):
- """
- Active la formule de calcul
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Student value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- if V1.size != V2.size:
- raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
- value = self.formula( V1, V2 )
- self.store( value = value, step = step)
-
-# ==============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
-
- print " Test d'égalite des moyennes au sens de Student pour deux vecteurs"
- print " indépendants supposés de variances égales"
- print
- #
- # Initialisation des inputs et appel du diagnostic
- # --------------------------------------------------------------------
- tolerance = 0.05
- D = ElementaryDiagnostic("ComputeMeanStudent_IndepVect_EgalVar", parameters = {
- "tolerance":tolerance,
- })
- #
- # Tirage de l'echantillon aleatoire
- # --------------------------------------------------------------------
- x1 = numpy.array(([-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]))
- x2 = numpy.array(([-0.23, 1.36, 0.32, 0.24, -0.66, -0.19, -0.31, 0.56, 1.21, 0.99]))
- #
- # Calcul
- # --------------------------------------------------------------------
- D.calculate(x1, x2)
- #
- if D.valueserie(0) :
- print " L'hypothèse d'égalité des moyennes est valide."
- print
- else :
- raise ValueError("The egality of the means is NOT valid")
-
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2009 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic qui compare les variances de 2 vecteurs au sens de Fisher à
- l'aide du calcul de la p-value pour le test de Fisher.
- - entrée : la tolérance (tolerance) sous forme de paramètres dans le
- dictionnaire Par, et les deux vecteurs d'échantillons.
- - sortie : le résultat du diagnostic est une réponse booléenne au test :
- True si l'égalite des variances est valide au sens du test de Fisher,
- False dans le cas contraire
-"""
-__author__ = "Sophie RICCI - Juillet 2008"
-
-import sys ; sys.path.insert(0, "../daCore")
-
-import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from ComputeFisher import ComputeFisher
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
- """
- Diagnostic qui compare les variances de 2 vecteurs au sens de Fisher à
- l'aide du calcul de la p-value pour le test de Fisher.
- - entrée : la tolérance (tolerance) sous forme de paramètres dans le
- dictionnaire Par, et les deux vecteurs d'échantillons.
- - sortie : le résultat du diagnostic est une réponse booléenne au test :
- True si l'égalite des variances est valide au sens du test de Fisher,
- False dans le cas contraire
- """
- def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
- if not self.parameters.has_key("tolerance"):
- raise ValueError("A parameter named \"tolerance\" is required.")
-
- def formula(self, V1, V2):
- """
- Effectue le test de Fisher avec la p-value pour 2 vecteurs
- """
- [aire, f, reponse, message] = ComputeFisher(
- vector1 = V1,
- vector2 = V2,
- tolerance = self.parameters["tolerance"],
- )
- answerKhisquareTest = False
- if (aire < (100.*self.parameters["tolerance"])) :
- answerKhisquareTest = False
- else:
- answerKhisquareTest = True
- logging.info( message )
- #
- return answerKhisquareTest
-
- def calculate(self, vector1 = None, vector2 = None, step = None):
- """
- Active la formule de calcul
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Fisher p-value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- if V1.size != V2.size:
- raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
- #
- value = self.formula( V1, V2 )
- #
- self.store( value = value, step = step)
-
-# ==============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
-
- print " Test d'égalite des variances pour deux vecteurs de taille 10"
- print
- #
- # Initialisation des inputs et appel du diagnostic
- # --------------------------------------------------------------------
- tolerance = 0.05
- D = ElementaryDiagnostic("CompareVarianceFisher", parameters = {
- "tolerance":tolerance,
- })
- #
- # Tirage de l'echantillon aleatoire
- # --------------------------------------------------------------------
- x1 = numpy.array(([-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]))
- x2 = numpy.array(([-0.23, 1.36, 0.32, 0.24, -0.66, -0.19, -0.31, 0.56, 1.21, 0.99]))
- #
- # Calcul
- # --------------------------------------------------------------------
- D.calculate(x1, x2)
- #
- if D.valueserie(0) :
- print " L'hypothèse d'égalité des deux variances est correcte."
- print
- else :
- raise ValueError("L'hypothèse d'égalité des deux variances est fausse.")
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sophie RICCI - Aout 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from AssimilationStudy import AssimilationStudy
+from daCore import BasicObjects, Persistence
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = float )
def _formula(self, V):
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sophie RICCI - Octobre 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from AssimilationStudy import AssimilationStudy
+from daCore import BasicObjects, Persistence
import logging
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name)
+ BasicObjects.Diagnostic.__init__(self, name)
Persistence.OneScalar.__init__( self, name, unit, basetype = float)
def _formula(self, X, HX, Xb, Y, R, B):
"""
Calcul de la fonction cout
"""
- Jb = 1./2. * (X - Xb).T * B.I * (X - Xb)
+# Jb = 1./2. * (X - Xb).T * B.I * (X - Xb)
+ Jb = 1./2. * numpy.dot((X - Xb) ,numpy.asarray(numpy.dot(B.I,(X - Xb)).A1))
logging.info( "Partial cost function : Jb = %s"%Jb )
#
- Jo = 1./2. * (Y - HX).T * R.I * (Y - HX)
+# Jo = 1./2. * (Y - HX).T * R.I * (Y - HX)
+ Jo = 1./2. * numpy.dot((Y - HX) ,numpy.asarray(numpy.dot(R.I,(Y - HX)).A1))
logging.info( "Partial cost function : Jo = %s"%Jo )
#
J = Jb + Jo
xb = numpy.array([2., 2.])
yo = numpy.array([5., 6.])
H = numpy.matrix(numpy.identity(2))
- Hx = H*x
+# Hx = H*x
+ Hx = numpy.dot(H,x)
Hx = Hx.T
B = numpy.matrix(numpy.identity(2))
R = numpy.matrix(numpy.identity(2))
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2009 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Calcul de la fonction coût avec Hlin
- HX = Hxb + Hlin dx
-"""
-__author__ = "Sophie RICCI - Octobre 2008"
-
-import sys ; sys.path.insert(0, "../daCore")
-
-import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from AssimilationStudy import AssimilationStudy
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
- def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name)
- Persistence.OneScalar.__init__( self, name, unit, basetype = float)
- self.__name = str( name )
-
- def _formula(self, X = None, dX = None, Hlin = None, Xb=None, HXb = None, Y=None, R=None, B=None):
-
- """
- Calcul de la fonction cout
- """
- HX = HXb + Hlin.T * dX
- if hasattr(HX, 'A1') :
- HX = HX.A1
- #
- Jb = 1./2. * (X - Xb).T * B.I * (X - Xb)
- logging.info( "Partial cost function : Jb = %s"%Jb )
- #
- Jo = 1./2. * (Y - HX).T * R.I * (Y - HX)
- logging.info( "Partial cost function : Jo = %s"%Jo )
- #
- J = Jb + Jo
- logging.info( "Total cost function : J = Jo + Jb = %s"%J )
- return J
-
- def calculate(self, x = None, dx = None, Hlin = None, xb = None, Hxb = None, yo = None, R = None, B = None , step = None):
- """
- Teste les arguments, active la formule de calcul et stocke le résultat
- """
- if (x is None) or (xb is None) or (yo is None) or (dx is None):
- raise ValueError("Vectors x, dx, xb and yo must be given to compute J")
- dX = dx
- if hasattr(numpy.matrix(x), 'A1') :
- X = numpy.matrix(x).A1
- if hasattr(numpy.matrix(xb), 'A1') :
- Xb = numpy.matrix(xb).A1
- if hasattr(numpy.matrix(yo), 'A1') :
- Y = numpy.matrix(yo).A1
- B = numpy.matrix(B)
- R = numpy.matrix(R)
- if (Hlin is None ) :
- raise ValueError("HlinT vector must be given")
- if (Hxb is None ) :
- raise ValueError("The given vector must be given")
- HXb = Hxb
- if (B is None ) or (R is None ):
- raise ValueError("The matrices B and R must be given")
- #
- value = self._formula(X, dX, Hlin, Xb, HXb, Y, R, B)
- #
- self.store( value = value, step = step )
-
-#===============================================================================
-if __name__ == "__main__":
- print "\nAUTOTEST\n"
- #
- D = ElementaryDiagnostic("Ma fonction cout")
- #
- # Vecteur de type array
- # ---------------------
- x = numpy.array([1., 2.])
- dx = numpy.array([0.1, 0.2])
- xb = numpy.array([2., 2.])
- yo = numpy.array([5., 6.])
- Hlin = numpy.matrix(numpy.identity(2))
- Hxb = Hlin *xb
- Hxb = Hxb.T
- Hxb = Hxb.A1
- B = numpy.matrix(numpy.identity(2))
- R = numpy.matrix(numpy.identity(2))
- #
- D.calculate( x = x, dx = dx, Hlin = Hlin, xb = xb, Hxb = Hxb, yo = yo, R = R, B = B)
- print "Le vecteur x choisi est...:", x
- print "L ebauche xb choisie est...:", xb
- print "Le vecteur d observation est...:", yo
- print "B = ", B
- print "R = ", R
- print "La fonction cout J vaut ...: %.2e"%D.valueserie(0)
- #
- if (abs(D.valueserie(0) - 11.925) > 1.e-6) :
- raise ValueError("The computation of the cost function is NOT correct")
- else :
- print "The computation of the cost function is OK"
- print
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from AssimilationStudy import AssimilationStudy
+from daCore import BasicObjects, Persistence
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = float)
def _formula(self, V):
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sophie RICCI - Juillet 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-from numpy import random
-import Persistence
-from BasicObjects import Diagnostic
+from daCore import BasicObjects, Persistence
from ComputeKhi2 import ComputeKhi2_Gauss
import logging
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
"""
"""
def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
for key in ["tolerance", "dxclasse", "nbclasses"]:
if not self.parameters.has_key(key):
nbclasses = self.parameters["nbclasses"],
SuppressEmptyClasses = True)
-
logging.info( message )
logging.info( "(si <%.2f %s on refuse effectivement l'adéquation)"%(100.*self.parameters["tolerance"],"%") )
logging.info("vecteur des classes=%s"%numpy.size(vectclasse) )
# Tirage de l'echantillon aleatoire
# ---------------------------------
numpy.random.seed(2490)
- x = random.normal(50.,1.5,1000)
+ x = numpy.random.normal(50.,1.5,1000)
#
# Calcul
# ------
# Tirage de l'echantillon aleatoire
# ---------------------------------
numpy.random.seed(2490)
- x = random.normal(50.,1.5,1000)
+ x = numpy.random.normal(50.,1.5,1000)
#
# Calcul
# ------
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sophie RICCI - Septembre 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from RMS import ElementaryDiagnostic as RMS
-from AssimilationStudy import AssimilationStudy
+from daCore import AssimilationStudy, BasicObjects, Persistence
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
if not self.parameters.has_key("tolerance"):
raise ValueError("A parameter named \"tolerance\" is required.")
#
# Calcul de l'ecart entre Hx1 et Hx et entre Hx2 et Hx
# ----------------------------------------------------
- ADD = AssimilationStudy()
+ ADD = AssimilationStudy.AssimilationStudy()
ADD.setDiagnostic("RMS",
name = "Calcul de la RMS entre Hx1 et Hx et entre Hx2 et Hx")
RMS = ADD.get("Calcul de la RMS entre Hx1 et Hx et entre Hx2 et Hx")
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sophie RICCI - Juillet 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-from numpy import random
-
-import Persistence
-from BasicObjects import Diagnostic
+from daCore import BasicObjects, Persistence
from ComputeKhi2 import ComputeKhi2_Homogen
import logging
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = bool )
for key in ["tolerance", "dxclasse", "nbclasses"]:
if not self.parameters.has_key(key):
# Tirage de l'echantillon aleatoire
# --------------------------------------------------------------------
numpy.random.seed(4000)
- x1 = random.normal(50.,1.5,10000)
+ x1 = numpy.random.normal(50.,1.5,10000)
numpy.random.seed(2490)
- x2 = random.normal(50.,1.5,10000)
+ x2 = numpy.random.normal(50.,1.5,10000)
#
# Calcul
# --------------------------------------------------------------------
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD - Juillet 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import os.path
import numpy
-from BasicObjects import Diagnostic
+from daCore import BasicObjects
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic):
+class ElementaryDiagnostic(BasicObjects.Diagnostic):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
try:
import Gnuplot
self.__gnuplot = Gnuplot
return 1
def calculate(self, vector = None, steps = None,
- title = "", xlabel = "", ylabel = "", ltitle = None,
+ title = "", xlabel = "", ylabel = "", ltitle = "",
geometry = "600x400",
filename = "",
persist = False,
"""
if vector is None:
raise ValueError("One vector must be given to plot it.")
- if ltitle is None:
- ltitle = ""
Vector = numpy.array(vector)
if Vector.size < 1:
raise ValueError("The given vector must not be empty")
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import os.path
import numpy
-from BasicObjects import Diagnostic
+from daCore import BasicObjects
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic):
+class ElementaryDiagnostic(BasicObjects.Diagnostic):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
try:
import Gnuplot
self.__gnuplot = Gnuplot
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD - Juillet 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
-import math
-import numpy
-import Persistence
-from BasicObjects import Diagnostic
+import math, numpy
+from daCore import BasicObjects, Persistence
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = float)
def _formula(self, V1, V2):
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sophie RICCI - Aout 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from AssimilationStudy import AssimilationStudy
+from daCore import BasicObjects, Persistence
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
def _formula(self, V1, V2):
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from AssimilationStudy import AssimilationStudy
+from daCore import BasicObjects, Persistence
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = bool )
def _formula(self, V1, V2):
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2009 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic sur l'arrêt (ou le ralentissement) de la réduction de la variance
- au fil des pas (ou itérations) de l'analyse.
- Ce diagnostic s'applique typiquement au vecteur de différence entre la
- variance de OMB et la variance de OMA au fil du temps ou des itérations:
- V[i] = vecteur des VAR(OMB)[i] - VAR(OMA)[i] au temps ou itération i.
-"""
-__author__ = "Sophie Ricci - Septembre 2008"
-
-import sys ; sys.path.insert(0, "../daCore")
-
-import numpy
-import Persistence
-from BasicObjects import Diagnostic
-from AssimilationStudy import AssimilationStudy
-
-# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
- def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = int )
-
- def _formula(self, V, CutOffSlope, MultiSlope0):
- """
- Recherche du pas de temps ou iteration pour laquelle la reduction
- de la variance est
- - inferieure a la valeur seuil CutOffSlope
- (si une valeure est donnee a CutOffSlope)
- - inferieure a MultiSlope0 * la pente a la premiere iteration
- (si une valeure est donnee a MultiSlope0)
- V[i] = vecteur des VAR(OMB)[i] - VAR(OMA)[i] au temps ou iteration i.
- """
- N = V.size
- pente = numpy.matrix(numpy.zeros((N,))).T
- iterstopreduction = 0.
- for i in range (1, N) :
- pente[i] = V[i]- V[i-1]
- if pente[i] > 0.0 :
- raise ValueError("The analysis is INCREASING the variance a l iteration ", i)
- if CutOffSlope is not None:
- if numpy.abs(pente[i]) < CutOffSlope :
- iterstopreduction = i
- break
- if MultiSlope0 is not None:
- if numpy.abs(pente[i]) < MultiSlope0 * numpy.abs(pente[1]) :
- iterstopreduction = i
- break
- #
- return iterstopreduction
-
- def calculate(self, vector = None, CutOffSlope = None, MultiSlope0 = None, step = None) :
- """
- Teste les arguments, active la formule de calcul et stocke le resultat
- Arguments :
- - vector : vecteur des VAR(OMB) - VAR(OMA) au fil des iterations
- - CutOffSlope : valeur minimale de la pente
- - MultiSlope0 : Facteur multiplicatif de la pente initiale pour comparaison
- """
- if (vector is None) :
- raise ValueError("One vector must be given to test the convergence of the variance after analysis")
- V = numpy.array(vector)
- if V.size < 1 :
- raise ValueError("The given vector must not be empty")
- if (MultiSlope0 is None) and (CutOffSlope is None) :
- raise ValueError("You must set the value of ONE of the CutOffSlope of MultiSlope0 key word")
- #
- value = self._formula( V, CutOffSlope, MultiSlope0 )
- #
- self.store( value = value, step = step )
-
-#===============================================================================
-if __name__ == "__main__":
- print "\n AUTODIAGNOSTIC \n"
-
- # Instanciation de l'objet diagnostic
- # ------------------------------------------------
- D = ElementaryDiagnostic("Mon StopReductionVariance")
-
- # Vecteur de reduction VAR(OMB)-VAR(OMA)
- # ------------------------------------------------
- x = numpy.array(([0.60898111, 0.30449056, 0.15224528, 0.07612264, 0.03806132, 0.01903066, 0.00951533, 0.00475766, 0.00237883, 0.00118942]))
- print " Le vecteur choisi est :", x
- print " Sur ce vecteur, la reduction a l iteration N = 7 est inferieure a 0.005"
- print " Sur ce vecteur, la reduction a l iteration N = 8 est inferieure a 0.01 * la reduction a l iteration 1"
-
- # Comparaison a la valeur seuil de la reduction
- # ------------------------------------------------
- D.calculate( vector = x, CutOffSlope = 0.005, MultiSlope0 = None)
- if (D.valueserie(0) - 7.) < 1.e-15 :
- print " Test : La comparaison a la valeur seuil de la reduction est juste"
- else :
- print " Test : La comparaison a la valeur seuil de la reduction est fausse"
-
- # Comparaison a alpha* la reduction a la premiere iteration
- # ------------------------------------------------
- D.calculate( vector = x, CutOffSlope = None, MultiSlope0 = 0.01)
- if (D.valueserie(1) - 8.) < 1.e-15 :
- print " Test : La comparaison a la reduction a la premiere iteration est juste"
- else :
- print " Test : La comparaison a la reduction a la premiere iteration est fausse"
- print
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Septembre 2008"
-import sys ; sys.path.insert(0, "../daCore")
-
import numpy
-import Persistence
-from BasicObjects import Diagnostic
from scipy.linalg import eig
+from daCore import BasicObjects, Persistence
import logging
# ==============================================================================
-class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- Diagnostic.__init__(self, name, parameters)
+ BasicObjects.Diagnostic.__init__(self, name, parameters)
Persistence.OneScalar.__init__( self, name, unit, basetype = bool )
def _formula(self, xb, B, yo, R):
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
+++ /dev/null
-<?xml version='1.0' encoding='iso-8859-1' ?>
-<proc>
- <property name="DefaultStudyID" value="1"/>
- <type name="bool" kind="bool"/>
- <sequence name="boolvec" content="bool"/>
- <type name="double" kind="double"/>
- <sequence name="dblevec" content="double"/>
- <objref name="file" id="file"/>
- <type name="int" kind="int"/>
- <sequence name="intvec" content="int"/>
- <objref name="pyobj" id="python:obj:1.0"/>
- <sequence name="seqdblevec" content="dblevec"/>
- <type name="string" kind="string"/>
- <sequence name="stringvec" content="string"/>
- <container name="DefaultContainer">
- <property name="container_name" value="FactoryServer"/>
- <property name="hostname" value="localhost"/>
- </container>
- <inline name="Building_B">
- <script><code><![CDATA[# Construction de la matrice de covariances d'erreurs d'ebauche
-import numpy
-
-dimension = len( parametres )
-
-B = numpy.matrix(numpy.core.identity(dimension))
-B = B.A1
-]]></code></script>
- <inport name="parametres" type="pyobj"/>
- <outport name="B" type="dblevec"/>
- </inline>
- <inline name="Building_R">
- <script><code><![CDATA[# Construction de la matrice de covariances d'erreurs de mesure
-import numpy
-
-dimension = len( experiences ) * len( experiences[0] )
-
-R = numpy.matrix(numpy.core.identity(dimension))
-R = R.A1
-]]></code></script>
- <inport name="experiences" type="pyobj"/>
- <outport name="R" type="dblevec"/>
- </inline>
- <inline name="Building_Xb">
- <script><code><![CDATA[# Script pour extraire l'ebauche comme une liste
-# et une liste de bornes [min,max]
-
-dx = 1.e-2 # Increment en pourcent de Xb
-
-Xb = []
-dXb = []
-Bornes = []
-
-for parametre in parametres:
- Xb.append( parametre[1] )
- dXb.append( dx * parametre[1] )
- Bornes.append( parametre[2:4] )
-]]></code></script>
- <inport name="parametres" type="pyobj"/>
- <outport name="Xb" type="dblevec"/>
- <outport name="Bornes" type="seqdblevec"/>
- <outport name="dXb" type="dblevec"/>
- </inline>
- <inline name="Building_Yo">
- <script><code><![CDATA[# Script pour extraire l'es mesures/observations et leur nom
-# comme une liste
-
-if len(calcul) != len(experiences):
- raise ValueError("Les nombres de variables calculees et observees doient etre les memes")
-
-if len(experiences) >=1:
- nb_observations_par_experience = len( experiences[0] )
-
-Yo = []
-Nom_Yo = []
-
-for i in range( len(experiences) ):
- Nom_Yo.append( calcul[i][2] )
- if len( experiences[i] ) != nb_observations_par_experience:
- raise ValueError("Le nombre de mesures par serie experimentale doit etre constant egal a %i mais la serie %i est longue de %i"%(nb_observations_par_experience,i,len(experiences[i])))
- for observation in experiences[i]:
- Yo.append( observation[1] )
-]]></code></script>
- <inport name="calcul" type="pyobj"/>
- <inport name="experiences" type="pyobj"/>
- <outport name="Yo" type="dblevec"/>
- <outport name="Nom_Yo" type="stringvec"/>
- </inline>
- <inline name="Entrees du calcul AD">
- <script><code><![CDATA[import numpy
-
-print
-print "### ============================================================="
-print "### Verification des arguments fabriques pour passer d'Aster a AD"
-print
-print " Tailles et types :"
-print " de Xb",len(Xb),type(Xb)
-print " de Yo",len(Yo),type(Yo)
-print " de B ",len(B),type(B)
-print " de R ",len(R),type(R)
-print " de H ",len(H),type(H)
-print
-print " Verification de remise a l'echelle :"
-dimensionXb = len( Xb )
-dimensionYo = len( Yo )
-print " de B"
-numpyB = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
-print " de R"
-numpyR = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
-print " de H"
-numpyH = numpy.matrix( H, numpy.float ).reshape((dimensionYo,dimensionXb))
-print "### ============================================================="
-print
-]]></code></script>
- <inport name="Xb" type="dblevec"/>
- <inport name="Yo" type="dblevec"/>
- <inport name="B" type="dblevec"/>
- <inport name="R" type="dblevec"/>
- <inport name="H" type="dblevec"/>
- <outport name="Xb" type="dblevec"/>
- <outport name="Yo" type="dblevec"/>
- <outport name="B" type="dblevec"/>
- <outport name="R" type="dblevec"/>
- <outport name="H" type="dblevec"/>
- </inline>
- <inline name="Sorties du calcul AD">
- <script><code><![CDATA[import numpy
-print "Diagnostics de sortie de test :"
-print " Remise en numpy des resultats"
-Xa = numpy.array(xa)
-Xb = numpy.array(xb)
-dimensionXb = len( Xb )
-dimensionYo = len( Yo )
-print "Verification de remise a l'echelle de B, R et H"
-B = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
-R = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
-H = numpy.matrix( H, numpy.float ).reshape((dimensionYo,dimensionXb))
-I = numpy.matrix(numpy.core.identity(dimensionYo))
-deltaB = Xa - Xb
-deltaA = (Yo - numpy.dot(H,Xa)).A1
-deltaI = (Yo - numpy.dot(H,Xb)).A1
-print "Calcul de LS"
-print "deltaA",deltaA
-LS = float(numpy.dot(deltaA, deltaA))
-print "Calcul de LSI"
-LSI = float(numpy.dot(deltaI, deltaI))
-print "Calcul de J"
-J = float( numpy.dot(deltaB, numpy.dot(B,deltaB).A1) + numpy.dot(deltaA, numpy.dot(R,deltaA).A1) )
-print "Calcul de JI"
-JI = float( numpy.dot(deltaI, numpy.dot(R,deltaI).A1) )
-
-print
-print "Sortie du test :"
-print " Xb =",Xb
-print " Xa =",Xa
-print "Difference Xa-Xb :"
-print " Xa-Xb =",deltaB
-print " max(Xa-Xb) =",max(deltaB)
-print " min(Xa-Xb) =",min(deltaB)
-print "Fonctionnelles d'ecarts :"
-print " Initialement : J =",JI
-print " LS =",LSI
-print " Analyse : J =",J,"(Baise de %i%s)"%(100.*(JI-J)/J,"%")
-print " LS =",LS,"(Baise de %i%s)"%(100.*(LSI-LS)/LSI,"%")
-print "Autres informations :"
-print " d =",Innovation
-print " A =",A
-print
-]]></code></script>
- <inport name="xa" type="dblevec"/>
- <inport name="A" type="dblevec"/>
- <inport name="Innovation" type="dblevec"/>
- <inport name="xb" type="dblevec"/>
- <inport name="Yo" type="dblevec"/>
- <inport name="B" type="dblevec"/>
- <inport name="R" type="dblevec"/>
- <inport name="H" type="dblevec"/>
- </inline>
- <parameter>
- <tonode>Building_B</tonode><toport>parametres</toport>
- <value><objref>(lp1
-.</objref></value>
- </parameter>
- <parameter>
- <tonode>Building_Xb</tonode><toport>parametres</toport>
- <value><objref>(lp1
-.</objref></value>
- </parameter>
- <parameter>
- <tonode>Building_Yo</tonode><toport>calcul</toport>
- <value><objref>(lp1
-.</objref></value>
- </parameter>
- <parameter>
- <tonode>Building_Yo</tonode><toport>experiences</toport>
- <value><objref>(lp1
-.</objref></value>
- </parameter>
- <parameter>
- <tonode>Sorties du calcul AD</tonode><toport>xa</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Sorties du calcul AD</tonode><toport>A</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Sorties du calcul AD</tonode><toport>Innovation</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Sorties du calcul AD</tonode><toport>xb</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Sorties du calcul AD</tonode><toport>Yo</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Sorties du calcul AD</tonode><toport>B</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Sorties du calcul AD</tonode><toport>R</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Sorties du calcul AD</tonode><toport>H</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Building_R</tonode><toport>experiences</toport>
- <value><objref>(lp1
-.</objref></value>
- </parameter>
- <parameter>
- <tonode>Entrees du calcul AD</tonode><toport>Xb</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Entrees du calcul AD</tonode><toport>Yo</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Entrees du calcul AD</tonode><toport>B</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Entrees du calcul AD</tonode><toport>R</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>Entrees du calcul AD</tonode><toport>H</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <presentation name="Building_B" x="9" y="81.5" width="227.5" height="110"/>
- <presentation name="Building_Xb" x="9" y="233" width="227.5" height="168"/>
- <presentation name="Building_Yo" x="249" y="231" width="227.5" height="139"/>
- <presentation name="Sorties du calcul AD" x="246" y="432" width="223" height="313"/>
- <presentation name="Building_R" x="250" y="84.5" width="227.5" height="110"/>
- <presentation name="Entrees du calcul AD" x="9" y="429" width="227.5" height="226"/>
- <presentation name="__ROOT__" x="0" y="0" width="487" height="754"/>
-</proc>
+++ /dev/null
-<?xml version='1.0' encoding='iso-8859-1' ?>
-<proc>
- <property name="DefaultStudyID" value="1"/>
- <type name="bool" kind="bool"/>
- <sequence name="boolvec" content="bool"/>
- <type name="double" kind="double"/>
- <sequence name="dblevec" content="double"/>
- <objref name="file" id="file"/>
- <type name="int" kind="int"/>
- <sequence name="intvec" content="int"/>
- <objref name="pyobj" id="python:obj:1.0"/>
- <sequence name="seqdblevec" content="dblevec"/>
- <sequence name="seqint" content="int"/>
- <sequence name="seqintvec" content="intvec"/>
- <type name="string" kind="string"/>
- <sequence name="stringvec" content="string"/>
- <container name="DefaultContainer">
- <property name="container_name" value="FactoryServer"/>
- <property name="hostname" value="localhost"/>
- </container>
- <bloc name="H_linearization">
- <foreach name="Finite_differences_derivation" type="int">
- <bloc name="Elementary_calculation">
- <inline name="Perturbated_point_X">
- <script><code><![CDATA[print iter
-print seq_X[iter]
-X = seq_X[iter]
-]]></code></script>
- <inport name="iter" type="int"/>
- <inport name="seq_X" type="seqdblevec"/>
- <outport name="X" type="dblevec"/>
- <outport name="iter" type="int"/>
- </inline>
- <inline name="ASTER">
- <script><code><![CDATA[print 'Debut ASTER_recal'
-
-import os
-execfile( os.path.join(SOURCES_ROOT, 'N_Code_Aster.py') )
-
-print 'RESU_CALC (%s): %s ' % (iter, RESU_CALC)
-print 'DIAG (%s): %s' % (iter, DIAG)
-]]></code></script>
- <inport name="X" type="dblevec"/>
- <inport name="iter" type="int"/>
- <inport name="ASTER_ROOT" type="string"/>
- <inport name="rcdir" type="string"/>
- <inport name="debug" type="bool"/>
- <inport name="DISPLAY" type="string"/>
- <inport name="SOURCES_ROOT" type="string"/>
- <inport name="export" type="string"/>
- <inport name="parametres" type="pyobj"/>
- <inport name="calcul" type="pyobj"/>
- <inport name="experience" type="pyobj"/>
- <inport name="fileparameters" type="string"/>
- <outport name="FX" type="dblevec"/>
- <outport name="FY" type="dblevec"/>
- <outport name="DIMS" type="intvec"/>
- <outport name="DIAG" type="string"/>
- <outport name="iter" type="int"/>
- </inline>
- <control> <fromnode>Perturbated_point_X</fromnode> <tonode>ASTER</tonode> </control>
- <datalink control="false">
- <fromnode>Perturbated_point_X</fromnode> <fromport>X</fromport>
- <tonode>ASTER</tonode> <toport>X</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Perturbated_point_X</fromnode> <fromport>iter</fromport>
- <tonode>ASTER</tonode> <toport>iter</toport>
- </datalink>
- </bloc>
- </foreach>
- <inline name="Gradient">
- <script><code><![CDATA[
-print "seq_FX=", seq_FX
-print "seq_FY=", seq_FY
-print "seq_DIMS=", seq_DIMS
-print "lst_DIAG=", lst_DIAG
-print "lst_iter=", lst_iter
-print "dX=", dX
-
-import os
-execfile( os.path.join(SOURCES_ROOT, 'N_Gradient.py') )
-
-print "fonctionnelle=", fonctionnelle
-print "gradient=", gradient
-]]></code></script>
- <inport name="seq_FX" type="seqdblevec"/>
- <inport name="seq_FY" type="seqdblevec"/>
- <inport name="seq_DIMS" type="seqintvec"/>
- <inport name="lst_DIAG" type="stringvec"/>
- <inport name="lst_iter" type="intvec"/>
- <inport name="dX" type="dblevec"/>
- <inport name="SOURCES_ROOT" type="string"/>
- <outport name="fonctionnelle" type="dblevec"/>
- <outport name="gradient" type="seqdblevec"/>
- </inline>
- <inline name="Input">
- <function name="inputctrl">
- <code><![CDATA[# debug : nb permet de limiter le nb de directions pour les differences finies
-nb = None # on calcule pour tous les parametres
-#nb = 1
-
-print
-print "XXXXXXXXXXXXX======== Demarrage de la procedure"
-print
-import copy
-
-def inputctrl( X, dX ):
-
- if len( X ) != len( dX ):
- raise ValueError("X and dX has to present the same lenght")
-
- nbparam = len(X)
- if nb: nbparam = min(nb, nbparam)
- seq_X = [ X ]
- for i in range(nbparam):
- Xplus = copy.copy(X)
- Xplus[i] = X[i] + dX[i]
- seq_X.append( Xplus )
-
- nb_core = 2
- itervect = range( len(seq_X) )
- nbBranches = min(nb_core, len( itervect ))
-
- print " Taille du vecteur :", nbparam
- print "seq_X:", seq_X
- print "nbBranches:", nbBranches
- print "itervect:", itervect
-
- # debug
- #nbBranches = 1
- #if nbBranches==1:
- # itervect = [0]
- # seq_X = [ X ]
-
- return nbBranches, itervect, seq_X, dX
-]]></code>
- </function>
- <inport name="X" type="dblevec"/>
- <inport name="dX" type="dblevec"/>
- <outport name="nbBranches" type="int"/>
- <outport name="itervect" type="intvec"/>
- <outport name="seq_X" type="seqdblevec"/>
- <outport name="dX" type="dblevec"/>
- </inline>
- <inline name="Temporary_Parameters">
- <script><code><![CDATA[# Boitier pour une interface simple de variables temporaires]]></code></script>
- <inport name="ASTER_ROOT" type="string"/>
- <inport name="rcdir" type="string"/>
- <inport name="debug" type="bool"/>
- <inport name="DISPLAY" type="string"/>
- <inport name="SOURCES_ROOT" type="string"/>
- <inport name="export" type="string"/>
- <inport name="parametres" type="pyobj"/>
- <inport name="calcul" type="pyobj"/>
- <inport name="experience" type="pyobj"/>
- <inport name="fileparameters" type="string"/>
- <outport name="ASTER_ROOT" type="string"/>
- <outport name="rcdir" type="string"/>
- <outport name="debug" type="bool"/>
- <outport name="DISPLAY" type="string"/>
- <outport name="SOURCES_ROOT" type="string"/>
- <outport name="export" type="string"/>
- <outport name="parametres" type="pyobj"/>
- <outport name="calcul" type="pyobj"/>
- <outport name="experience" type="pyobj"/>
- <outport name="fileparameters" type="string"/>
- </inline>
- <control> <fromnode>Finite_differences_derivation</fromnode> <tonode>Gradient</tonode> </control>
- <control> <fromnode>Input</fromnode> <tonode>Finite_differences_derivation</tonode> </control>
- <control> <fromnode>Input</fromnode> <tonode>Gradient</tonode> </control>
- <control> <fromnode>Temporary_Parameters</fromnode> <tonode>Finite_differences_derivation</tonode> </control>
- <control> <fromnode>Temporary_Parameters</fromnode> <tonode>Gradient</tonode> </control>
- <datalink control="false">
- <fromnode>Finite_differences_derivation</fromnode> <fromport>SmplPrt</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.Perturbated_point_X</tonode> <toport>iter</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Input</fromnode> <fromport>nbBranches</fromport>
- <tonode>Finite_differences_derivation</tonode> <toport>nbBranches</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Input</fromnode> <fromport>itervect</fromport>
- <tonode>Finite_differences_derivation</tonode> <toport>SmplsCollection</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Input</fromnode> <fromport>seq_X</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.Perturbated_point_X</tonode> <toport>seq_X</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Input</fromnode> <fromport>dX</fromport>
- <tonode>Gradient</tonode> <toport>dX</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>ASTER_ROOT</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>ASTER_ROOT</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>rcdir</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>rcdir</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>debug</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>debug</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>DISPLAY</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>DISPLAY</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>SOURCES_ROOT</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>SOURCES_ROOT</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>SOURCES_ROOT</fromport>
- <tonode>Gradient</tonode> <toport>SOURCES_ROOT</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>export</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>export</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>parametres</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>parametres</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>calcul</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>calcul</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>experience</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>experience</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Temporary_Parameters</fromnode> <fromport>fileparameters</fromport>
- <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>fileparameters</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>FX</fromport>
- <tonode>Gradient</tonode> <toport>seq_FX</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>FY</fromport>
- <tonode>Gradient</tonode> <toport>seq_FY</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>DIMS</fromport>
- <tonode>Gradient</tonode> <toport>seq_DIMS</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>DIAG</fromport>
- <tonode>Gradient</tonode> <toport>lst_DIAG</toport>
- </datalink>
- <datalink control="false">
- <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>iter</fromport>
- <tonode>Gradient</tonode> <toport>lst_iter</toport>
- </datalink>
- </bloc>
- <parameter>
- <tonode>H_linearization.Finite_differences_derivation.Elementary_calculation.ASTER</tonode><toport>X</toport>
- <value><array><data>
-<value><double>80000</double></value>
-<value><double>1000</double></value>
-<value><double>30</double></value>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>ASTER_ROOT</toport>
- <value><string>''</string></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>rcdir</toport>
- <value><string>''</string></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>debug</toport>
- <value><boolean>0</boolean></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>DISPLAY</toport>
- <value><string>:0.0</string></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>SOURCES_ROOT</toport>
- <value><string>.</string></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>export</toport>
- <value><string>''</string></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>parametres</toport>
- <value><objref><![CDATA[(lp1
-.]]></objref></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>calcul</toport>
- <value><objref><![CDATA[(lp1
-.]]></objref></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>experience</toport>
- <value><objref><![CDATA[(lp1
-.]]></objref></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Temporary_Parameters</tonode><toport>fileparameters</toport>
- <value><string>[]</string></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Input</tonode><toport>X</toport>
- <value><array><data>
-<value><double>80000</double></value>
-<value><double>1000</double></value>
-<value><double>30</double></value>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>H_linearization.Input</tonode><toport>dX</toport>
- <value><array><data>
-<value><double>0.001</double></value>
-<value><double>0.001</double></value>
-<value><double>0.0001</double></value>
-</data></array></value>
- </parameter>
- <presentation name="H_linearization.Finite_differences_derivation.Elementary_calculation.ASTER" x="396" y="80.5" width="227.5" height="429"/>
- <presentation name="H_linearization.Finite_differences_derivation" x="368.856" y="80.5" width="651.5" height="657.5"/>
- <presentation name="H_linearization.Finite_differences_derivation.Elementary_calculation.Perturbated_point_X" x="9" y="81" width="227.5" height="139"/>
- <presentation name="H_linearization.Gradient" x="1064.38" y="338.546" width="227.5" height="284"/>
- <presentation name="H_linearization.Finite_differences_derivation.Elementary_calculation" x="9" y="130" width="633" height="518.5"/>
- <presentation name="H_linearization.Temporary_Parameters" x="10.3499" y="364.775" width="258.5" height="371"/>
- <presentation name="H_linearization" x="10.06" y="80.75" width="1301.38" height="747"/>
- <presentation name="H_linearization.Input" x="9" y="83.5" width="227.5" height="197"/>
- <presentation name="__ROOT__" x="0" y="0" width="1320.44" height="836.75"/>
-</proc>
+++ /dev/null
-<?xml version='1.0' encoding='iso-8859-1' ?>
-<proc>
- <property name="DefaultStudyID" value="1"/>
- <type name="bool" kind="bool"/>
- <sequence name="boolvec" content="bool"/>
- <type name="double" kind="double"/>
- <sequence name="dblevec" content="double"/>
- <objref name="file" id="file"/>
- <type name="int" kind="int"/>
- <sequence name="intvec" content="int"/>
- <objref name="pyobj" id="python:obj:1.0"/>
- <type name="string" kind="string"/>
- <sequence name="stringvec" content="string"/>
- <container name="DefaultContainer">
- <property name="container_name" value="FactoryServer"/>
- <property name="hostname" value="localhost"/>
- </container>
- <inline name="BLUE par matrices">
- <function name="algorithm">
- <code><![CDATA[import sys, os
-sys.path.insert(0, "../../Sources/daCore")
-sys.path.insert(0, "../../ComposantAD/daCore")
-#sys.path.insert(0, os.path.join(os.environ["HOME"],"SALOME5/supplements_JPA/ComposantAD/daCore"))
-
-import numpy
-from AssimilationStudy import AssimilationStudy
-
-def algorithm(Yo, B, R, H, Xb):
- #
- # Remise en place des matrices
- # -------------------
- dimensionXb = len( Xb )
- dimensionYo = len( Yo )
- B = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
- R = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
- H = numpy.matrix( H, numpy.float ).reshape((dimensionYo,dimensionXb))
- #
- # Analyse
- # -------
- ADD = AssimilationStudy()
- ADD.setBackground (asVector = Xb )
- ADD.setBackgroundError (asCovariance = B )
- ADD.setObservation (asVector = Yo )
- ADD.setObservationError (asCovariance = R )
- ADD.setObservationOperator(asMatrix = H )
- #
- ADD.setAlgorithm(choice="Blue")
- #
- ADD.analyze()
- #
- Xa = ADD.get("Analysis").valueserie(0)
- Innovation = ADD.get("Innovation").valueserie(0)
- A = []
- #
- return Xa, A, Innovation
-]]></code>
- </function>
- <inport name="Yo" type="dblevec"/>
- <inport name="B" type="dblevec"/>
- <inport name="R" type="dblevec"/>
- <inport name="H" type="dblevec"/>
- <inport name="Xb" type="dblevec"/>
- <outport name="Xa" type="dblevec"/>
- <outport name="A" type="dblevec"/>
- <outport name="Innovation" type="dblevec"/>
- </inline>
- <inline name="3D-VAR par matrices">
- <function name="algorithm">
- <code><![CDATA[import sys, os
-sys.path.insert(0, "../../Sources/daCore")
-sys.path.insert(0, "../../ComposantAD/daCore")
-#sys.path.insert(0, os.path.join(os.environ["HOME"],"SALOME5/supplements_JPA/ComposantAD/daCore"))
-
-import numpy
-from AssimilationStudy import AssimilationStudy
-
-def algorithm(Yo, B, R, H, Xb):
- #
- # Remise en place des matrices
- # -------------------
- dimensionXb = len( Xb )
- dimensionYo = len( Yo )
- B = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
- R = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
- H = numpy.matrix( H, numpy.float ).reshape((dimensionYo,dimensionXb))
- #
- # Analyse
- # -------
- ADD = AssimilationStudy()
- ADD.setBackground (asVector = Xb )
- ADD.setBackgroundError (asCovariance = B )
- ADD.setObservation (asVector = Yo )
- ADD.setObservationError (asCovariance = R )
- ADD.setObservationOperator(asMatrix = H )
- #
- ADD.setAlgorithm(choice="3DVAR")
- #
- ADD.analyze()
- #
- Xa = ADD.get("Analysis").valueserie(0)
- Innovation = ADD.get("Innovation").valueserie(0)
- A = []
- #
- return Xa, A, Innovation
-]]></code>
- </function>
- <inport name="Yo" type="dblevec"/>
- <inport name="B" type="dblevec"/>
- <inport name="R" type="dblevec"/>
- <inport name="H" type="dblevec"/>
- <inport name="Xb" type="dblevec"/>
- <outport name="Xa" type="dblevec"/>
- <outport name="A" type="dblevec"/>
- <outport name="Innovation" type="dblevec"/>
- </inline>
- <inline name="3D-VAR par fonctions">
- <function name="algorithm">
- <code><![CDATA[import sys, os
-sys.path.insert(0, "../../Sources/daCore")
-sys.path.insert(0, "../../ComposantAD/daCore")
-#sys.path.insert(0, os.path.join(os.environ["HOME"],"SALOME5/supplements_JPA/ComposantAD/daCore"))
-
-import numpy
-from AssimilationStudy import AssimilationStudy
-
-def algorithm( Yo, B, R, FunctionH, TangentH, AdjointH, Xb, Bounds ):
- #
- # Remise en place des matrices
- # -------------------
- dimensionXb = len( Xb )
- dimensionYo = len( Yo )
- B = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
- R = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
- #
- # 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":TangentH,
- "Adjoint":AdjointH} )
- #
- ADD.setAlgorithm(choice="3DVAR")
- ADD.setAlgorithmParameters(asDico={
- "Minimizer":"LBFGSB",
- "Bounds" :Bounds,
- })
- #
- ADD.analyze()
- #
- Xa = ADD.get("Analysis").valueserie(0)
- Innovation = ADD.get("Innovation").valueserie(0)
- A = []
- #
- return Xa, A, Innovation
-]]></code>
- </function>
- <inport name="Yo" type="dblevec"/>
- <inport name="B" type="dblevec"/>
- <inport name="R" type="dblevec"/>
- <inport name="FunctionH" type="pyobj"/>
- <inport name="TangentH" type="pyobj"/>
- <inport name="AdjointH" type="pyobj"/>
- <inport name="Xb" type="dblevec"/>
- <inport name="Bounds" type="pyobj"/>
- <outport name="Xa" type="dblevec"/>
- <outport name="A" type="dblevec"/>
- <outport name="Innovation" type="dblevec"/>
- </inline>
- <parameter>
- <tonode>BLUE par matrices</tonode><toport>Yo</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>BLUE par matrices</tonode><toport>B</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>BLUE par matrices</tonode><toport>R</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>BLUE par matrices</tonode><toport>H</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>BLUE par matrices</tonode><toport>Xb</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par matrices</tonode><toport>Yo</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par matrices</tonode><toport>B</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par matrices</tonode><toport>R</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par matrices</tonode><toport>H</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par matrices</tonode><toport>Xb</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par fonctions</tonode><toport>Yo</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par fonctions</tonode><toport>B</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par fonctions</tonode><toport>R</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par fonctions</tonode><toport>Xb</toport>
- <value><array><data>
-</data></array></value>
- </parameter>
- <parameter>
- <tonode>3D-VAR par fonctions</tonode><toport>Bounds</toport>
- <value><objref>(lp1
-.</objref></value>
- </parameter>
- <presentation name="BLUE par matrices" x="9" y="80.5" width="227.5" height="226"/>
- <presentation name="3D-VAR par matrices" x="237.5" y="80.5" width="227.5" height="226"/>
- <presentation name="3D-VAR par fonctions" x="465.5" y="80.5" width="227.5" height="313"/>
- <presentation name="__ROOT__" x="0" y="0" width="702" height="402.5"/>
-</proc>
+++ /dev/null
-#
-# Copyright (C) 2008-2009 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2009 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Outil numérique de calcul de la variable de Fisher pour comparer les
- variances de 2 échantillons
-
- Ce calcul nécessite :
- - en input :
- - les deux vecteurs (comme liste, array ou matrix) d'échantillons
- dont on veut comparer la variance,
- - la tolérance
- - en output :
- - la p-value,
- - la valeur de la variable aléatoire,
- - la réponse au test ainsi que
- - le message qui interprete la reponse du test.
-"""
-__author__ = "Sophie RICCI - Juillet 2008"
-
-import numpy
-from scipy.stats import betai
-
-# ==============================================================================
-def ComputeFisher(vector1 = None, vector2 = None, tolerance = 0.05 ):
- """
- Outil numérique de calcul de la variable de Fisher pour comparer les
- variances de 2 échantillons
-
- Ce calcul nécessite :
- - en input : les deux vecteurs (comme liste, array ou matrix)
- d'échantillons dont on veut comparer la variance, la
- tolérance
- - en output : la p-value, la valeur de la variable aléatoire,
- la réponse au test ainsi que le message qui interprete
- la reponse du test.
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Fisher value value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- #
- # Calcul des variances des echantillons
- # -------------------------------------
- # où var est calculee comme : var = somme (xi -xmean)**2 /(n-1)
- n1 = V1.size
- n2 = V2.size
- var1 = V1.std() * V1.std()
- var2 = V2.std() * V2.std()
- if (var1 > var2):
- f = var1/var2
- df1 = n1-1
- df2 = n2-1
- else:
- f= var2/var1
- df1 = n2-1
- df2 = n1-1
- prob1= betai(0.5*df2,0.5*df1,float(df2)/float(df2+df1*f))
- prob2= (1. - betai(0.5*df1, 0.5*df2, float(df1)/float(df1+df2/f)))
- prob = prob1 + prob2
- #
- # Calcul de la p-value
- # --------------------
- areafisher = 100 * prob
- #
- # Test
- # ----
- message = "Il y a %.2f%s de chance de se tromper en refusant l'hypothèse d'égalité des variances des 2 échantillons (si <%.2f%s, on refuse effectivement l'égalité)"%(areafisher,"%",100.*tolerance,"%")
- if (areafisher < (100.*tolerance)) :
- answerTestFisher = False
- else:
- answerTestFisher = True
- # print "La reponse au test est", answerTestFisher
-
- return areafisher, f, answerTestFisher, message
-
-# ==============================================================================
-if __name__ == "__main__":
- print "\nAUTOTEST\n"
- #
- # Echantillons
- # ------------
- x1 = [-1., 0., 4., 2., -1., 3.]
- x2 = [-1., 0., 4., 2., -1., 3.]
- #
- # Appel du calcul
- # ---------------
- [aire, f, reponse, message] = ComputeFisher(
- vector1 = x1,
- vector2 = x2,
- tolerance = 0.05 )
- #
- print " aire.....:", aire
- print " f........:", f
- print " reponse..:", reponse
- print " message..:", message
- print
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
from numpy import random
from scipy import arange, asarray, stats
from scipy.stats import histogram2, chisquare, chisqprob, norm
+from scipy.stats import distributions
+import scipy, sys, time, os
import logging
# ==============================================================================
logging.debug("Valeur du Khi2=%s"%self.__Q)
return self.__Q
+ def ComputeValue_ks(self) :
+ """
+ Calcul de la valeur ks pour le test de kolmogorov smirnov
+ """
+ if self.__testHomogen :
+ kobs,ktheo=self.__kobs.tolist(),self.__kobsHomogen.tolist()
+ else :
+ kobs,ktheo=self.__kobs.tolist(),self.__ktheo.tolist()
+
+ self.sortobs=asarray([scipy.sum(kobs[:i]) for i in range(1,1+(scipy.shape(kobs)[0]))])
+ self.sortobs=self.sortobs/float(scipy.sum(kobs))
+ self.sortthe=self.__cdf(self.__subdiv)
+ DKS1=max(abs(self.sortthe-self.sortobs))
+ DKS2=max(abs(self.sortthe[1:]-self.sortobs[0:-1]))
+ self.__ks=max(DKS1,DKS2)
+ logging.debug("Valeur du ks test=%s"%self.__ks)
+ return self.__ks
+
def ComputeArea(self):
"""
Calcul de la p-value
self.__areakhi2 = 100 * chisqprob(self.__Q, self.__ddl)
return self.__areakhi2
+ def ComputeArea_ks(self):
+ """
+ Calcul de la p-value ks
+ """
+ self.__areaks = 100 * distributions.ksone.sf(self.__ks,self.__modalites)
+ return self.__areaks
+
def WriteMessage(self):
"""
Interpretation du test
"""
message = "Il y a %.2f%s de chance de se tromper en refusant l'adequation"%(self.__areakhi2,"%")
return message
+
+ def WriteMessage_ks(self):
+ """
+ Interpretation du test
+ """
+ message = "Il y a %.2f%s de chance de se tromper en refusant l'adequation"%(self.__areaks,"%")
+ return message
def WriteMessageHomogen(self):
"""
logging.debug("message %s"%message)
return classes, eftheo, efobs, valeurKhi2, areaKhi2, message
+def Compute_ks(
+ vectorV = None,
+ dx = 0.1,
+ SuppressEmptyClasses = True,
+ nbclasses = None
+ ):
+ """
+ Test du KS adequation entre un vecteur donne et une gaussienne theo de mean et std celles du vecteur
+ """
+ essai = StatspourTests( cdftheo=norm, pdftest = None, obs = vectorV, use_mean_std_exp=True,dxmin=dx, obsHomogen = None, nbclasses = nbclasses)
+ essai.MakeClasses()
+ essai.ComputeObs()
+ essai.ComputeTheo()
+ classes,eftheo, efobs = essai.Computepdfs()
+ essai.Computeddl()
+ valeur_ks= essai.ComputeValue_ks()
+ area_ks = essai.ComputeArea_ks()
+ message = essai.WriteMessage_ks()
+ logging.debug("message %s"%message)
+ return classes, eftheo, efobs, valeur_ks, area_ks, message
+
+
def ComputeKhi2_Homogen(
vectorV1 = None,
vectorV2 = None,
print '\n AUTODIAGNOSTIC \n'
#
numpy.random.seed(100)
-
+#-------------------------------------------------------------------------
# Test de verification d adequation entre une gaussienne et un tirage gaussien
print ''
print 'Test de verification d adequation entre une gaussienne centree normale et un tirage gaussien'
raise ValueError("The computation of the khisquare value is WRONG")
numpy.random.seed(2490)
-
+#-------------------------------------------------------------------------
# Test de verification d adequation entre une gaussienne et un vecteur donne
print ''
print 'Test de verification d adequation entre une gaussienne et un vecteur donne'
print "The computation of the khisquare value is OK"
else :
raise ValueError("The computation of the khisquare value is WRONG")
+#-------------------------------------------------------------------------
+ # Test de verification d adequation entre une gaussienne et un vecteur donne par KS
+ print ''
+ print 'Test de verification d adequation entre une gaussienne et un vecteur donne par KS'
+# V = random.normal(50.,1.5,1000)
+ classes, eftheo, efobs, valeur_ks, area_ks, message = Compute_ks(dx = 0.1, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
+ print ' valeur_ks=',valeur_ks
+ print ' area_ks=',area_ks
+ print ' ',message
+
+ if (numpy.abs(area_ks - 82.17)< 1.e-2) :
+ print "The computation of the ks value is OK"
+ else :
+ raise ValueError("The computation of the ks value is WRONG")
+#-------------------------------------------------------------------------
# Test de d homogeneite entre 2 vecteurs donnes
print ''
print 'Test d homogeneite entre 2 vecteurs donnes'
print "The computation of the khisquare value is OK"
else :
raise ValueError("The computation of the khisquare value is WRONG")
-
+#-------------------------------------------------------------------------
# Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 10000
print ''
print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 10000'
print "The computation of the khisquare value is OK"
else :
raise ValueError("The computation of the khisquare value is WRONG")
-
+#-------------------------------------------------------------------------
# Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 1000
print ''
print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 1000'
print "The computation of the khisquare value is OK"
else :
raise ValueError("The computation of the khisquare value is WRONG")
+#-------------------------------------------------------------------------
# Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 100
print ''
#-*-coding:iso-8859-1-*-
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
#
-# Copyright (C) 2008-2009 EDF R&D
+# Copyright (C) 2008-2010 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
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
+ test103_PlotVectors.py
-GEN_DATA_INST = context.py
+testsplateforme_DATA = ${DATA_INST}
-testsplateforme_DATA = ${DATA_INST} ${GEN_DATA_INST}
-
-EXTRA_DIST = ${DATA_INST} context.py.in
+EXTRA_DIST = ${DATA_INST}
--- /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"
+
+import numpy
+from daCore.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"
+
+import numpy
+from daCore.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"
+
+import numpy
+from daCore.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"
+
+import numpy
+from daCore.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"
+
+import numpy
+from daCore.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"
+
+import numpy
+from daCore.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__ = """
+ 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"
+
+import numpy
+from daCore.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"
+
+import numpy
+from daCore.AssimilationStudy import AssimilationStudy
+from daCore.Persistence import OneScalar
+
+#===============================================================================
+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 = 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 = 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__ = """
+ Analyse moindre carres sans ebauche
+"""
+__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Septembre 2008"
+
+import numpy
+from daCore.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"
+
+import sys ; sys.path.insert(0, "../../Sources/daCore")
+
+import numpy
+from daCore.AssimilationStudy import AssimilationStudy
+from daCore.Persistence import OneVector
+
+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 = 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 sur le 3D-VAR que si l'erreur est nulle, l'incrément
+ d'analyse est nul.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2009"
+
+import numpy
+from daCore.AssimilationStudy import AssimilationStudy
+
+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 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"
+
+import numpy
+from daCore.AssimilationStudy import AssimilationStudy
+
+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"
+
+import numpy
+from daCore.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"
+
+from daCore.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"
+
+from daCore.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()
