From 357f29863c7d674883b65836b15af671839e8f79 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Andr=C3=A9=20Ribes?= Date: Wed, 26 Jan 2011 15:40:09 +0100 Subject: [PATCH] =?utf8?q?Mise=20=C3=A0=20jour=20de=20la=20biblioth=C3=A8q?= =?utf8?q?ue?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit --- src/daComposant/daAlgorithms/3DVAR.py | 20 +- src/daComposant/daAlgorithms/Blue.py | 8 +- src/daComposant/daAlgorithms/EnsembleBlue.py | 2 +- src/daComposant/daAlgorithms/Kalman.py | 2 +- .../daAlgorithms/LinearLeastSquares.py | 2 +- src/daComposant/daCore/AssimilationStudy.py | 140 +++-- src/daComposant/daCore/BasicObjects.py | 25 +- src/daComposant/daCore/Persistence.py | 237 ++++++++- src/daComposant/daCore/version.py | 4 +- src/daComposant/daDiagnostics/ComputeBiais.py | 85 --- .../daDiagnostics/ComputeCostFunction.py | 140 ----- .../daDiagnostics/ComputeVariance.py | 86 --- .../daDiagnostics/GaussianAdequation.py | 154 ------ src/daComposant/daDiagnostics/HLinearity.py | 138 ----- .../daDiagnostics/HomogeneiteKhi2.py | 116 ----- src/daComposant/daDiagnostics/ReduceBiais.py | 107 ---- .../daDiagnostics/VarianceOrder.py | 126 ----- src/daComposant/daNumerics/ComputeKhi2.py | 492 ------------------ src/daComposant/daNumerics/ComputeStudent.py | 260 --------- 19 files changed, 375 insertions(+), 1769 deletions(-) delete mode 100644 src/daComposant/daDiagnostics/ComputeBiais.py delete mode 100644 src/daComposant/daDiagnostics/ComputeCostFunction.py delete mode 100644 src/daComposant/daDiagnostics/ComputeVariance.py delete mode 100644 src/daComposant/daDiagnostics/GaussianAdequation.py delete mode 100644 src/daComposant/daDiagnostics/HLinearity.py delete mode 100644 src/daComposant/daDiagnostics/HomogeneiteKhi2.py delete mode 100644 src/daComposant/daDiagnostics/ReduceBiais.py delete mode 100644 src/daComposant/daDiagnostics/VarianceOrder.py delete mode 100644 src/daComposant/daNumerics/ComputeKhi2.py delete mode 100644 src/daComposant/daNumerics/ComputeStudent.py diff --git a/src/daComposant/daAlgorithms/3DVAR.py b/src/daComposant/daAlgorithms/3DVAR.py index fca8d40..d42ce65 100644 --- a/src/daComposant/daAlgorithms/3DVAR.py +++ b/src/daComposant/daAlgorithms/3DVAR.py @@ -46,7 +46,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): self._name = "3DVAR" logging.debug("%s Initialisation"%self._name) - def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None): + def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None): """ Calcul de l'estimateur 3D-VAR """ @@ -64,6 +64,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): else: logging.debug("%s Calcul de Hm(Xb)"%self._name) HXb = Hm( Xb ) + HXb = numpy.asmatrix(HXb).flatten().T # # Calcul du préconditionnement # ---------------------------- @@ -92,6 +93,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): logging.info("%s CostFunction Jb = %s"%(self._name, Jb)) logging.info("%s CostFunction Jo = %s"%(self._name, Jo)) logging.info("%s CostFunction J = %s"%(self._name, J)) + self.StoredVariables["CurrentState"].store( _X.A1 ) self.StoredVariables["CostFunctionJb"].store( Jb ) self.StoredVariables["CostFunctionJo"].store( Jo ) self.StoredVariables["CostFunctionJ" ].store( J ) @@ -108,9 +110,6 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): logging.debug("%s GradientOfCostFunction GradJb = %s"%(self._name, numpy.asmatrix( GradJb ).flatten())) logging.debug("%s GradientOfCostFunction GradJo = %s"%(self._name, numpy.asmatrix( GradJo ).flatten())) logging.debug("%s GradientOfCostFunction GradJ = %s"%(self._name, numpy.asmatrix( GradJ ).flatten())) - # self.StoredVariables["GradientOfCostFunctionJb"].store( Jb ) - # self.StoredVariables["GradientOfCostFunctionJo"].store( Jo ) - # self.StoredVariables["GradientOfCostFunctionJ" ].store( J ) return GradJ.A1 # # Point de démarrage de l'optimisation : Xini = Xb @@ -123,18 +122,18 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): # # Paramètres de pilotage # ---------------------- - if Par.has_key("Bounds") and (type(Par["Bounds"]) is type([]) or type(Par["Bounds"]) is type(())) and (len(Par["Bounds"]) > 0): - Bounds = Par["Bounds"] + if Parameters.has_key("Bounds") and (type(Parameters["Bounds"]) is type([]) or type(Parameters["Bounds"]) is type(())) and (len(Parameters["Bounds"]) > 0): + Bounds = Parameters["Bounds"] else: Bounds = None MinimizerList = ["LBFGSB","TNC", "CG", "BFGS"] - if Par.has_key("Minimizer") and (Par["Minimizer"] in MinimizerList): - Minimizer = str( Par["Minimizer"] ) + if Parameters.has_key("Minimizer") and (Parameters["Minimizer"] in MinimizerList): + Minimizer = str( Parameters["Minimizer"] ) else: Minimizer = "LBFGSB" logging.debug("%s Minimiseur utilisé = %s"%(self._name, Minimizer)) - if Par.has_key("MaximumNumberOfSteps") and (Par["MaximumNumberOfSteps"] > -1): - maxiter = int( Par["MaximumNumberOfSteps"] ) + if Parameters.has_key("MaximumNumberOfSteps") and (Parameters["MaximumNumberOfSteps"] > -1): + maxiter = int( Parameters["MaximumNumberOfSteps"] ) else: maxiter = 15000 logging.debug("%s Nombre maximal de pas d'optimisation = %s"%(self._name, maxiter)) @@ -150,6 +149,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): bounds = Bounds, maxfun = maxiter, iprint = iprint, + factr = 1., ) logging.debug("%s %s Minimum = %s"%(self._name, Minimizer, Minimum)) logging.debug("%s %s Nb of F = %s"%(self._name, Minimizer, Informations['funcalls'])) diff --git a/src/daComposant/daAlgorithms/Blue.py b/src/daComposant/daAlgorithms/Blue.py index af7f44f..04efa6f 100644 --- a/src/daComposant/daAlgorithms/Blue.py +++ b/src/daComposant/daAlgorithms/Blue.py @@ -27,6 +27,8 @@ import logging from daCore import BasicObjects, PlatformInfo m = PlatformInfo.SystemUsage() +import numpy + # ============================================================================== class ElementaryAlgorithm(BasicObjects.Algorithm): def __init__(self): @@ -34,7 +36,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): self._name = "BLUE" logging.debug("%s Initialisation"%self._name) - def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None): + def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None): """ Calcul de l'estimateur BLUE (ou Kalman simple, ou Interpolation Optimale) """ @@ -52,8 +54,10 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): else: logging.debug("%s Calcul de Hm * Xb"%self._name) HXb = Hm * Xb - + HXb = numpy.asmatrix(HXb).flatten().T + # # Calcul de la matrice de gain dans l'espace le plus petit + # -------------------------------------------------------- if Y.size <= Xb.size: logging.debug("%s Calcul de K dans l'espace des observations"%self._name) K = B * Ht * (Hm * B * Ht + R).I diff --git a/src/daComposant/daAlgorithms/EnsembleBlue.py b/src/daComposant/daAlgorithms/EnsembleBlue.py index 8e27046..f09c2a0 100644 --- a/src/daComposant/daAlgorithms/EnsembleBlue.py +++ b/src/daComposant/daAlgorithms/EnsembleBlue.py @@ -36,7 +36,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): self._name = "ENSEMBLEBLUE" logging.debug("%s Initialisation"%self._name) - def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None ): + def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None ): """ Calcul d'une estimation BLUE d'ensemble : - génération d'un ensemble d'observations, de même taille que le diff --git a/src/daComposant/daAlgorithms/Kalman.py b/src/daComposant/daAlgorithms/Kalman.py index 06875d2..cccfd35 100644 --- a/src/daComposant/daAlgorithms/Kalman.py +++ b/src/daComposant/daAlgorithms/Kalman.py @@ -39,7 +39,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): self._name = "KALMAN" logging.debug("%s Initialisation"%self._name) - def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None): + def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None): """ Calcul de l'estimateur de Kalman """ diff --git a/src/daComposant/daAlgorithms/LinearLeastSquares.py b/src/daComposant/daAlgorithms/LinearLeastSquares.py index a4470be..7d6fa98 100644 --- a/src/daComposant/daAlgorithms/LinearLeastSquares.py +++ b/src/daComposant/daAlgorithms/LinearLeastSquares.py @@ -33,7 +33,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): 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): + def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None): """ Calcul de l'estimateur au sens des moindres carres sans ebauche """ diff --git a/src/daComposant/daCore/AssimilationStudy.py b/src/daComposant/daCore/AssimilationStudy.py index f7b1cd2..5587dd5 100644 --- a/src/daComposant/daCore/AssimilationStudy.py +++ b/src/daComposant/daCore/AssimilationStudy.py @@ -71,14 +71,15 @@ class AssimilationStudy: self.__M = {} # self.__X = Persistence.OneVector() - self.__Parameters = {} + self.__Parameters = {} self.__StoredDiagnostics = {} + self.__StoredInputs = {} # # Variables temporaires - self.__algorithm = {} - self.__algorithmFile = None - self.__algorithmName = None - self.__diagnosticFile = None + self.__algorithm = {} + self.__algorithmFile = None + self.__algorithmName = None + self.__diagnosticFile = None # # Récupère le chemin du répertoire parent et l'ajoute au path # (Cela complète l'action de la classe PathManagement dans PlatformInfo, @@ -92,6 +93,7 @@ class AssimilationStudy: asVector = None, asPersistentVector = None, Scheduler = None, + toBeStored = False, ): """ Permet de définir l'estimation a priori : @@ -100,6 +102,8 @@ class AssimilationStudy: - asPersistentVector : entrée des données, comme un vecteur de type persistent contruit avec la classe ad-hoc "Persistence" - Scheduler est le contrôle temporel des données + - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour + être rendue disponible au même titre que les variables de calcul """ if asVector is not None: if type( asVector ) is type( numpy.matrix([]) ): @@ -110,15 +114,24 @@ class AssimilationStudy: self.__Xb = asPersistentVector else: raise ValueError("Error: improperly defined background") + if toBeStored: + self.__StoredInputs["Background"] = self.__Xb return 0 - def setBackgroundError(self, asCovariance=None): + def setBackgroundError(self, + asCovariance = None, + toBeStored = False, + ): """ Permet de définir la covariance des erreurs d'ébauche : - asCovariance : entrée des données, comme une matrice compatible avec le constructeur de numpy.matrix + - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour + être rendue disponible au même titre que les variables de calcul """ self.__B = numpy.matrix( asCovariance, numpy.float ) + if toBeStored: + self.__StoredInputs["BackgroundError"] = self.__B return 0 # ----------------------------------------------------------- @@ -126,6 +139,7 @@ class AssimilationStudy: asVector = None, asPersistentVector = None, Scheduler = None, + toBeStored = False, ): """ Permet de définir les observations : @@ -134,6 +148,8 @@ class AssimilationStudy: - asPersistentVector : entrée des données, comme un vecteur de type persistent contruit avec la classe ad-hoc "Persistence" - Scheduler est le contrôle temporel des données disponibles + - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour + être rendue disponible au même titre que les variables de calcul """ if asVector is not None: if type( asVector ) is type( numpy.matrix([]) ): @@ -144,21 +160,31 @@ class AssimilationStudy: self.__Y = asPersistentVector else: raise ValueError("Error: improperly defined observations") + if toBeStored: + self.__StoredInputs["Observation"] = self.__Y return 0 - def setObservationError(self, asCovariance=None): + def setObservationError(self, + asCovariance = None, + toBeStored = False, + ): """ Permet de définir la covariance des erreurs d'observations : - asCovariance : entrée des données, comme une matrice compatible avec le constructeur de numpy.matrix + - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour + être rendue disponible au même titre que les variables de calcul """ self.__R = numpy.matrix( asCovariance, numpy.float ) + if toBeStored: + self.__StoredInputs["ObservationError"] = self.__R return 0 def setObservationOperator(self, asFunction = {"Direct":None, "Tangent":None, "Adjoint":None}, asMatrix = None, appliedToX = None, + toBeStored = False, ): """ Permet de définir un opérateur d'observation H. L'ordre de priorité des @@ -175,6 +201,8 @@ class AssimilationStudy: X divers, l'opérateur par sa valeur appliquée à cet X particulier, sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom. L'opérateur doit néanmoins déjà avoir été défini comme d'habitude. + - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour + être rendue disponible au même titre que les variables de calcul """ if (type(asFunction) is type({})) and (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None): if not asFunction.has_key("Direct") or (asFunction["Direct"] is None): @@ -207,6 +235,8 @@ class AssimilationStudy: else: self.__H["AppliedToX"] = None # + if toBeStored: + self.__StoredInputs["ObservationOperator"] = self.__H return 0 # ----------------------------------------------------------- @@ -214,6 +244,7 @@ class AssimilationStudy: asFunction = {"Direct":None, "Tangent":None, "Adjoint":None}, asMatrix = None, Scheduler = None, + toBeStored = False, ): """ Permet de définir un opérateur d'évolution M. L'ordre de priorité des @@ -226,6 +257,8 @@ class AssimilationStudy: la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La matrice fournie doit être sous une forme compatible avec le constructeur de numpy.matrix. + - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour + être rendue disponible au même titre que les variables de calcul """ if (type(asFunction) is type({})) and (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None): if not asFunction.has_key("Direct") or (asFunction["Direct"] is None): @@ -241,19 +274,32 @@ class AssimilationStudy: self.__M["Adjoint"] = Operator( fromMatrix = matrice.T ) else: raise ValueError("Error: improperly defined evolution operator") + # + if toBeStored: + self.__StoredInputs["EvolutionModel"] = self.__M return 0 - def setEvolutionError(self, asCovariance=None): + def setEvolutionError(self, + asCovariance = None, + toBeStored = False, + ): """ Permet de définir la covariance des erreurs de modèle : - asCovariance : entrée des données, comme une matrice compatible avec le constructeur de numpy.matrix + - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour + être rendue disponible au même titre que les variables de calcul """ self.__Q = numpy.matrix( asCovariance, numpy.float ) + if toBeStored: + self.__StoredInputs["EvolutionError"] = self.__Q return 0 # ----------------------------------------------------------- - def setControls (self, asVector = None ): + def setControls (self, + asVector = None, + toBeStored = False, + ): """ Permet de définir la valeur initiale du vecteur X contenant toutes les variables de contrôle, i.e. les paramètres ou l'état dont on veut @@ -261,9 +307,13 @@ class AssimilationStudy: algorithme itératif/incrémental - asVector : entrée des données, comme un vecteur compatible avec le constructeur de numpy.matrix. + - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour + être rendue disponible au même titre que les variables de calcul """ if asVector is not None: self.__X.store( asVector ) + if toBeStored: + self.__StoredInputs["Controls"] = self.__X return 0 # ----------------------------------------------------------- @@ -302,6 +352,7 @@ class AssimilationStudy: # Instancie un objet du type élémentaire du fichier # ------------------------------------------------- self.__algorithm = self.__algorithmFile.ElementaryAlgorithm() + self.__StoredInputs["AlgorithmName"] = str(choice) return 0 def setAlgorithmParameters(self, asDico=None): @@ -309,7 +360,11 @@ class AssimilationStudy: Permet de définir les paramètres de l'algorithme, sous la forme d'un dictionnaire. """ - self.__Parameters = dict( asDico ) + if asDico is not None: + self.__Parameters = dict( asDico ) + else: + self.__Parameters = {} + self.__StoredInputs["AlgorithmParameters"] = self.__Parameters return 0 # ----------------------------------------------------------- @@ -341,8 +396,10 @@ class AssimilationStudy: # # Instancie un objet du type élémentaire du fichier # ------------------------------------------------- - if self.__StoredDiagnostics.has_key(name): - raise ValueError("A diagnostic with the same name already exists") + if self.__StoredInputs.has_key(name): + raise ValueError("A default input with the same name \"%s\" already exists."%str(name)) + elif self.__StoredDiagnostics.has_key(name): + raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name)) else: self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic( name = name, @@ -450,14 +507,14 @@ class AssimilationStudy: self.shape_validate() # self.__algorithm.run( - Xb = self.__Xb, - Y = self.__Y, - H = self.__H, - M = self.__M, - R = self.__R, - B = self.__B, - Q = self.__Q, - Par = self.__Parameters, + Xb = self.__Xb, + Y = self.__Y, + H = self.__H, + M = self.__M, + R = self.__R, + B = self.__B, + Q = self.__Q, + Parameters = self.__Parameters, ) return 0 @@ -472,31 +529,33 @@ class AssimilationStudy: if key is not None: if self.__algorithm.has_key(key): return self.__algorithm.get( key ) + elif self.__StoredInputs.has_key(key): + return self.__StoredInputs[key] elif self.__StoredDiagnostics.has_key(key): return self.__StoredDiagnostics[key] else: - raise ValueError("The requested key \"%s\" does not exists as a diagnostic or as a stored variable."%key) + raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key) else: allvariables = self.__algorithm.get() allvariables.update( self.__StoredDiagnostics ) + allvariables.update( self.__StoredInputs ) return allvariables - def get_available_variables(self, key=None): + def get_available_variables(self): """ Renvoie les variables potentiellement utilisables pour l'étude, + initialement stockées comme données d'entrées ou dans les algorithmes, 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( {} ): + if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0: 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 = [] + if len( self.__algorithm.keys()) > 0: + variables.extend( self.__algorithm.get().keys() ) + if len( self.__StoredInputs.keys() ) > 0: + variables.extend( self.__StoredInputs.keys() ) variables.sort() return variables @@ -616,10 +675,11 @@ class AssimilationStudy: HookFunction = HookFunction, HookParameters = HookParameters, ) + def prepare_to_pickle(self): - self.__algorithmFile = None - self.__diagnosticFile = None - self.__H = {} + self.__algorithmFile = None + self.__diagnosticFile = None + self.__H = {} # ============================================================================== if __name__ == "__main__": @@ -646,19 +706,21 @@ if __name__ == "__main__": print "Innovation :", ADD.get("Innovation").valueserie(0) print - print "Algorithmes disponibles :", ADD.get_available_algorithms() - # print " Chemin des algorithmes :", ADD.get_algorithms_main_path() - print "Diagnostics disponibles :", ADD.get_available_diagnostics() - # print " Chemin des diagnostics :", ADD.get_diagnostics_main_path() + print "Algorithmes disponibles.......................:", ADD.get_available_algorithms() + # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path() + print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics() + # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path() + print "Variables disponibles.........................:", ADD.get_available_variables() print ADD.setDiagnostic("RMS", "Ma RMS") liste = ADD.get().keys() liste.sort() - print "Variables et diagnostics disponibles :", liste - print + print "Variables et diagnostics nommés disponibles...:", liste + print + print "Exemple de mise en place d'un observeur :" def obs(var=None,info=None): print " ---> Mise en oeuvre de l'observer" print " var =",var.valueserie(-1) diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index eee46b5..df117ad 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -121,6 +121,13 @@ class Algorithm: interne à l'objet, mais auquel on accède par la méthode "get". Les variables prévues sont : + - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes + - CostFunctionJb : partie ébauche ou background de la fonction-cout + - CostFunctionJo : partie observations de la fonction-cout + - GradientOfCostFunctionJ : gradient de la fonction-cout globale + - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout + - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout + - CurrentState : état courant lors d'itérations - Analysis : l'analyse - Innovation : l'innovation : d = Y - H Xb - SigmaObs2 : correction optimale des erreurs d'observation @@ -132,13 +139,16 @@ class Algorithm: On peut rajouter des variables à stocker dans l'initialisation de l'algorithme élémentaire qui va hériter de cette classe """ + self._name = None self.StoredVariables = {} + # self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ") self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb") self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo") - self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneScalar(name = "GradientOfCostFunctionJ") - self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneScalar(name = "GradientOfCostFunctionJb") - self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneScalar(name = "GradientOfCostFunctionJo") + self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ") + self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb") + self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo") + self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState") self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis") self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation") self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2") @@ -147,7 +157,6 @@ class Algorithm: self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB") self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA") self.StoredVariables["CovarianceAPosteriori"] = Persistence.OneMatrix(name = "CovarianceAPosteriori") - self._name = None def get(self, key=None): """ @@ -168,7 +177,13 @@ class Algorithm: """ return self.StoredVariables.has_key(key) - def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None): + def keys(self): + """ + Renvoie la liste des clés de variables stockées. + """ + return self.StoredVariables.keys() + + def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None): """ Doit implémenter l'opération élémentaire de calcul d'assimilation sous sa forme mathématique la plus naturelle possible. diff --git a/src/daComposant/daCore/Persistence.py b/src/daComposant/daCore/Persistence.py index 3ac7425..79dae99 100644 --- a/src/daComposant/daCore/Persistence.py +++ b/src/daComposant/daCore/Persistence.py @@ -153,6 +153,20 @@ class Persistence: """ return len(self.__steps) + # --------------------------------------------------------- + # Méthodes d'accès de type dictionnaire + def keys(self): + return self.stepserie() + + def values(self): + return self.valueserie() + + def items(self): + pairs = [] + for i in xrange(self.stepnumber()): + pairs.append( (self.stepserie(item=i), self.valueserie(item=i)) ) + return pairs + # --------------------------------------------------------- def mean(self): """ @@ -564,6 +578,180 @@ class OneList(Persistence): def __init__(self, name="", unit="", basetype = list): Persistence.__init__(self, name, unit, basetype) +def NoType( value ): return value + +class OneNoType(Persistence): + """ + Classe définissant le stockage d'un objet sans modification (cast) de type. + Attention, selon le véritable type de l'objet stocké à chaque pas, les + opérations arithmétiques à base de numpy peuvent être invalides ou donner + des résultats inatendus. Cette classe n'est donc à utiliser qu'à bon escient + volontairement, et pas du tout par défaut. + """ + def __init__(self, name="", unit="", basetype = NoType): + Persistence.__init__(self, name, unit, basetype) + +# ============================================================================== +class CompositePersistence: + """ + Structure de stockage permettant de rassembler plusieurs objets de + persistence. + + Des objets par défaut sont prévus, et des objets supplémentaires peuvent + être ajoutés. + """ + def __init__(self, name=""): + """ + name : nom courant + + La gestion interne des données est exclusivement basée sur les variables + initialisées ici (qui ne sont pas accessibles depuis l'extérieur des + objets comme des attributs) : + __StoredObjects : objets de type persistence collectés dans cet objet + """ + self.__name = str(name) + # + self.__StoredObjects = {} + # + # Definition des objets par defaut + # -------------------------------- + self.__StoredObjects["Background"] = OneVector("Background", basetype=numpy.array) + self.__StoredObjects["BackgroundError"] = OneMatrix("BackgroundError") + self.__StoredObjects["Observation"] = OneVector("Observation", basetype=numpy.array) + self.__StoredObjects["ObservationError"] = OneMatrix("ObservationError") + self.__StoredObjects["Analysis"] = OneVector("Analysis", basetype=numpy.array) + self.__StoredObjects["AnalysisError"] = OneMatrix("AnalysisError") + self.__StoredObjects["Innovation"] = OneVector("Innovation", basetype=numpy.array) + self.__StoredObjects["KalmanGainK"] = OneMatrix("KalmanGainK") + self.__StoredObjects["OperatorH"] = OneMatrix("OperatorH") + self.__StoredObjects["RmsOMA"] = OneScalar("RmsOMA") + self.__StoredObjects["RmsOMB"] = OneScalar("RmsOMB") + self.__StoredObjects["RmsBMA"] = OneScalar("RmsBMA") + # + + def store(self, name=None, value=None, step=None): + """ + Stockage d'une valeur "value" pour le "step" dans la variable "name". + """ + if name is None: raise ValueError("Storable object name is required for storage.") + if name not in self.__StoredObjects.keys(): + raise ValueError("No such name '%s' exists in storable objects."%name) + self.__StoredObjects[name].store( value=value, step=step ) + + def add_object(self, name=None, persistenceType=Persistence, basetype=numpy.array ): + """ + Ajoute dans les objets stockables un nouvel objet défini par son nom, son + type de Persistence et son type de base à chaque pas. + """ + if name is None: raise ValueError("Object name is required for adding an object.") + if name in self.__StoredObjects.keys(): + raise ValueError("An object with the same name '%s' already exists in storable objects. Choose another one."%name) + self.__StoredObjects[name] = persistenceType( name=str(name), basetype=basetype ) + + def get_object(self, name=None ): + """ + Renvoie l'objet de type Persistence qui porte le nom demandé. + """ + if name is None: raise ValueError("Object name is required for retrieving an object.") + if name not in self.__StoredObjects.keys(): + raise ValueError("No such name '%s' exists in stored objects."%name) + return self.__StoredObjects[name] + + def set_object(self, name=None, objet=None ): + """ + Affecte directement un 'objet' qui porte le nom 'name' demandé. + Attention, il n'est pas effectué de vérification sur le type, qui doit + comporter les méthodes habituelles de Persistence pour que cela + fonctionne. + """ + if name is None: raise ValueError("Object name is required for setting an object.") + if name in self.__StoredObjects.keys(): + raise ValueError("An object with the same name '%s' already exists in storable objects. Choose another one."%name) + self.__StoredObjects[name] = objet + + def del_object(self, name=None ): + """ + Supprime un objet de la liste des objets stockables. + """ + if name is None: raise ValueError("Object name is required for retrieving an object.") + if name not in self.__StoredObjects.keys(): + raise ValueError("No such name '%s' exists in stored objects."%name) + del self.__StoredObjects[name] + + # --------------------------------------------------------- + # Méthodes d'accès de type dictionnaire + def __getitem__(self, name=None ): + return self.get_object( name ) + + def __setitem__(self, name=None, objet=None ): + self.set_object( name, objet ) + + def keys(self): + return self.get_stored_objects(hideVoidObjects = False) + + def values(self): + return self.__StoredObjects.values() + + def items(self): + return self.__StoredObjects.items() + + # --------------------------------------------------------- + def get_stored_objects(self, hideVoidObjects = False): + objs = self.__StoredObjects.keys() + if hideVoidObjects: + usedObjs = [] + for k in objs: + try: + if len(self.__StoredObjects[k]) > 0: usedObjs.append( k ) + except: + pass + objs = usedObjs + objs.sort() + return objs + + # --------------------------------------------------------- + def save_composite(self, filename=None, mode="pickle"): + """ + Enregistre l'objet dans le fichier indiqué selon le "mode" demandé, + et renvoi le nom du fichier + """ + import os + if filename is None: + filename = os.tempnam( os.getcwd(), 'dacp' ) + ".pkl" + else: + filename = os.path.abspath( filename ) + # + import cPickle + if mode is "pickle": + output = open( filename, 'wb') + cPickle.dump(self, output) + output.close() + else: + raise ValueError("Save mode '%s' unknown. Choose another one."%mode) + # + return filename + + def load_composite(self, filename=None, mode="pickle"): + """ + Recharge un objet composite sauvé en fichier + """ + import os + if filename is None: + raise ValueError("A file name if requested to load a composite.") + else: + filename = os.path.abspath( filename ) + # + import cPickle + if mode is "pickle": + pkl_file = open(filename, 'rb') + output = cPickle.load(pkl_file) + for k in output.keys(): + self[k] = output[k] + else: + raise ValueError("Load mode '%s' unknown. Choose another one."%mode) + # + return filename + # ============================================================================== if __name__ == "__main__": print '\n AUTODIAGNOSTIC \n' @@ -704,7 +892,48 @@ if __name__ == "__main__": del OBJET_DE_TEST print - print "======> Affichage d'objets stockés" + print "======> Utilisation des méthodes d'accès de type dictionnaire" + OBJET_DE_TEST = OneScalar("My int", unit="cm", basetype=int) + for i in range(5): + OBJET_DE_TEST.store( 7+i ) + print "Taille \"len\" :", len(OBJET_DE_TEST) + print "Les pas de stockage :", OBJET_DE_TEST.keys() + print "Les valeurs :", OBJET_DE_TEST.values() + print "Les paires :", OBJET_DE_TEST.items() + del OBJET_DE_TEST + print + + print "======> Persistence composite" + OBJET_DE_TEST = CompositePersistence("My CompositePersistence") + print "Objets stockables :", OBJET_DE_TEST.get_stored_objects() + print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True ) + print "--> Stockage d'une valeur de Background" + OBJET_DE_TEST.store("Background",numpy.zeros(5)) + print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True ) + print "--> Ajout d'un objet nouveau par defaut, de type vecteur numpy par pas" + OBJET_DE_TEST.add_object("ValeursVectorielles") + OBJET_DE_TEST.store("ValeursVectorielles",numpy.zeros(5)) + print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True ) + print "--> Ajout d'un objet nouveau de type liste par pas" + OBJET_DE_TEST.add_object("ValeursList", persistenceType=OneList ) + OBJET_DE_TEST.store("ValeursList",range(5)) + print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True ) + print "--> Ajout d'un objet nouveau, de type vecteur string par pas" + OBJET_DE_TEST.add_object("ValeursStr", persistenceType=Persistence, basetype=str ) + OBJET_DE_TEST.store("ValeursStr","c020") + OBJET_DE_TEST.store("ValeursStr","c021") + print "Les valeurs :", OBJET_DE_TEST.get_object("ValeursStr").valueserie() + print "Acces comme dict :", OBJET_DE_TEST["ValeursStr"].stepserie() + print "Acces comme dict :", OBJET_DE_TEST["ValeursStr"].valueserie() + print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True ) + print "--> Suppression d'un objet" + OBJET_DE_TEST.del_object("ValeursVectorielles") + print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True ) + print "--> Enregistrement de l'objet complet de Persistence composite" + OBJET_DE_TEST.save_composite("composite.pkl") + print + + print "======> Affichage graphique d'objets stockés" OBJET_DE_TEST = Persistence("My object", unit="", basetype=numpy.array) D = OBJET_DE_TEST vect1 = [1, 2, 1, 2, 1] @@ -715,7 +944,7 @@ if __name__ == "__main__": D.store(vect1) D.store(vect2) D.store(vect3) - print "Affichage de l'ensemble du stockage sur une même image" + print "Affichage graphique de l'ensemble du stockage sur une même image" D.stepplot( title = "Tous les vecteurs", filename="vecteurs.ps", @@ -724,7 +953,7 @@ if __name__ == "__main__": pause = False ) print "Stockage d'un quatrième vecteur de longueur différente" D.store(vect4) - print "Affichage séparé du dernier stockage" + print "Affichage graphique séparé du dernier stockage" D.plot( item = 3, title = "Vecteurs", @@ -738,7 +967,7 @@ if __name__ == "__main__": del OBJET_DE_TEST print - print "======> Affichage dynamique d'objets" + print "======> Affichage graphique dynamique d'objets" OBJET_DE_TEST = Persistence("My object", unit="", basetype=float) D = OBJET_DE_TEST D.plot( diff --git a/src/daComposant/daCore/version.py b/src/daComposant/daCore/version.py index 8bb6b99..487717a 100644 --- a/src/daComposant/daCore/version.py +++ b/src/daComposant/daCore/version.py @@ -19,5 +19,5 @@ # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com # name = "Data Assimilation Package" -version = "0.2.0" -date = "lundi 23 septembre 2009, 11:11:11 (UTC+0200)" +version = "0.3.0" +date = "lundi 15 octobre 2010, 11:11:11 (UTC+0200)" diff --git a/src/daComposant/daDiagnostics/ComputeBiais.py b/src/daComposant/daDiagnostics/ComputeBiais.py deleted file mode 100644 index 5ee216f..0000000 --- a/src/daComposant/daDiagnostics/ComputeBiais.py +++ /dev/null @@ -1,85 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Calcul du biais (i.e. la moyenne) à chaque pas. Ce diagnostic très simple - est présent pour rappeller à l'utilisateur de l'assimilation qu'il faut - qu'il vérifie le biais de ses erreurs en particulier. -""" -__author__ = "Sophie RICCI - Aout 2008" - -import numpy -from daCore import BasicObjects, Persistence - -# ============================================================================== -class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar): - def __init__(self, name = "", unit = "", basetype = None, parameters = {}): - BasicObjects.Diagnostic.__init__(self, name, parameters) - Persistence.OneScalar.__init__( self, name, unit, basetype = float ) - - def _formula(self, V): - """ - Calcul du biais, qui est simplement la moyenne du vecteur - """ - biais = V.mean() - # - return biais - - def calculate(self, vector = None, step = None): - """ - Teste les arguments, active la formule de calcul et stocke le résultat - """ - if vector is None: - raise ValueError("One vector must be given to compute biais") - V = numpy.array(vector) - if V.size < 1: - raise ValueError("The given vector must not be empty") - # - value = self._formula( V) - # - self.store( value = value, step = step ) - -#=============================================================================== -if __name__ == "__main__": - print '\n AUTODIAGNOSTIC \n' - # - # Instanciation de l'objet diagnostic - # ----------------------------------- - D = ElementaryDiagnostic("Mon ComputeBiais") - # - # Tirage d un vecteur choisi - # -------------------------- - x = numpy.matrix(([3., 4., 5.])) - print " Le vecteur de type 'matrix' choisi est..:", x - print " Le biais attendu de ce vecteur est......:", x.mean() - # - D.calculate( vector = x) - print " Le biais obtenu de ce vecteur est.......:", D.valueserie(0) - print - # - # Tirage d un vecteur choisi - # -------------------------- - x = numpy.array(range(11)) - print " Le vecteur de type 'array' choisi est...:", x - print " Le biais attendu de ce vecteur est......:", x.mean() - # - D.calculate( vector = x) - print " Le biais obtenu de ce vecteur est.......:", D.valueserie(1) - print diff --git a/src/daComposant/daDiagnostics/ComputeCostFunction.py b/src/daComposant/daDiagnostics/ComputeCostFunction.py deleted file mode 100644 index 2db1acc..0000000 --- a/src/daComposant/daDiagnostics/ComputeCostFunction.py +++ /dev/null @@ -1,140 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Calcul de la fonction coût -""" -__author__ = "Sophie RICCI - Octobre 2008" - -import numpy -from daCore import BasicObjects, Persistence -import logging - -# ============================================================================== -class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar): - def __init__(self, name = "", unit = "", basetype = None, parameters = {}): - 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. * 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. * numpy.dot((Y - HX) ,numpy.asarray(numpy.dot(R.I,(Y - HX)).A1)) - 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, Hx = None, xb = 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) : - raise ValueError("Vectors x, xb and yo must be given to compute J") -# if (type(x) is not float) and (type(x) is not numpy.float64) : -# if (x.size < 1) or (xb.size < 1) or (yo.size < 1): -# raise ValueError("Vectors x, xb and yo must not be empty") - 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 (Hx is None ) : - raise ValueError("The given vector must be given") -# if (Hx.size < 1) : -# raise ValueError("The given vector must not be empty") - HX = Hx.A1 - if (B is None ) or (R is None ): - raise ValueError("The matrices B and R must be given") -# if (B.size < 1) or (R.size < 1) : -# raise ValueError("The matrices B and R must not be empty") - # - value = self._formula(X, HX, Xb, 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.]) - xb = numpy.array([2., 2.]) - yo = numpy.array([5., 6.]) - H = numpy.matrix(numpy.identity(2)) -# Hx = H*x - Hx = numpy.dot(H,x) - Hx = Hx.T - B = numpy.matrix(numpy.identity(2)) - R = numpy.matrix(numpy.identity(2)) - # - D.calculate( x = x, Hx = Hx, xb = xb, 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) - print "La fonction cout J vaut ...: ",D.valueserie(0) - - if (abs(D.valueserie(0) - 16.5) > 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 - # - # float simple - # ------------ - x = 1. - print type(x) - xb = 2. - yo = 5. - H = numpy.matrix(numpy.identity(1)) - Hx = numpy.dot(H,x) - Hx = Hx.T - B = 1. - R = 1. - # - D.calculate( x = x, Hx = Hx, xb = xb, 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(1) - if (abs(D.valueserie(1) - 8.5) > 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 - diff --git a/src/daComposant/daDiagnostics/ComputeVariance.py b/src/daComposant/daDiagnostics/ComputeVariance.py deleted file mode 100644 index fefe24c..0000000 --- a/src/daComposant/daDiagnostics/ComputeVariance.py +++ /dev/null @@ -1,86 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Calcul de la variance d'un vecteur à chaque pas. Ce diagnostic très simple - est présent pour rappeller à l'utilisateur de l'assimilation qu'il faut - qu'il vérifie les variances de ses écarts en particulier. -""" -__author__ = "Jean-Philippe ARGAUD - Septembre 2008" - -import numpy -from daCore import BasicObjects, Persistence - -# ============================================================================== -class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar): - def __init__(self, name = "", unit = "", basetype = None, parameters = {}): - BasicObjects.Diagnostic.__init__(self, name, parameters) - Persistence.OneScalar.__init__( self, name, unit, basetype = float) - - def _formula(self, V): - """ - Calcul de la variance du vecteur en argument. Elle est faite avec une - division par la taille du vecteur. - """ - variance = V.var() - # - return variance - - def calculate(self, vector = None, step = None): - """ - Teste les arguments, active la formule de calcul et stocke le résultat - """ - if vector is None: - raise ValueError("One vector must be given to compute biais") - V = numpy.array(vector) - if V.size < 1: - raise ValueError("The given vector must not be empty") - # - value = self._formula( V) - # - self.store( value = value, step = step ) - -#=============================================================================== -if __name__ == "__main__": - print '\n AUTODIAGNOSTIC \n' - # - D = ElementaryDiagnostic("Ma variance") - # - # Vecteur de type matrix - # ---------------------- - x = numpy.matrix(([3., 4., 5.])) - print " Le vecteur de type 'matrix' choisi est..:", x - print " Le moyenne de ce vecteur est............:", x.mean() - print " La variance attendue de ce vecteur est..:", x.var() - # - D.calculate( vector = x) - print " La variance obtenue de ce vecteur est...:", D.valueserie(0) - print - # - # Vecteur de type array - # --------------------- - x = numpy.array(range(11)) - print " Le vecteur de type 'array' choisi est...:", x - print " Le moyenne de ce vecteur est............:", x.mean() - print " La variance attendue de ce vecteur est..:", x.var() - # - D.calculate( vector = x) - print " La variance obtenue de ce vecteur est...:", D.valueserie(1) - print diff --git a/src/daComposant/daDiagnostics/GaussianAdequation.py b/src/daComposant/daDiagnostics/GaussianAdequation.py deleted file mode 100644 index 660ed82..0000000 --- a/src/daComposant/daDiagnostics/GaussianAdequation.py +++ /dev/null @@ -1,154 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Diagnostic qui effectue le test du Khi2 pour juger de l'adéquation entre - la distribution d'un échantillon et une distribution gaussienne dont la - moyenne et l'écart-type sont calculés sur l'échantillon. - En input : la tolerance(tolerance) et le nombre de classes(nbclasse) - En output : Le resultat du diagnostic est une reponse booleenne au test : - True si l adequation a une distribution gaussienne est valide - au sens du test du Khi2, - False dans le cas contraire. -""" -__author__ = "Sophie RICCI - Juillet 2008" - -import numpy -from daCore import BasicObjects, Persistence -from ComputeKhi2 import ComputeKhi2_Gauss -import logging - -# ============================================================================== -class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar): - """ - """ - def __init__(self, name="", unit="", basetype = None, 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): - raise ValueError("A parameter named \"%s\" is required."%key) - - def formula(self, V): - """ - Effectue le calcul de la p-value pour un vecteur et une distribution - gaussienne et un nombre de classes donne en parametre du diagnostic. - """ - - [vectclasse, eftho, efobs, valeurKhi2, areaKhi2, message] = ComputeKhi2_Gauss( - vectorV = V, - dx = self.parameters["dxclasse"], - 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) ) - logging.info("valeurKhi2=%s"%valeurKhi2) - logging.info("areaKhi2=%s"%areaKhi2) - logging.info("tolerance=%s"%self.parameters["tolerance"]) - - if (areaKhi2 < (100.*self.parameters["tolerance"])) : - answerKhisquareTest = False - else: - answerKhisquareTest = True - logging.info( "La réponse au test est donc est %s"%answerKhisquareTest ) - return answerKhisquareTest - - def calculate(self, vector = None, step = None): - """ - Active la formule de calcul - """ - if vector is None: - raise ValueError("One vector must be given to calculate the Khi2 test") - V = numpy.array(vector) - if V.size < 1: - raise ValueError("The given vector must not be empty") - # - value = self.formula( V ) - # - self.store( value = value, step = step) - -# ============================================================================== -if __name__ == "__main__": - print "\n AUTODIAGNOSTIC \n" - - print " Test d adequation du khi-2 a une gaussienne pour un vecteur x" - print " connu de taille 1000, issu d'une distribution gaussienne normale" - print " en fixant la largeur des classes" - print - # - # Initialisation des inputs et appel du diagnostic - # ------------------------------------------------ - tolerance = 0.05 - dxclasse = 0.1 - D = ElementaryDiagnostic("AdequationGaussKhi2", parameters = { - "tolerance":tolerance, - "dxclasse":dxclasse, - "nbclasses":None, - }) - # - # Tirage de l'echantillon aleatoire - # --------------------------------- - numpy.random.seed(2490) - x = numpy.random.normal(50.,1.5,1000) - # - # Calcul - # ------ - D.calculate(x) - # - if D.valueserie(0) : - print " L'adequation a une distribution gaussienne est valide." - print - else : - raise ValueError("L'adéquation a une distribution gaussienne n'est pas valide.") - - - print " Test d adequation du khi-2 a une gaussienne pour u:n vecteur x" - print " connu de taille 1000, issu d'une distribution gaussienne normale" - print " en fixant le nombre de classes" - print - # - # Initialisation des inputs et appel du diagnostic - # ------------------------------------------------ - tolerance = 0.05 - nbclasses = 70. - D = ElementaryDiagnostic("AdequationGaussKhi2", parameters = { - "tolerance":tolerance, - "dxclasse":None, - "nbclasses":nbclasses - }) - # - # Tirage de l'echantillon aleatoire - # --------------------------------- - numpy.random.seed(2490) - x = numpy.random.normal(50.,1.5,1000) - # - # Calcul - # ------ - D.calculate(x) - # - if D.valueserie(0) : - print " L'adequation a une distribution gaussienne est valide." - print - else : - raise ValueError("L'adequation a une distribution gaussienne n'est pas valide.") - - diff --git a/src/daComposant/daDiagnostics/HLinearity.py b/src/daComposant/daDiagnostics/HLinearity.py deleted file mode 100644 index 0b3fe21..0000000 --- a/src/daComposant/daDiagnostics/HLinearity.py +++ /dev/null @@ -1,138 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Diagnotic de test sur la validité de l'hypothèse de linéarité de l'opérateur - H entre xp et xm - - Pour calculer Hlin on utilise un schéma différences finies centrées 2 - points. On définit un dxparam tel que : - xp = xb + dxparam - et - xm = xb - dxparam - On calcule Hxp et Hxm pour obtenir Hlin. Hlin est utilise dans le Blue pour - caler un paramêtre. La question importante est de choisir un dxparam pas - trop grand. - - On veut vérifier ici que l'hypothèse de linéarite du modèle par rapport au - paramêtre est valide sur l'intervalle du paramêtre [xm, xp]. Pour cela on - s'assure que l'on peut retrouver la valeur Hxb par les développemenents de - Taylor en xp et xm. Ainsi on calcule 2 estimations de Hxb, l'une à partir de - Hxp (notee Hx1) et l'autre à partir de Hxm (notee Hx2), que l'on compare à - la valeur calculée de Hxb. On s'intèresse ensuite a la distance entre Hxb et - ses estimés Hx1 et Hx2. Si la distance est inférieure a un seuil de - tolerance, l hypothese est valide. -""" -__author__ = "Sophie RICCI - Septembre 2008" - -import numpy -from daCore import AssimilationStudy, BasicObjects, Persistence - -# ============================================================================== -class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar): - def __init__(self, name="", unit="", basetype = None, 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.") - - def formula(self, H, dxparam, Hxp, Hxm, Hx): - """ - Test sur la validite de l hypothese de linearite de H entre xp et xm - """ - dimension = numpy.size(Hx) - # - # Reconstruit les valeurs Hx1 et Hx2 de Hx a partir de Hxm et Hxp - # --------------------------------------------------------------- - Hx1 = Hxm + H.T * dxparam - Hx2 = Hxp - H.T * dxparam - # - # Calcul de l'ecart entre Hx1 et Hx et entre Hx2 et Hx - # ---------------------------------------------------- - 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") - RMS.calculate(Hx1,Hx) - std1 = RMS.valueserie(0) - RMS.calculate(Hx2,Hx) - std2 = RMS.valueserie(1) - # - # Normalisation des écarts par Hx pour comparer a un pourcentage - # -------------------------------------------------------------- - RMS.calculate(Hx,Hx-Hx) - std = RMS.valueserie(2) - err1=std1/std - err2=std2/std - # - # Comparaison - # ----------- - if ( (err1 < self.parameters["tolerance"]) and (err2 < self.parameters["tolerance"]) ): - reponse = True - else: - reponse = False - return reponse - - def calculate(self, Hlin = None, deltaparam = None, Hxp = None, Hxm = None, Hx = None, step = None): - """ - Arguments : - - Hlin : Operateur d obsevation lineaire - - deltaparam : pas sur le parametre param - - Hxp : calcul en xp = xb + deltaparam - - Hxm : calcul en xm = xb - deltaparam - - Hx : calcul en x (generalement xb) - """ - value = self.formula( Hlin, deltaparam, Hxp, Hxm, Hx ) - # - self.store( value = value, step = step) - -#=============================================================================== -if __name__ == "__main__": - print '\n AUTODIAGNOSTIC \n' - - print " Diagnotic de test sur la validité de l'hypothèse de linéarité de" - print " l'opérateur H entre xp et xm." - print - # - dimension = 3 - # - # Définition des données - # ---------------------- - Hx = numpy.array(([ 2., 4., 6.])) - Hxp = numpy.array(([ 3., 5., 7.])) - Hxm = numpy.array(([ 1., 3., 5.])) - H = (Hxp - Hxm)/(2.) - dxparam = 1. - # - # Instanciation de l'objet diagnostic - # ----------------------------------- - D = ElementaryDiagnostic("Mon TestHlin", parameters = {"tolerance": 0.1}) - # - # Calcul - # ------ - D.calculate( Hlin = H, deltaparam = dxparam, Hxp = Hxp, Hxm = Hxm, Hx = Hx) - - # Validation du calcul - # -------------------- - if not D.valueserie(0) : - raise ValueError("La linearisation de H autour de x entre xm et xp est fausse pour ce cas test lineaire") - else : - print " La linéarisation de H autour de x entre xm et xp est valide pour ce cas-test linéaire." - print diff --git a/src/daComposant/daDiagnostics/HomogeneiteKhi2.py b/src/daComposant/daDiagnostics/HomogeneiteKhi2.py deleted file mode 100644 index fa7098d..0000000 --- a/src/daComposant/daDiagnostics/HomogeneiteKhi2.py +++ /dev/null @@ -1,116 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Diagnostic qui effectue le test du Khi2 pour juger de l'homogénéite entre - les distributions de 2 vecteurs quelconques. - - entrée : la tolerance (tolerance) et le nombre de classes (nbclasse), - sous forme de paramètres dans le dictionnaire Par - - sortie : le resultat du diagnostic est une reponse booleenne au test : - True si l homogeneite est valide au sens du test du Khi2, - False dans le cas contraire. -""" -__author__ = "Sophie RICCI - Juillet 2008" - -import numpy -from daCore import BasicObjects, Persistence -from ComputeKhi2 import ComputeKhi2_Homogen -import logging - -# ============================================================================== -class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar): - def __init__(self, name="", unit="", basetype = None, 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): - raise ValueError("A parameter named \"%s\" is required."%key) - - def _formula(self, V1, V2): - """ - Effectue le calcul de la p-value pour deux vecteurs et un nombre de - classes donne en parametre du diagnostic. - """ - [classes, eftheo, efobs, valeurKhi2, areaKhi2, message] = ComputeKhi2_Homogen( - vectorV1 = V1, - vectorV2 = V2, - dx = self.parameters["dxclasse"], - nbclasses = self.parameters["nbclasses"], - SuppressEmptyClasses = True) - # - logging.info( message ) - logging.info( "(si <%.2f %s on refuse effectivement l'homogeneite)"%(100.*self.parameters["tolerance"],"%") ) - # - answerKhisquareTest = False - if (areaKhi2 < (100.*self.parameters["tolerance"])) : - answerKhisquareTest = False - else: - answerKhisquareTest = True - # - 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 Khi2 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") - # - value = self._formula( V1, V2 ) - # - self.store( value = value, step = step ) - -# ============================================================================== -if __name__ == "__main__": - print "\n AUTODIAGNOSTIC \n" - - print " Test d'homogeneite du Khi-2 pour deux vecteurs de taille 10," - print " issus d'une distribution gaussienne normale" - print - # - # Initialisation des inputs et appel du diagnostic - # -------------------------------------------------------------------- - tolerance = 0.05 - dxclasse = 0.5 - D = ElementaryDiagnostic("HomogeneiteKhi2", parameters = { - "tolerance":tolerance, - "dxclasse":dxclasse, - "nbclasses":None, - }) - # - # Tirage de l'echantillon aleatoire - # -------------------------------------------------------------------- - numpy.random.seed(4000) - x1 = numpy.random.normal(50.,1.5,10000) - numpy.random.seed(2490) - x2 = numpy.random.normal(50.,1.5,10000) - # - # Calcul - # -------------------------------------------------------------------- - D.calculate(x1, x2) - # - print " La reponse du test est \"%s\" pour une tolerance de %.2e et une largeur de classe de %.2e "%(D.valueserie(0), tolerance, dxclasse) - print diff --git a/src/daComposant/daDiagnostics/ReduceBiais.py b/src/daComposant/daDiagnostics/ReduceBiais.py deleted file mode 100644 index f29610d..0000000 --- a/src/daComposant/daDiagnostics/ReduceBiais.py +++ /dev/null @@ -1,107 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Diagnostic sur la reduction du biais lors de l'analyse -""" -__author__ = "Sophie RICCI - Aout 2008" - -import numpy -from daCore import BasicObjects, Persistence - -# ============================================================================== -class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar): - def __init__(self, name = "", unit = "", basetype = None, parameters = {}): - BasicObjects.Diagnostic.__init__(self, name, parameters) - Persistence.OneScalar.__init__( self, name, unit, basetype = bool) - - def _formula(self, V1, V2): - """ - Vérification de la reduction du biais entre OMB et OMA lors de l'analyse - """ - biaisOMB = V1.mean() - biaisOMA = V2.mean() - # - if biaisOMA > biaisOMB: - reducebiais = False - else : - reducebiais = True - # - return reducebiais - - def calculate(self, vectorOMB = None, vectorOMA = None, step = None): - """ - Teste les arguments, active la formule de calcul et stocke le résultat - Arguments : - - vectorOMB : vecteur d'écart entre les observations et l'ébauche - - vectorOMA : vecteur d'écart entre les observations et l'analyse - """ - if ( (vectorOMB is None) or (vectorOMA is None) ): - raise ValueError("Two vectors must be given to test the reduction of the biais after analysis") - V1 = numpy.array(vectorOMB) - V2 = numpy.array(vectorOMA) - 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") - # - value = self._formula( V1, V2 ) - # - self.store( value = value, step = step ) - -#=============================================================================== -if __name__ == "__main__": - print '\n AUTODIAGNOSTIC \n' - # - # Instanciation de l'objet diagnostic - # ----------------------------------- - D = ElementaryDiagnostic("Mon ReduceBiais") - # - # Tirage des 2 vecteurs choisis - # ------------------------------- - x1 = numpy.matrix(([3. , 4., 5. ])) - x2 = numpy.matrix(([1.5, 2., 2.5])) - print " L'écart entre les observations et l'ébauche est OMB :", x1 - print " La moyenne de OMB (i.e. le biais) est de............:", x1.mean() - print " L'écart entre les observations et l'analyse est OMA :", x2 - print " La moyenne de OMA (i.e. le biais) est de............:", x2.mean() - # - D.calculate( vectorOMB = x1, vectorOMA = x2) - if not D.valueserie(0) : - print " Résultat : l'analyse NE RÉDUIT PAS le biais" - else : - print " Résultat : l'analyse RÉDUIT le biais" - print - # - # Tirage des 2 vecteurs choisis - # ------------------------------- - x1 = numpy.matrix(range(-5,6)) - x2 = numpy.array(range(11)) - print " L'écart entre les observations et l'ébauche est OMB :", x1 - print " La moyenne de OMB (i.e. le biais) est de............:", x1.mean() - print " L'écart entre les observations et l'analyse est OMA :", x2 - print " La moyenne de OMA (i.e. le biais) est de............:", x2.mean() - # - D.calculate( vectorOMB = x1, vectorOMA = x2) - if not D.valueserie(1) : - print " Résultat : l'analyse NE RÉDUIT PAS le biais" - else : - print " Résultat : l'analyse RÉDUIT le biais" - print diff --git a/src/daComposant/daDiagnostics/VarianceOrder.py b/src/daComposant/daDiagnostics/VarianceOrder.py deleted file mode 100644 index 49cbf5b..0000000 --- a/src/daComposant/daDiagnostics/VarianceOrder.py +++ /dev/null @@ -1,126 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Diagnostic sur les variances dans B et R par rapport à l'ébauche Xb et aux - observations Y. On teste si on a les conditions : - 1%*xb < sigma_b < 10%*xb - et - 1%*yo < sigma_o < 10%*yo - Le diagnostic renvoie True si les deux conditions sont simultanément - vérifiées, False dans les autres cas. -""" -__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Septembre 2008" - -import numpy -from scipy.linalg import eig -from daCore import BasicObjects, Persistence -import logging - -# ============================================================================== -class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar): - def __init__(self, name = "", unit = "", basetype = None, parameters = {}): - BasicObjects.Diagnostic.__init__(self, name, parameters) - Persistence.OneScalar.__init__( self, name, unit, basetype = bool ) - - def _formula(self, xb, B, yo, R): - """ - Comparaison des variables et de leur variance relative - """ - valpB = eig(B, left = False, right = False) - valpR = eig(R, left = False, right = False) - logging.info(" Si l on souhaite 1%s*xb < sigma_b < 10%s*xb, les valeurs propres de B doivent etre comprises dans l intervalle [%.3e,%.3e]"%("%","%",1.e-4*xb.mean()*xb.mean(),1.e-2*xb.mean()*xb.mean())) - logging.info(" Si l on souhaite 1%s*yo < sigma_o < 10%s*yo, les valeurs propres de R doivent etre comprises dans l intervalle [%.3e,%.3e]"%("%","%",1.e-4*yo.mean()*yo.mean(),1.e-2*yo.mean()*yo.mean())) - # - limite_inf_valp = 1.e-4*xb.mean()*xb.mean() - limite_sup_valp = 1.e-2*xb.mean()*xb.mean() - variancexb = (valpB >= limite_inf_valp).all() and (valpB <= limite_sup_valp).all() - logging.info(" La condition empirique sur la variance de Xb est....: %s"%variancexb) - # - limite_inf_valp = 1.e-4*yo.mean()*yo.mean() - limite_sup_valp = 1.e-2*yo.mean()*yo.mean() - varianceyo = (valpR >= limite_inf_valp).all() and (valpR <= limite_sup_valp).all() - logging.info(" La condition empirique sur la variance de Y est.....: %s",varianceyo) - # - variance = variancexb and varianceyo - logging.info(" La condition empirique sur la variance globale est..: %s"%variance) - # - return variance - - def calculate(self, Xb = None, B = None, Y = None, R = None, step = None): - """ - Teste les arguments, active la formule de calcul et stocke le résultat - Arguments : - - Xb : valeur d'ébauche du paramêtre - - B : matrice de covariances d'erreur d'ébauche - - yo : vecteur d'observation - - R : matrice de covariances d'erreur d'observation - """ - if (Xb is None) or (B is None) or (Y is None) or (R is None): - raise ValueError("You must specify Xb, B, Y, R") - yo = numpy.array(Y) - BB = numpy.matrix(B) - xb = numpy.array(Xb) - RR = numpy.matrix(R) - if (RR.size < 1 ) or (BB.size < 1) : - raise ValueError("The background and the observation covariance matrices must not be empty") - if ( yo.size < 1 ) or ( xb.size < 1 ): - raise ValueError("The Xb background and the Y observation vectors must not be empty") - if xb.size*xb.size != BB.size: - raise ValueError("Xb background vector and B covariance matrix sizes are not consistent") - if yo.size*yo.size != RR.size: - raise ValueError("Y observation vector and R covariance matrix sizes are not consistent") - if yo.all() == 0. or xb.all() == 0. : - raise ValueError("The diagnostic can not be applied to zero vectors") - # - value = self._formula( xb, BB, yo, RR) - # - self.store( value = value, step = step ) - -#=============================================================================== -if __name__ == "__main__": - print '\n AUTODIAGNOSTIC \n' - # - # Instanciation de l'objet diagnostic - # ----------------------------------- - D = ElementaryDiagnostic("Mon OrdreVariance") - # - # Vecteur de type matrix - # ---------------------- - xb = numpy.array([11000.]) - yo = numpy.array([1.e12 , 2.e12, 3.e12 ]) - B = 1.e06 * numpy.matrix(numpy.identity(1)) - R = 1.e22 * numpy.matrix(numpy.identity(3)) - # - D.calculate( Xb = xb, B = B, Y = yo, R = R) - print " L'ébauche est.......................................:",xb - print " Les observations sont...............................:",yo - print " La valeur moyenne des observations est..............: %.2e"%yo.mean() - print " La valeur moyenne de l'ebauche est..................: %.2e"%xb.mean() - print " La variance d'ébauche specifiée est.................: %.2e"%1.e6 - print " La variance d'observation spécifiée est.............: %.2e"%1.e22 - # - if D.valueserie(0) : - print " Les variances specifiées sont de l'ordre de 1% a 10% de l'ébauche et des observations" - else : - print " Les variances specifiées ne sont pas de l'ordre de 1% a 10% de l'ébauche et des observations" - print - - diff --git a/src/daComposant/daNumerics/ComputeKhi2.py b/src/daComposant/daNumerics/ComputeKhi2.py deleted file mode 100644 index d5d8846..0000000 --- a/src/daComposant/daNumerics/ComputeKhi2.py +++ /dev/null @@ -1,492 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Outil numerique de calcul de la variable Khi2 - - On peut realiser deux types de test du Khi2 : - - test d'adequation : comparer la distribution d'un echantillon a une - distribution theorique, - - test d'homogeneite : comparer les distributions de 2 vecteurs. - - Pour le test d'adequation, on travaille sur une gaussienne - dont la moyenne et l'ecart type sont calcules sur - l'echantillon, soit donnes. - - Ce fichier contient une classe "StatspourTests" de methodes qui realisent - differentes etapes utiles aux calculs des tests du Khi2. - - Ce fichier contient de plus 3 methodes : ComputeKhi2_testGauss, - ComputeKhi2_Gauss et ComputeKhi2_Homogen. - - ComputeKhi2_testGauss : calcul la distance du Khi2 entre un vecteur - aleatoire issu d un gaussienne et une distribution theorique gaussienne - dont on specifie la moyenne et l ecart type - - ComputeKhi2_Gauss : calcul la distance du Khi2 entre un vecteur donne et - une distribution theorique gaussienne dont la moyenne et l ecart type sont - calcules sur l echantillon - - ComputeKhi2_Homogen : calcul la distance du Khi2 entre deux vecteurs donnes - - Ces methodes necessitent et fournissent : - - en input : - - le ou les vecteurs dont on etudie la distribution, - - la distribution theorique et eventuellement la moyenne et ecart type, - - la largeur des classes, - - un booleen traduisant la suppression des classes vides - - en output : - - le vecteur des classes, - - les pdf theorique et donnee, - - la valeur du Khi2, - - la p-value qui represent l'aire de la queue de la distribution du - Khi2 et - - le message qui interprete le test. -""" -__author__ = "Sophie RICCI - Mars 2010" - -import numpy -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 - -# ============================================================================== -class StatspourTests : - """ - Classe de methodes pour la preparation du test de Khi2 - """ - def __init__(self, cdftheo=None, meantheo = None, stdtheo = None, pdftest=None,obs=None,use_mean_std_exp=True, dxmin=0.01, obsHomogen = None, nbclasses = None) : - - - if (pdftest is None and obs is None) : - raise ValueError('Donner soit une pdf de test soit un vecteur obs') - if not obs is None : - if pdftest is None : - self.__obs=asarray(obs) - if not pdftest is None : - if obs is None : - if len(pdftest) == 3: - niter=eval(pdftest[2]) - obs=[eval(" ".join(pdftest[:2])) for z in range(niter)] - self.__obs=asarray(obs) - else : - self.__obs=asarray(eval(" ".join(pdftest[:2]))) - if not (obsHomogen is None) : - self.__obsHomogen = asarray(obsHomogen) - self.__testHomogen = True - else : - self.__testHomogen = False - - - self.__mean_exp = self.__obs.mean() - self.__std_exp = self.__obs.std() - - if cdftheo is None : raiseValueError(" ... Definir le parametre cdftheo ...") - if use_mean_std_exp : - self.__cdf=cdftheo( self.__mean_exp, self.__std_exp).cdf - else : - self.__cdf=cdftheo( meantheo, stdtheo).cdf - - self.__min=min(self.__obs) - self.__max=max(self.__obs) - self.__N=len(self.__obs) - self.__use_mean_std_exp=use_mean_std_exp - self.__dxmin=dxmin - self.__nbclasses = nbclasses - if not (dxmin is None) and not (nbclasses is None) : - raise ValueError("... Specifier soit le nombre de classes, soit la largeur des classes") - if (dxmin is None) and (nbclasses is None) : - raise ValueError("... Specifier soit le nombre de classes, soit la largeur des classes") - if not (nbclasses is None) and (dxmin is None) : - self.__dxmin = (self.__max - self.__min ) / float(self.__nbclasses) - return None - - def MakeClasses(self) : - """ - Classification en classes - """ - self.__subdiv=arange(self.__min,self.__max+self.__dxmin,self.__dxmin) - self.__modalites=len(self.__subdiv) - return None - - def ComputeObs(self): - """ - Calcul de la probabilite observee de chaque classe - """ - self.__kobs=histogram2(self.__obs,self.__subdiv)[1:] - return self.__kobs - - def ComputeObsHomogen(self): - """ - Calcul de la probabilite observee pour le test homogeneite de chaque classe - """ - self.__kobsHomogen=histogram2(self.__obsHomogen,self.__subdiv)[1:] - return self.__kobsHomogen - - def ComputeTheo(self): - """ - Calcul de la probabilite theorique de chaque classe - """ - self.__ktheo=[self.__cdf(self.__subdiv[i+1])-self.__cdf(self.__subdiv[i]) for i in range(self.__modalites-1)] - self.__ktheo=asarray(self.__ktheo) - self.__ktheo=(sum(self.__kobs)/sum(self.__ktheo))*self.__ktheo - - def Computepdfs(self) : - - self.__subdiv=self.__subdiv[1:] - self.__pdfobs=[self.__kobs[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)] - - if self.__testHomogen : - self.__pdftheo=[self.__kobsHomogen[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)] - else : - self.__pdftheo=[self.__ktheo[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)] - - return self.__subdiv, self.__pdftheo, self.__pdfobs - - def Computeddl(self): - """ - Calcul du nombre de degres de liberte - """ - self.__ddl = self.__modalites - 1. - if self.__use_mean_std_exp : - self.__ddl = self.__ddl - 2. - if (self.__ddl < 1.): - raise ValueError("The ddl is 0, you must increase the number of classes nbclasse ") - logging.debug("Nombre de degres de liberte=%s"%self.__ddl) - - def ComputeValue(self) : - """ - Calcul de la variable Q qui suit une loi Khi-2 - """ - if self.__testHomogen : - kobs,ktheo=self.__kobs.tolist(),self.__kobsHomogen.tolist() - else : - kobs,ktheo=self.__kobs.tolist(),self.__ktheo.tolist() - - # on supprime les classes theoriques qui ont moins d'un element (sinon la distance khi2 tendrait vers l'infini) - ko,kt=[],[] - self.__count0 = 0. - for k,val in enumerate(ktheo): - if val > 1.0: - kt.append(val) - ko.append(kobs[k]) - else : - self.__count0 = self.__count0 +1. - logging.debug("WARNING : nombre de classes vides supprimees (effectif theorique inferieur a 1.) pour le calcul de la valeur du Khi2 = %s"%self.__count0) - ef1,ef2=asarray(ko),asarray(kt) - count = 0. - for el in ef1.tolist() : - if el < 5. : - count = count +1. - for el in ef2.tolist() : - if el < 5. : - count = count +1. - pourcent_nbclasse_effecinf = count /(2.*len(ef1.tolist())) *100. - if (pourcent_nbclasse_effecinf > 20.) : - logging.debug("WARNING : nombre de classes dont l effectif est inferieur a 5 elements %s"%pourcent_nbclasse_effecinf) - k,p = chisquare(ef1, ef2) - k2, p2 = [k],[p] - for shift in range(1,6): - k,p=chisquare(ef1[shift:],ef2[:-shift]) - k2.append(k) - p2.append(p) - k,p=chisquare(ef1[:-shift],ef2[shift:]) - k2.append(k) - p2.append(p) - logging.debug("Liste des valeurs du Khi2 = %s"%k2) - self.__khi2=min(k2) - self.__Q=self.__khi2 - - 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): - """ - Interpretation du test - """ - message = "Il y a %.2f%s de chance de se tromper en refusant l'homogeneite"%(self.__areakhi2,"%") - return message - -# ============================================================================== -def ComputeKhi2_testGauss( - meantheo = 0., - stdtheo = 1., - nech = 10, - dx = 0.1, - nbclasses = None, - SuppressEmptyClasses = True, - ): - """ - Test du Khi2 d adequation entre tirage aleatoire dans gaussienne et une gaussienne theo - """ - essai = StatspourTests( cdftheo=norm, meantheo = meantheo, stdtheo = stdtheo, pdftest = ["random.normal","(%.3f,%.2f,%d)"%(meantheo,stdtheo,nech)], obs = None, use_mean_std_exp=False,dxmin=dx, obsHomogen = None, nbclasses = nbclasses) - essai.MakeClasses() - essai.ComputeObs() - essai.ComputeTheo() - classes,eftheo, efobs = essai.Computepdfs() - essai.Computeddl() - valeurKhi2= essai.ComputeValue() - areaKhi2 = essai.ComputeArea() - message = essai.WriteMessage() - logging.debug("message %s"%message) - return classes, eftheo, efobs, valeurKhi2, areaKhi2, message - -def ComputeKhi2_Gauss( - vectorV = None, - dx = 0.1, - SuppressEmptyClasses = True, - nbclasses = None - ): - """ - Test du Khi2 d 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() - valeurKhi2= essai.ComputeValue() - areaKhi2 = essai.ComputeArea() - message = essai.WriteMessage() - 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, - dx = 0.1, - SuppressEmptyClasses = True, - nbclasses = None - ): - """ - Test du Khi2 d homogeniete entre 2 vecteurs - """ - essai = StatspourTests( cdftheo=norm, pdftest = None, obs = vectorV1, use_mean_std_exp=True,dxmin=dx, obsHomogen = vectorV2, nbclasses = nbclasses) - essai.MakeClasses() - essai.ComputeObs() - essai.ComputeObsHomogen() - classes,eftheo, efobs = essai.Computepdfs() - essai.Computeddl() - valeurKhi2= essai.ComputeValue() - areaKhi2 = essai.ComputeArea() - message = essai.WriteMessageHomogen() - logging.debug("message %s"%message) - return classes, eftheo, efobs, valeurKhi2, areaKhi2, message - -# ============================================================================== -if __name__ == "__main__": - 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' - classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_testGauss(meantheo = 0., stdtheo = 1., nech = 1000., dx = 0.1, SuppressEmptyClasses = True, nbclasses = None) - print ' valeurKhi2=',valeurKhi2 - print ' areaKhi2=',areaKhi2 - print ' ',message - - if (numpy.abs(areaKhi2 - 99.91)< 1.e-2) : - print "The computation of the khisquare value is OK" - else : - 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' - V = random.normal(50.,1.5,1000) - classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = 0.1, vectorV = V, SuppressEmptyClasses = True, nbclasses = None) - print ' valeurKhi2=',valeurKhi2 - print ' areaKhi2=',areaKhi2 - print ' ',message - - if (numpy.abs(areaKhi2 - 99.60)< 1.e-2) : - 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' - V1 = random.normal(50.,1.5,10000) - numpy.random.seed(2490) - V2 = random.normal(50.,1.5,10000) - classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Homogen(dx = 0.5, vectorV1 = V1, vectorV2 = V2, SuppressEmptyClasses = True, nbclasses = None) - print ' valeurKhi2=',valeurKhi2 - print ' areaKhi2=',areaKhi2 - print ' ',message - - if (numpy.abs(areaKhi2 - 99.98)< 1.e-2) : - 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' -# file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech10000.gnu' -# fid = open(file, "w") -# lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' ) - numpy.random.seed(4000) - V = random.normal(0., 1.,10000) - aire = [] - for dx in arange(0.01, 1., 0.001) : - classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None) -# lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2) - aire.append(areaKhi2) - meanaire = numpy.asarray(aire) -# fid.writelines(lines) - - print " En moyenne, il y a ", meanaire.mean(),"% de chance de se tromper en refusant l adequation a une loi gaussienne pour un echantillon de taille 10000" - print - if (numpy.abs( meanaire.mean() - 71.79)< 1.e-2) : - 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' -# file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech1000.gnu' -# fid = open(file, "w") -# lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' ) - numpy.random.seed(4000) - V = random.normal(0., 1.,1000) - aire = [] - for dx in arange(0.05, 1., 0.001) : - classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None) -# lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2) - aire.append(areaKhi2) - meanaire = numpy.asarray(aire) -# fid.writelines(lines) - - print " En moyenne, il y a ", meanaire.mean(),"% de chance de se tromper en refusant l adequation a une loi gaussienne pour un echantillon de taille 1000" - print - if (numpy.abs( meanaire.mean() - 90.60)< 1.e-2) : - 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 '' - print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 100' -# file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech100.gnu' -# fid = open(file, "w") -# lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' ) - numpy.random.seed(4000) - V = random.normal(0., 1.,100) - aire = [] - for dx in arange(0.1, 1., 0.01) : - classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None) -# lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2) - aire.append(areaKhi2) - meanaire = numpy.asarray(aire) -# fid.writelines(lines) - - print " En moyenne, il y a ", meanaire.mean(),"% de chance de se tromper en refusant l adequation a une loi gaussienne pour un echantillon de taille 100" - print diff --git a/src/daComposant/daNumerics/ComputeStudent.py b/src/daComposant/daNumerics/ComputeStudent.py deleted file mode 100644 index 9d2fd10..0000000 --- a/src/daComposant/daNumerics/ComputeStudent.py +++ /dev/null @@ -1,260 +0,0 @@ -#-*-coding:iso-8859-1-*- -# -# 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 -# -__doc__ = """ - Outil numérique de calcul des variables de Student pour 2 vecteurs - dépendants ou indépendants, avec variances supposées égales ou différentes -""" -__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Octobre 2008" - -import sys ; sys.path.insert(0, "../daCore") - -import numpy -from scipy.stats import ttest_rel, ttest_ind, betai -import logging - -# ============================================================================== -def DependantVectors(vector1 = None, vector2 = None, tolerance = 0.05 ): - """ - Outil numérique de calcul de la variable de Student pour 2 vecteurs - dépendants - 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 reponse au test pour une tolerance 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 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") - # - # Calcul de la p-value du Test de Student - # -------------------------------------------------------------------- - [t, prob] = ttest_rel(V1, V2) - areastudent = 100. * prob - # - logging.debug("DEPENDANTVECTORS t = %.3f, areastudent = %.3f"%(t, areastudent)) - # - # Tests - # -------------------------------------------------------------------- - message = "DEPENDANTVECTORS Il y a %.2f %s de chance de se tromper en refusant l'hypothèse d'égalité des moyennes des 2 échantillons dépendants (si <%.2f %s on refuse effectivement l'égalité)"%(areastudent, "%", 100.*tolerance,"%") - logging.debug(message) - # - if (areastudent < (100.*tolerance)) : - answerTestStudent = False - else: - answerTestStudent = True - # - return areastudent, t, answerTestStudent, message - -# ============================================================================== -def IndependantVectorsDifferentVariance(vector1 = None, vector2 = None, tolerance = 0.05 ): - """ - Outil numerique de calcul de la variable de Student pour 2 vecteurs independants supposes de variances vraies differentes - En input : la tolerance - En output : la p-value, la valeur de la variable aleatoire, la reponse au test pour une tolerance 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 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") - # - # Calcul de la p-value du Test de Student - # -------------------------------------------------------------------- - # t = (m1 - m2)/ sqrt[ (var1/n1 + var2/n2) ] - # ou var est calcule comme var = somme (xi -xmena)**2 /(n-1) - n1 = V1.size - n2 = V2.size - mean1 = V1.mean() - mean2 = V2.mean() - var1 = numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std() * numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std() - var2 = numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std() * numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std() - t = (mean1 - mean2)/ numpy.sqrt( var1/n1 + var2/n2 ) - df = ( (var1/n1 + var2/n2) * (var1/n1 + var2/n2) ) / ( (var1/n1)*(var1/n1)/(n1-1) + (var2/n2)*(var2/n2)/(n2-1) ) - zerodivproblem = var1/n1 + var2/n2 == 0 - t = numpy.where(zerodivproblem, 1.0, t) # replace NaN t-values with 1.0 - prob = betai(0.5*df,0.5,float(df)/(df+t*t)) - areastudent = 100. * prob - # - logging.debug("IndependantVectorsDifferentVariance t = %.3f, areastudent = %.3f"%(t, areastudent)) - # - # Tests - # -------------------------------------------------------------------- - message = "IndependantVectorsDifferentVariance Il y a %.2f %s de chance de se tromper en refusant l'hypothèse d'égalité des moyennes des 2 échantillons indépendants supposés de variances différentes (si <%.2f %s on refuse effectivement l'égalité)"%(areastudent, "%", 100.* tolerance,"%") - logging.debug(message) - if (areastudent < (100.*tolerance)) : - answerTestStudent = False - else: - answerTestStudent = True - # - return areastudent, t, answerTestStudent, message - -# ============================================================================== -def IndependantVectorsEqualVariance(vector1 = None, vector2 = None, tolerance = 0.05 ): - """ - Outil numerique de calcul de la variable de Student pour 2 vecteurs independants supposes de meme variance vraie - En input : la tolerance - En output : la p-value, la valeur de la variable aleatoire, la reponse au test pour une tolerance 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 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") - # - # Calcul de la p-value du Test de Student - # -------------------------------------------------------------------- - # t = sqrt(n1+n2-2) * (m1 - m2)/ sqrt[ (1/n1 +1/n2) * ( (n1-1)var1 + (n2-1)var2 )] - # ou var est calcule comme var = somme (xi -xmena)**2 /(n-1) - [t, prob] = ttest_ind(V1, V2) - areastudent = 100. * prob - # - logging.debug("IndependantVectorsEqualVariance t = %.3f, areastudent = %.3f"%(t, areastudent)) - # Tests - # -------------------------------------------------------------------- - message = "IndependantVectorsEqualVariance Il y a %.2f %s de chance de se tromper en refusant l'hypothèse d'égalité des moyennes des 2 échantillons indépendants supposés de même variance (si <%.2f %s on refuse effectivement l'égalité)"%(areastudent, "%", 100.* tolerance,"%") - logging.debug(message) - if (areastudent < (100.*tolerance)) : - answerTestStudent = False - else: - answerTestStudent = True - - return areastudent, t, answerTestStudent, message - -# ============================================================================== -if __name__ == "__main__": - print '\n AUTODIAGNOSTIC \n' - # logging.getLogger().setLevel(logging.DEBUG) - - print - print " Test de Student pour des vecteurs dépendants" - print " --------------------------------------------" - # Tirage de l'echantillon - V1 = numpy.matrix(([-1., 0., 4.])).T - V2 = numpy.matrix(([-2., 0., 8.])).T - V1 = V1.A1 - V2 = V2.A1 - # - # Appel de l outil DependantVectors et initialisation des inputs - [aire, Q, reponse, message] = DependantVectors( - vector1 = V1, - vector2 = V2, - tolerance = 0.05) - # - # Verification par les calculs sans les routines de scipy.stats - # (ref numerical recipes) - n = V1.size - df= n -1 - # Les routines de scipy.stats utilisent une variance calculee avec n-1 et non n comme dans std - # t = (m1 - m2)/ sqrt[(varx1 + varx2 - 2 cov(x1, x2))/n ] - # ou var est calcule comme var = somme (xi -xmean)**2 /(n-1) - var1 = numpy.sqrt(n)/numpy.sqrt(n-1)* V1.std() * numpy.sqrt(n)/numpy.sqrt(n-1) * V1.std() - var2 = numpy.sqrt(n)/numpy.sqrt(n-1)* V2.std() * numpy.sqrt(n)/numpy.sqrt(n-1) * V2.std() - m1 = V1.mean() - m2 = V2.mean() - cov = 0. - for j in range(0, n) : - cov = cov + (V1[j] - m1)*(V2[j] - m2) - cov = cov /df - sd = numpy.sqrt((var1 + var2 - 2. *cov) / n) - tverif = (m1 -m2) /sd - aireverif = 100. * betai(0.5*df,0.5,float(df)/(df+tverif*tverif)) - if (aireverif - aire < 1.e-5) : - print " Le calcul est conforme à celui de l'algorithme du Numerical Recipes" - else : - raise ValueError("Le calcul n'est pas conforme à celui de l'algorithme Numerical Recipes") - - if (numpy.abs(aire - 57.99159)< 1.e-5) : - print " Le calcul est JUSTE sur cet exemple." - else : - raise ValueError("Le calcul est FAUX sur cet exemple.") - - print - print " Test de Student pour des vecteurs independants supposés de même variance" - print " ------------------------------------------------------------------------" - # Tirage de l'echantillon - V1 = numpy.matrix(([-1., 0., 4.])).T - V2 = numpy.matrix(([-2., 0., 8.])).T - V1 = V1.A1 - V2 = V2.A1 - # - # Appel de l outil IndependantVectorsDifferentVariance et initialisation des inputs - [aire, Q, reponse, message] = IndependantVectorsDifferentVariance( - vector1 = V1, - vector2 = V2, - tolerance = 0.05) - # - if (numpy.abs(aire - 78.91339)< 1.e-5) : - print " Le calcul est JUSTE sur cet exemple." - else : - raise ValueError("Le calcul est FAUX sur cet exemple.") - - print - print " Test de Student pour des vecteurs indépendants supposés de même variance" - print " ------------------------------------------------------------------------" - # Tirage de l'echantillon - V1 = numpy.matrix(([-1., 0., 4.])).T - V2 = numpy.matrix(([-2., 0., 8.])).T - V1 = V1.A1 - V2 = V2.A1 - # - # Appel de l outil IndependantVectorsEqualVariance et initialisation des inputs - [aire, Q, reponse, message] = IndependantVectorsEqualVariance( - vector1 = V1, - vector2 = V2, - tolerance = 0.05) - # - # Verification par les calculs sans les routines de scipy.stats (ref numerical recipes) - n1 = V1.size - n2 = V2.size - df= n1 + n2 -2 - # Les routines de scipy.stats utilisent une variance calculee avec n-1 et non n comme dans std - var1 = numpy.sqrt(n1)/numpy.sqrt(n1-1)* V1.std() * numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std() - var2 = numpy.sqrt(n2)/numpy.sqrt(n2-1)* V2.std() * numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std() - m1 = V1.mean() - m2 = V2.mean() - var = ((n1 -1.) *var1 + (n2 -1.) *var2 ) /df - tverif = (m1 -m2) /numpy.sqrt(var*(1./n1 + 1./n2)) - aireverif = 100. * betai(0.5*df,0.5,float(df)/(df+tverif*tverif)) - # - if (aireverif - aire < 1.e-5) : - print " Le calcul est conforme à celui de l'algorithme du Numerical Recipes" - else : - raise ValueError("Le calcul n'est pas conforme à celui de l'algorithme Numerical Recipes") - - if (numpy.abs(aire - 78.42572)< 1.e-5) : - print " Le calcul est JUSTE sur cet exemple." - else : - raise ValueError("Le calcul est FAUX sur cet exemple.") - - print -- 2.39.2