idl/Makefile
resources/Makefile
src/Makefile
-
+ src/daComposant/Makefile
])
AC_OUTPUT
-SUBDIRS=
+SUBDIRS= daComposant
--- /dev/null
+include $(top_srcdir)/adm_local/make_common_starter.am
+
+EXTRA_DIST = daCore
+
+DIR = $(top_srcdir)/src/daComposant/
+
+install-data-local:
+ ${mkinstalldirs} $(DESTDIR)$(salomepythondir)
+ cp -R $(DIR)daCore $(DESTDIR)$(salomepythondir)
+
+uninstall-local:
+ chmod -R +w $(DESTDIR)$(salomepythondir)/daCore
+ rm -rf $(DESTDIR)$(salomepythondir)/daCore
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Algorithme variationnel statique (3D-VAR)
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2009"
+
+import sys ; sys.path.insert(0, "../daCore")
+import logging
+import Persistence
+from BasicObjects import Algorithm
+import PlatformInfo ; m = PlatformInfo.SystemUsage()
+
+import numpy
+import scipy.optimize
+
+if logging.getLogger().level < 30:
+ iprint = 1
+ message = scipy.optimize.tnc.MSG_ALL
+ disp = 1
+else:
+ iprint = -1
+ message = scipy.optimize.tnc.MSG_NONE
+ disp = 0
+
+# ==============================================================================
+class ElementaryAlgorithm(Algorithm):
+ def __init__(self):
+ Algorithm.__init__(self)
+ 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):
+ """
+ Calcul de l'estimateur 3D-VAR
+ """
+ logging.debug("%s Lancement"%self._name)
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+ #
+ Hm = H["Direct"].appliedTo
+ Ht = H["Adjoint"].appliedInXTo
+ #
+ # Utilisation éventuelle d'un vecteur H(Xb) précalculé
+ # ----------------------------------------------------
+ if H["AppliedToX"] is not None and H["AppliedToX"].has_key("HXb"):
+ logging.debug("%s Utilisation de HXb"%self._name)
+ HXb = H["AppliedToX"]["HXb"]
+ else:
+ logging.debug("%s Calcul de Hm(Xb)"%self._name)
+ HXb = Hm( Xb )
+ #
+ # Calcul du préconditionnement
+ # ----------------------------
+ # Bdemi = numpy.linalg.cholesky(B)
+ #
+ # Calcul de l'innovation
+ # ----------------------
+ d = Y - HXb
+ logging.debug("%s Innovation d = %s"%(self._name, d))
+ #
+ # Précalcul des inversion appellée dans les fonction-coût et gradient
+ # -------------------------------------------------------------------
+ BI = B.I
+ RI = R.I
+ #
+ # Définition de la fonction-coût
+ # ------------------------------
+ def CostFunction(x):
+ _X = numpy.asmatrix(x).flatten().T
+ logging.info("%s CostFunction X = %s"%(self._name, numpy.asmatrix( _X ).flatten()))
+ _HX = Hm( _X )
+ _HX = numpy.asmatrix(_HX).flatten().T
+ Jb = 0.5 * (_X - Xb).T * BI * (_X - Xb)
+ Jo = 0.5 * (Y - _HX).T * RI * (Y - _HX)
+ J = float( Jb ) + float( Jo )
+ 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["CostFunctionJb"].store( Jb )
+ self.StoredVariables["CostFunctionJo"].store( Jo )
+ self.StoredVariables["CostFunctionJ" ].store( J )
+ return float( J )
+ #
+ def GradientOfCostFunction(x):
+ _X = numpy.asmatrix(x).flatten().T
+ logging.info("%s GradientOfCostFunction X = %s"%(self._name, numpy.asmatrix( _X ).flatten()))
+ _HX = Hm( _X )
+ _HX = numpy.asmatrix(_HX).flatten().T
+ GradJb = BI * (_X - Xb)
+ GradJo = - Ht( (_X, RI * (Y - _HX)) )
+ GradJ = numpy.asmatrix( GradJb ).flatten().T + numpy.asmatrix( GradJo ).flatten().T
+ 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
+ # ------------------------------------
+ if type(Xb) is type(numpy.matrix([])):
+ Xini = Xb.A1.tolist()
+ else:
+ Xini = list(Xb)
+ logging.debug("%s Point de démarrage Xini = %s"%(self._name, Xini))
+ #
+ # 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"]
+ else:
+ Bounds = None
+ MinimizerList = ["LBFGSB","TNC", "CG", "BFGS"]
+ if Par.has_key("Minimizer") and (Par["Minimizer"] in MinimizerList):
+ Minimizer = str( Par["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"] )
+ else:
+ maxiter = 15000
+ logging.debug("%s Nombre maximal de pas d'optimisation = %s"%(self._name, maxiter))
+ #
+ # Minimisation de la fonctionnelle
+ # --------------------------------
+ if Minimizer == "LBFGSB":
+ Minimum, J_optimal, Informations = scipy.optimize.fmin_l_bfgs_b(
+ func = CostFunction,
+ x0 = Xini,
+ fprime = GradientOfCostFunction,
+ args = (),
+ bounds = Bounds,
+ maxfun = maxiter,
+ iprint = iprint,
+ )
+ logging.debug("%s %s Minimum = %s"%(self._name, Minimizer, Minimum))
+ logging.debug("%s %s Nb of F = %s"%(self._name, Minimizer, Informations['funcalls']))
+ logging.debug("%s %s RetCode = %s"%(self._name, Minimizer, Informations['warnflag']))
+ elif Minimizer == "TNC":
+ Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
+ func = CostFunction,
+ x0 = Xini,
+ fprime = GradientOfCostFunction,
+ args = (),
+ bounds = Bounds,
+ maxfun = maxiter,
+ messages = message,
+ )
+ logging.debug("%s %s Minimum = %s"%(self._name, Minimizer, Minimum))
+ logging.debug("%s %s Nb of F = %s"%(self._name, Minimizer, nfeval))
+ logging.debug("%s %s RetCode = %s"%(self._name, Minimizer, rc))
+ elif Minimizer == "CG":
+ Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
+ f = CostFunction,
+ x0 = Xini,
+ fprime = GradientOfCostFunction,
+ args = (),
+ maxiter = maxiter,
+ disp = disp,
+ full_output = True,
+ )
+ logging.debug("%s %s Minimum = %s"%(self._name, Minimizer, Minimum))
+ logging.debug("%s %s Nb of F = %s"%(self._name, Minimizer, nfeval))
+ logging.debug("%s %s RetCode = %s"%(self._name, Minimizer, rc))
+ elif Minimizer == "BFGS":
+ Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
+ f = CostFunction,
+ x0 = Xini,
+ fprime = GradientOfCostFunction,
+ args = (),
+ maxiter = maxiter,
+ disp = disp,
+ full_output = True,
+ )
+ logging.debug("%s %s Minimum = %s"%(self._name, Minimizer, Minimum))
+ logging.debug("%s %s Nb of F = %s"%(self._name, Minimizer, nfeval))
+ logging.debug("%s %s RetCode = %s"%(self._name, Minimizer, rc))
+ else:
+ raise ValueError("Error in Minimizer name: %s"%Minimizer)
+ #
+ # Calcul de l'analyse
+ # --------------------
+ Xa = numpy.asmatrix(Minimum).T
+ logging.debug("%s Analyse Xa = %s"%(self._name, Xa))
+ #
+ self.StoredVariables["Analysis"].store( Xa.A1 )
+ self.StoredVariables["Innovation"].store( d.A1 )
+ #
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("MB")))
+ logging.debug("%s Terminé"%self._name)
+ #
+ return 0
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Algorithme de Kalman simple (BLUE)
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+import logging
+import Persistence
+from BasicObjects import Algorithm
+import PlatformInfo ; m = PlatformInfo.SystemUsage()
+
+# ==============================================================================
+class ElementaryAlgorithm(Algorithm):
+ def __init__(self):
+ Algorithm.__init__(self)
+ 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):
+ """
+ Calcul de l'estimateur BLUE (ou Kalman simple, ou Interpolation Optimale)
+ """
+ logging.debug("%s Lancement"%self._name)
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+ #
+ Hm = H["Direct"].asMatrix()
+ Ht = H["Adjoint"].asMatrix()
+ #
+ # Utilisation éventuelle d'un vecteur H(Xb) précalculé
+ # ----------------------------------------------------
+ if H["AppliedToX"] is not None and H["AppliedToX"].has_key("HXb"):
+ logging.debug("%s Utilisation de HXb"%self._name)
+ HXb = H["AppliedToX"]["HXb"]
+ else:
+ logging.debug("%s Calcul de Hm * Xb"%self._name)
+ HXb = Hm * Xb
+
+ # 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
+ else:
+ logging.debug("%s Calcul de K dans l'espace d'ébauche"%self._name)
+ K = (Ht * R.I * Hm + B.I).I * Ht * R.I
+ #
+ # Calcul de l'innovation et de l'analyse
+ # --------------------------------------
+ d = Y - HXb
+ logging.debug("%s Innovation d = %s"%(self._name, d))
+ Xa = Xb + K*d
+ logging.debug("%s Analyse Xa = %s"%(self._name, Xa))
+ #
+ self.StoredVariables["Analysis"].store( Xa.A1 )
+ self.StoredVariables["Innovation"].store( d.A1 )
+ #
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("MB")))
+ logging.debug("%s Terminé"%self._name)
+ #
+ return 0
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Algorithme de methode d'ensemble simple
+"""
+__author__ = "Sebastien MASSART, Jean-Philippe ARGAUD - Novembre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+import logging
+import numpy
+import Persistence
+from BasicObjects import Algorithm
+import PlatformInfo ; m = PlatformInfo.SystemUsage()
+
+# ==============================================================================
+class ElementaryAlgorithm(Algorithm):
+ def __init__(self):
+ Algorithm.__init__(self)
+ 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 ):
+ """
+ Calcul d'une estimation BLUE d'ensemble :
+ - génération d'un ensemble d'observations, de même taille que le
+ nombre d'ébauches
+ - calcul de l'estimateur BLUE pour chaque membre de l'ensemble
+ """
+ logging.debug("%s Lancement"%self._name)
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+ #
+ # Nombre d'ensemble pour l'ébauche
+ # --------------------------------
+ nb_ens = Xb.stepnumber()
+ #
+ # Construction de l'ensemble des observations, par génération a partir
+ # de la diagonale de R
+ # --------------------------------------------------------------------
+ DiagonaleR = numpy.diag(R)
+ EnsembleY = numpy.zeros([len(Y),nb_ens])
+ for npar in range(len(DiagonaleR)) :
+ bruit = numpy.random.normal(0,DiagonaleR[npar],nb_ens)
+ EnsembleY[npar,:] = Y[npar] + bruit
+ EnsembleY = numpy.matrix(EnsembleY)
+ #
+ # Initialisation des opérateurs d'observation et de la matrice gain
+ # -----------------------------------------------------------------
+ Hm = H["Direct"].asMatrix()
+ Ht = H["Adjoint"].asMatrix()
+
+ K = B * Ht * (Hm * B * Ht + R).I
+
+ # Calcul du BLUE pour chaque membre de l'ensemble
+ # -----------------------------------------------
+ for iens in range(nb_ens):
+ d = EnsembleY[:,iens] - Hm * Xb.valueserie(iens)
+ Xa = Xb.valueserie(iens) + K*d
+
+ self.StoredVariables["Analysis"].store( Xa.A1 )
+ self.StoredVariables["Innovation"].store( d.A1 )
+ #
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+ logging.debug("%s Terminé"%self._name)
+ return 0
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Algorithme de Kalman pour un système discret
+
+ Remarque : les observations sont exploitées à partir du pas de temps 1, et
+ sont utilisées dans Yo comme rangées selon ces indices. Donc le pas 0 n'est
+ pas utilisé puisque la première étape de Kalman passe de 0 à 1 avec
+ l'observation du pas 1.
+"""
+__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+import logging
+import Persistence
+from BasicObjects import Algorithm
+import PlatformInfo ; m = PlatformInfo.SystemUsage()
+
+# ==============================================================================
+class ElementaryAlgorithm(Algorithm):
+ def __init__(self):
+ Algorithm.__init__(self)
+ 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):
+ """
+ Calcul de l'estimateur de Kalman
+ """
+ logging.debug("%s Lancement"%self._name)
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+ #
+ # Opérateur d'observation
+ # -----------------------
+ Hm = H["Direct"].asMatrix()
+ Ht = H["Adjoint"].asMatrix()
+ #
+ # Opérateur d'évolution
+ # ---------------------
+ Mm = M["Direct"].asMatrix()
+ Mt = M["Adjoint"].asMatrix()
+ #
+ duration = Y.stepnumber()
+ #
+ # Initialisation
+ # --------------
+ Xn = Xb
+ Pn = B
+ self.StoredVariables["Analysis"].store( Xn.A1 )
+ self.StoredVariables["CovarianceAPosteriori"].store( Pn )
+ #
+ for step in range(duration-1):
+ logging.debug("%s Etape de Kalman %i (i.e. %i->%i) sur un total de %i"%(self._name, step+1, step,step+1, duration-1))
+ #
+ # Etape de prédiction
+ # -------------------
+ Xn_predicted = Mm * Xn
+ Pn_predicted = Mm * Pn * Mt + Q
+ #
+ # Etape de correction
+ # -------------------
+ d = Y.valueserie(step+1) - Hm * Xn_predicted
+ K = Pn_predicted * Ht * (Hm * Pn_predicted * Ht + R).I
+ Xn = Xn_predicted + K * d
+ Pn = Pn_predicted - K * Hm * Pn_predicted
+ #
+ self.StoredVariables["Analysis"].store( Xn.A1 )
+ self.StoredVariables["CovarianceAPosteriori"].store( Pn )
+ self.StoredVariables["Innovation"].store( d.A1 )
+ #
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+ logging.debug("%s Terminé"%self._name)
+ #
+ return 0
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Algorithme de moindre carres pondérés (analyse sans ebauche)
+"""
+__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Septembre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+import logging
+import Persistence
+from BasicObjects import Algorithm
+import PlatformInfo ; m = PlatformInfo.SystemUsage()
+
+# ==============================================================================
+class ElementaryAlgorithm(Algorithm):
+ def __init__(self):
+ 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):
+ """
+ Calcul de l'estimateur au sens des moindres carres sans ebauche
+ """
+ logging.debug("%s Lancement"%self._name)
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+ #
+ Hm = H["Direct"].asMatrix()
+ Ht = H["Adjoint"].asMatrix()
+ #
+ K = (Ht * R.I * Hm ).I * Ht * R.I
+ Xa = K * Y
+ #
+ self.StoredVariables["Analysis"].store( Xa.A1 )
+ #
+ logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+ logging.debug("%s Terminé"%self._name)
+ #
+ return 0
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+
--- /dev/null
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Définit les outils généraux élémentaires.
+
+ Ce module est destiné à etre appelée par AssimilationStudy pour constituer
+ les objets élémentaires de l'algorithme.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+import os, sys
+import numpy
+import Logging ; Logging.Logging() # A importer en premier
+import Persistence
+from BasicObjects import Operator
+
+# ==============================================================================
+class AssimilationStudy:
+ """
+ Cette classe sert d'interface pour l'utilisation de l'assimilation de
+ données. Elle contient les méthodes ou accesseurs nécessaires à la
+ construction d'un calcul d'assimilation.
+ """
+ def __init__(self, name=""):
+ """
+ Prévoit de conserver l'ensemble des variables nécssaires à un algorithme
+ élémentaire. Ces variables sont ensuite disponibles pour implémenter un
+ algorithme élémentaire particulier.
+
+ Background............: vecteur Xb
+ Observation...........: vecteur Y (potentiellement temporel)
+ d'observations
+ State.................: vecteur d'état dont une partie est le vecteur de
+ contrôle. Cette information n'est utile que si l'on veut faire des
+ calculs sur l'état complet, mais elle n'est pas indispensable pour
+ l'assimilation.
+ Control...............: vecteur X contenant toutes les variables de
+ contrôle, i.e. les paramètres ou l'état dont on veut estimer la
+ valeur pour obtenir les observations
+ ObservationOperator...: opérateur d'observation H
+
+ Les observations présentent une erreur dont la matrice de covariance est
+ R. L'ébauche du vecteur de contrôle présente une erreur dont la matrice
+ de covariance est B.
+ """
+ self.__name = str(name)
+ self.__Xb = None
+ self.__Y = None
+ self.__B = None
+ self.__R = None
+ self.__Q = None
+ self.__H = {}
+ self.__M = {}
+ #
+ self.__X = Persistence.OneVector()
+ self.__Parameters = {}
+ self.__StoredDiagnostics = {}
+ #
+ # Variables temporaires
+ 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,
+ # qui est activée dans Persistence)
+ self.__parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
+ sys.path.insert(0, self.__parent)
+ sys.path = list(set(sys.path)) # Conserve en unique exemplaire chaque chemin
+
+ # ---------------------------------------------------------
+ def setBackground(self,
+ asVector = None,
+ asPersistentVector = None,
+ Scheduler = None,
+ ):
+ """
+ Permet de définir l'estimation a priori :
+ - asVector : entrée des données, comme un vecteur compatible avec le
+ constructeur de numpy.matrix
+ - 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
+ """
+ if asVector is not None:
+ if type( asVector ) is type( numpy.matrix([]) ):
+ self.__Xb = numpy.matrix( asVector.A1, numpy.float ).T
+ else:
+ self.__Xb = numpy.matrix( asVector, numpy.float ).T
+ elif asPersistentVector is not None:
+ self.__Xb = asPersistentVector
+ else:
+ raise ValueError("Error: improperly defined background")
+ return 0
+
+ def setBackgroundError(self, asCovariance=None):
+ """
+ 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
+ """
+ self.__B = numpy.matrix( asCovariance, numpy.float )
+ return 0
+
+ # -----------------------------------------------------------
+ def setObservation(self,
+ asVector = None,
+ asPersistentVector = None,
+ Scheduler = None,
+ ):
+ """
+ Permet de définir les observations :
+ - asVector : entrée des données, comme un vecteur compatible avec le
+ constructeur de numpy.matrix
+ - 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
+ """
+ if asVector is not None:
+ if type( asVector ) is type( numpy.matrix([]) ):
+ self.__Y = numpy.matrix( asVector.A1, numpy.float ).T
+ else:
+ self.__Y = numpy.matrix( asVector, numpy.float ).T
+ elif asPersistentVector is not None:
+ self.__Y = asPersistentVector
+ else:
+ raise ValueError("Error: improperly defined observations")
+ return 0
+
+ def setObservationError(self, asCovariance=None):
+ """
+ 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
+ """
+ self.__R = numpy.matrix( asCovariance, numpy.float )
+ return 0
+
+ def setObservationOperator(self,
+ asFunction = {"Direct":None, "Tangent":None, "Adjoint":None},
+ asMatrix = None,
+ appliedToX = None,
+ ):
+ """
+ Permet de définir un opérateur d'observation H. L'ordre de priorité des
+ définitions et leur sens sont les suivants :
+ - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
+ alors on définit l'opérateur à l'aide de fonctions. Si la fonction
+ "Direct" n'est pas définie, on prend la fonction "Tangent".
+ - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
+ None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
+ 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.
+ - si l'argument "appliedToX" n'est pas None, alors on définit, pour des
+ 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.
+ """
+ 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):
+ self.__H["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
+ else:
+ self.__H["Direct"] = Operator( fromMethod = asFunction["Direct"] )
+ self.__H["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
+ self.__H["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
+ elif asMatrix is not None:
+ mat = numpy.matrix( asMatrix, numpy.float )
+ self.__H["Direct"] = Operator( fromMatrix = mat )
+ self.__H["Tangent"] = Operator( fromMatrix = mat )
+ self.__H["Adjoint"] = Operator( fromMatrix = mat.T )
+ else:
+ raise ValueError("Error: improperly defined observation operator")
+ #
+ if appliedToX is not None:
+ self.__H["AppliedToX"] = {}
+ if type(appliedToX) is not dict:
+ raise ValueError("Error: observation operator defined by \"appliedToX\" need a dictionary as argument.")
+ for key in appliedToX.keys():
+ if type( appliedToX[key] ) is type( numpy.matrix([]) ):
+ # Pour le cas où l'on a une vraie matrice
+ self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].A1, numpy.float ).T
+ elif type( appliedToX[key] ) is type( numpy.array([]) ) and len(appliedToX[key].shape) > 1:
+ # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
+ self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key].reshape(len(appliedToX[key]),), numpy.float ).T
+ else:
+ self.__H["AppliedToX"][key] = numpy.matrix( appliedToX[key], numpy.float ).T
+ else:
+ self.__H["AppliedToX"] = None
+ #
+ return 0
+
+ # -----------------------------------------------------------
+ def setEvolutionModel(self,
+ asFunction = {"Direct":None, "Tangent":None, "Adjoint":None},
+ asMatrix = None,
+ Scheduler = None,
+ ):
+ """
+ Permet de définir un opérateur d'évolution M. L'ordre de priorité des
+ définitions et leur sens sont les suivants :
+ - si asFunction["Tangent"] et asFunction["Adjoint"] ne sont pas None
+ alors on définit l'opérateur à l'aide de fonctions. Si la fonction
+ "Direct" n'est pas définie, on prend la fonction "Tangent".
+ - si les fonctions ne sont pas disponibles et si asMatrix n'est pas
+ None, alors on définit l'opérateur "Direct" et "Tangent" à l'aide de
+ 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.
+ """
+ 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):
+ self.__M["Direct"] = Operator( fromMethod = asFunction["Tangent"] )
+ else:
+ self.__M["Direct"] = Operator( fromMethod = asFunction["Direct"] )
+ self.__M["Tangent"] = Operator( fromMethod = asFunction["Tangent"] )
+ self.__M["Adjoint"] = Operator( fromMethod = asFunction["Adjoint"] )
+ elif asMatrix is not None:
+ matrice = numpy.matrix( asMatrix, numpy.float )
+ self.__M["Direct"] = Operator( fromMatrix = matrice )
+ self.__M["Tangent"] = Operator( fromMatrix = matrice )
+ self.__M["Adjoint"] = Operator( fromMatrix = matrice.T )
+ else:
+ raise ValueError("Error: improperly defined evolution operator")
+ return 0
+
+ def setEvolutionError(self, asCovariance=None):
+ """
+ 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
+ """
+ self.__Q = numpy.matrix( asCovariance, numpy.float )
+ return 0
+
+ # -----------------------------------------------------------
+ def setControls (self, asVector = None ):
+ """
+ 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
+ estimer la valeur pour obtenir les observations. C'est utile pour un
+ algorithme itératif/incrémental
+ - asVector : entrée des données, comme un vecteur compatible avec le
+ constructeur de numpy.matrix.
+ """
+ if asVector is not None:
+ self.__X.store( asVector )
+ return 0
+
+ # -----------------------------------------------------------
+ def setAlgorithm(self, choice = None ):
+ """
+ Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
+ d'assimilation. L'argument est un champ caractère se rapportant au nom
+ d'un fichier contenu dans "../daAlgorithms" et réalisant l'opération
+ d'assimilation sur les arguments (Xb,Y,H,R,B,Xa).
+ """
+ if choice is None:
+ raise ValueError("Error: algorithm choice has to be given")
+ if self.__algorithmName is not None:
+ raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
+ daDirectory = "daAlgorithms"
+ #
+ # Recherche explicitement le fichier complet
+ # ------------------------------------------
+ module_path = None
+ for directory in sys.path:
+ if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
+ module_path = os.path.abspath(os.path.join(directory, daDirectory))
+ if module_path is None:
+ raise ImportError("No algorithm module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
+ #
+ # Importe le fichier complet comme un module
+ # ------------------------------------------
+ try:
+ sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
+ self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
+ self.__algorithmName = str(choice)
+ sys.path = sys_path_tmp ; del sys_path_tmp
+ except ImportError, e:
+ raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
+ #
+ # Instancie un objet du type élémentaire du fichier
+ # -------------------------------------------------
+ self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
+ return 0
+
+ def setAlgorithmParameters(self, asDico=None):
+ """
+ Permet de définir les paramètres de l'algorithme, sous la forme d'un
+ dictionnaire.
+ """
+ self.__Parameters = dict( asDico )
+ return 0
+
+ # -----------------------------------------------------------
+ def setDiagnostic(self, choice = None, name = "", unit = "", basetype = None, parameters = {} ):
+ """
+ Permet de sélectionner un diagnostic a effectuer.
+ """
+ if choice is None:
+ raise ValueError("Error: diagnostic choice has to be given")
+ daDirectory = "daDiagnostics"
+ #
+ # Recherche explicitement le fichier complet
+ # ------------------------------------------
+ module_path = None
+ for directory in sys.path:
+ if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
+ module_path = os.path.abspath(os.path.join(directory, daDirectory))
+ if module_path is None:
+ raise ImportError("No diagnostic module named \"%s\" was found in a \"%s\" subdirectory\n The search path is %s"%(choice, daDirectory, sys.path))
+ #
+ # Importe le fichier complet comme un module
+ # ------------------------------------------
+ try:
+ sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
+ self.__diagnosticFile = __import__(str(choice), globals(), locals(), [])
+ sys.path = sys_path_tmp ; del sys_path_tmp
+ except ImportError, e:
+ raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n The import error message is: %s"%(choice,e))
+ #
+ # 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")
+ else:
+ self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
+ name = name,
+ unit = unit,
+ basetype = basetype,
+ parameters = parameters )
+ return 0
+
+ # -----------------------------------------------------------
+ def shape_validate(self):
+ """
+ Validation de la correspondance correcte des tailles des variables et
+ des matrices s'il y en a.
+ """
+ if self.__Xb is None: __Xb_shape = (0,)
+ elif hasattr(self.__Xb,"shape"):
+ if type(self.__Xb.shape) is tuple: __Xb_shape = self.__Xb.shape
+ else: __Xb_shape = self.__Xb.shape()
+ else: raise TypeError("Xb has no attribute of shape: problem !")
+ #
+ if self.__Y is None: __Y_shape = (0,)
+ elif hasattr(self.__Y,"shape"):
+ if type(self.__Y.shape) is tuple: __Y_shape = self.__Y.shape
+ else: __Y_shape = self.__Y.shape()
+ else: raise TypeError("Y has no attribute of shape: problem !")
+ #
+ if self.__B is None: __B_shape = (0,0)
+ elif hasattr(self.__B,"shape"):
+ if type(self.__B.shape) is tuple: __B_shape = self.__B.shape
+ else: __B_shape = self.__B.shape()
+ else: raise TypeError("B has no attribute of shape: problem !")
+ #
+ if self.__R is None: __R_shape = (0,0)
+ elif hasattr(self.__R,"shape"):
+ if type(self.__R.shape) is tuple: __R_shape = self.__R.shape
+ else: __R_shape = self.__R.shape()
+ else: raise TypeError("R has no attribute of shape: problem !")
+ #
+ if self.__Q is None: __Q_shape = (0,0)
+ elif hasattr(self.__Q,"shape"):
+ if type(self.__Q.shape) is tuple: __Q_shape = self.__Q.shape
+ else: __Q_shape = self.__Q.shape()
+ else: raise TypeError("Q has no attribute of shape: problem !")
+ #
+ if len(self.__H) == 0: __H_shape = (0,0)
+ elif type(self.__H) is type({}): __H_shape = (0,0)
+ elif hasattr(self.__H["Direct"],"shape"):
+ if type(self.__H["Direct"].shape) is tuple: __H_shape = self.__H["Direct"].shape
+ else: __H_shape = self.__H["Direct"].shape()
+ else: raise TypeError("H has no attribute of shape: problem !")
+ #
+ if len(self.__M) == 0: __M_shape = (0,0)
+ elif type(self.__M) is type({}): __M_shape = (0,0)
+ elif hasattr(self.__M["Direct"],"shape"):
+ if type(self.__M["Direct"].shape) is tuple: __M_shape = self.__M["Direct"].shape
+ else: __M_shape = self.__M["Direct"].shape()
+ else: raise TypeError("M has no attribute of shape: problem !")
+ #
+ # Vérification des conditions
+ # ---------------------------
+ if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
+ raise ValueError("Shape characteristic of Xb is incorrect: \"%s\""%(__Xb_shape,))
+ if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
+ raise ValueError("Shape characteristic of Y is incorrect: \"%s\""%(__Y_shape,))
+ #
+ if not( min(__B_shape) == max(__B_shape) ):
+ raise ValueError("Shape characteristic of B is incorrect: \"%s\""%(__B_shape,))
+ if not( min(__R_shape) == max(__R_shape) ):
+ raise ValueError("Shape characteristic of R is incorrect: \"%s\""%(__R_shape,))
+ if not( min(__Q_shape) == max(__Q_shape) ):
+ raise ValueError("Shape characteristic of Q is incorrect: \"%s\""%(__Q_shape,))
+ if not( min(__M_shape) == max(__M_shape) ):
+ raise ValueError("Shape characteristic of M is incorrect: \"%s\""%(__M_shape,))
+ #
+ if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[1] == max(__Xb_shape) ):
+ raise ValueError("Shape characteristic of H \"%s\" and X \"%s\" are incompatible"%(__H_shape,__Xb_shape))
+ if len(self.__H) > 0 and not(type(self.__H) is type({})) and not( __H_shape[0] == max(__Y_shape) ):
+ raise ValueError("Shape characteristic of H \"%s\" and Y \"%s\" are incompatible"%(__H_shape,__Y_shape))
+ if len(self.__H) > 0 and not(type(self.__H) is type({})) and len(self.__B) > 0 and not( __H_shape[1] == __B_shape[0] ):
+ raise ValueError("Shape characteristic of H \"%s\" and B \"%s\" are incompatible"%(__H_shape,__B_shape))
+ if len(self.__H) > 0 and not(type(self.__H) is type({})) and len(self.__R) > 0 and not( __H_shape[0] == __R_shape[1] ):
+ raise ValueError("Shape characteristic of H \"%s\" and R \"%s\" are incompatible"%(__H_shape,__R_shape))
+ #
+ if len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
+ raise ValueError("Shape characteristic of B \"%s\" and Xb \"%s\" are incompatible"%(__B_shape,__Xb_shape))
+ #
+ if len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
+ raise ValueError("Shape characteristic of R \"%s\" and Y \"%s\" are incompatible"%(__R_shape,__Y_shape))
+ #
+ if len(self.__M) > 0 and not(type(self.__M) is type({})) and not( __M_shape[1] == max(__Xb_shape) ):
+ raise ValueError("Shape characteristic of M \"%s\" and X \"%s\" are incompatible"%(__M_shape,__Xb_shape))
+ #
+ return 1
+
+ # -----------------------------------------------------------
+ def analyze(self):
+ """
+ Permet de lancer le calcul d'assimilation.
+
+ Le nom de la méthode à activer est toujours "run". Les paramètres en
+ arguments de la méthode sont fixés. En sortie, on obtient les résultats
+ dans la variable de type dictionnaire "StoredVariables", qui contient en
+ particulier des objets de Persistence pour les analyses, OMA...
+ """
+ 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,
+ )
+ return 0
+
+ # -----------------------------------------------------------
+ def get(self, key=None):
+ """
+ Renvoie les résultats disponibles après l'exécution de la méthode
+ d'assimilation, ou les diagnostics disponibles. Attention, quand un
+ diagnostic porte le même nom qu'un variable stockée, c'est la variable
+ stockée qui est renvoyée, et le diagnostic est inatteignable.
+ """
+ if key is not None:
+ if self.__algorithm.has_key(key):
+ return self.__algorithm.get( 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)
+ else:
+ allvariables = self.__algorithm.get()
+ allvariables.update( self.__StoredDiagnostics )
+ return allvariables
+
+ def get_available_algorithms(self):
+ """
+ Renvoie la liste des algorithmes identifiés par les chaînes de
+ caractères
+ """
+ files = []
+ for directory in sys.path:
+ if os.path.isdir(os.path.join(directory,"daAlgorithms")):
+ for fname in os.listdir(os.path.join(directory,"daAlgorithms")):
+ root, ext = os.path.splitext(fname)
+ if ext == '.py' and root != '__init__':
+ files.append(root)
+ files.sort()
+ return files
+
+ def get_available_diagnostics(self):
+ """
+ Renvoie la liste des diagnostics identifiés par les chaînes de
+ caractères
+ """
+ files = []
+ for directory in sys.path:
+ if os.path.isdir(os.path.join(directory,"daDiagnostics")):
+ for fname in os.listdir(os.path.join(directory,"daDiagnostics")):
+ root, ext = os.path.splitext(fname)
+ if ext == '.py' and root != '__init__':
+ files.append(root)
+ files.sort()
+ return files
+
+ # -----------------------------------------------------------
+ def get_algorithms_main_path(self):
+ """
+ Renvoie le chemin pour le répertoire principal contenant les algorithmes
+ dans un sous-répertoire "daAlgorithms"
+ """
+ return self.__parent
+
+ def add_algorithms_path(self, asPath=None):
+ """
+ Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
+ se trouve un sous-répertoire "daAlgorithms"
+
+ Remarque : si le chemin a déjà été ajouté pour les diagnostics, il n'est
+ pas indispensable de le rajouter ici.
+ """
+ if not os.path.isdir(asPath):
+ raise ValueError("The given "+asPath+" argument must exist as a directory")
+ if not os.path.isdir(os.path.join(asPath,"daAlgorithms")):
+ raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daAlgorithms\"")
+ if not os.path.isfile(os.path.join(asPath,"daAlgorithms","__init__.py")):
+ raise ValueError("The given \""+asPath+"/daAlgorithms\" path must contain a file named \"__init__.py\"")
+ sys.path.insert(0, os.path.abspath(asPath))
+ sys.path = list(set(sys.path)) # Conserve en unique exemplaire chaque chemin
+ return 1
+
+ def get_diagnostics_main_path(self):
+ """
+ Renvoie le chemin pour le répertoire principal contenant les diagnostics
+ dans un sous-répertoire "daDiagnostics"
+ """
+ return self.__parent
+
+ def add_diagnostics_path(self, asPath=None):
+ """
+ Ajoute au chemin de recherche des algorithmes un répertoire dans lequel
+ se trouve un sous-répertoire "daDiagnostics"
+
+ Remarque : si le chemin a déjà été ajouté pour les algorithmes, il n'est
+ pas indispensable de le rajouter ici.
+ """
+ if not os.path.isdir(asPath):
+ raise ValueError("The given "+asPath+" argument must exist as a directory")
+ if not os.path.isdir(os.path.join(asPath,"daDiagnostics")):
+ raise ValueError("The given \""+asPath+"\" argument must contain a subdirectory named \"daDiagnostics\"")
+ if not os.path.isfile(os.path.join(asPath,"daDiagnostics","__init__.py")):
+ raise ValueError("The given \""+asPath+"/daDiagnostics\" path must contain a file named \"__init__.py\"")
+ sys.path.insert(0, os.path.abspath(asPath))
+ sys.path = list(set(sys.path)) # Conserve en unique exemplaire chaque chemin
+ return 1
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ ADD = AssimilationStudy("Ma premiere etude BLUE")
+
+ ADD.setBackground (asVector = [0, 1, 2])
+ ADD.setBackgroundError (asCovariance = "1 0 0;0 1 0;0 0 1")
+ ADD.setObservation (asVector = [0.5, 1.5, 2.5])
+ ADD.setObservationError (asCovariance = "1 0 0;0 1 0;0 0 1")
+ ADD.setObservationOperator(asMatrix = "1 0 0;0 1 0;0 0 1")
+
+ ADD.setAlgorithm(choice="Blue")
+
+ ADD.analyze()
+
+ print "Nombre d'analyses :", ADD.get("Analysis").stepnumber()
+ print "Analyse résultante :", ADD.get("Analysis").valueserie(0)
+ 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
+
+ ADD.setDiagnostic("RMS", "Ma RMS")
+
+ liste = ADD.get().keys()
+ liste.sort()
+ print "Variables et diagnostics disponibles :", liste
+ print
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Définit les outils généraux élémentaires.
+
+ Ce module est destiné à etre appelée par AssimilationStudy pour constituer
+ les objets élémentaires de l'algorithme.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+import numpy
+import Persistence
+
+# ==============================================================================
+class Operator:
+ """
+ Classe générale d'interface de type opérateur
+ """
+ def __init__(self, fromMethod=None, fromMatrix=None):
+ """
+ On construit un objet de ce type en fournissant à l'aide de l'un des
+ deux mots-clé, soit une fonction python, soit matrice.
+ Arguments :
+ - fromMethod : argument de type fonction Python
+ - fromMatrix : argument adapté au constructeur numpy.matrix
+ """
+ if fromMethod is not None:
+ self.__Method = fromMethod
+ self.__Matrix = None
+ elif fromMatrix is not None:
+ self.__Method = None
+ self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
+ else:
+ self.__Method = None
+ self.__Matrix = None
+
+ def appliedTo(self, xValue):
+ """
+ Permet de restituer le résultat de l'application de l'opérateur à un
+ argument xValue. Cette méthode se contente d'appliquer, son argument
+ devant a priori être du bon type.
+ Arguments :
+ - xValue : argument adapté pour appliquer l'opérateur
+ """
+ if self.__Matrix is not None:
+ return self.__Matrix * xValue
+ else:
+ return self.__Method( xValue )
+
+ def appliedInXTo(self, (xNominal, xValue) ):
+ """
+ Permet de restituer le résultat de l'application de l'opérateur à un
+ argument xValue, sachant que l'opérateur est valable en xNominal.
+ Cette méthode se contente d'appliquer, son argument devant a priori
+ être du bon type. Si l'opérateur est linéaire car c'est une matrice,
+ alors il est valable en tout point nominal et il n'est pas nécessaire
+ d'utiliser xNominal.
+ Arguments : une liste contenant
+ - xNominal : argument permettant de donner le point où l'opérateur
+ est construit pour etre ensuite appliqué
+ - xValue : argument adapté pour appliquer l'opérateur
+ """
+ if self.__Matrix is not None:
+ return self.__Matrix * xValue
+ else:
+ return self.__Method( (xNominal, xValue) )
+
+ def asMatrix(self):
+ """
+ Permet de renvoyer l'opérateur sous la forme d'une matrice
+ """
+ if self.__Matrix is not None:
+ return self.__Matrix
+ else:
+ raise ValueError("Matrix form of the operator is not available")
+
+ def shape(self):
+ """
+ Renvoie la taille sous forme numpy si l'opérateur est disponible sous
+ la forme d'une matrice
+ """
+ if self.__Matrix is not None:
+ return self.__Matrix.shape
+ else:
+ raise ValueError("Matrix form of the operator is not available, nor the shape")
+
+# ==============================================================================
+class Algorithm:
+ """
+ Classe générale d'interface de type algorithme
+
+ Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
+ d'assimilation, en fournissant un container (dictionnaire) de variables
+ persistantes initialisées, et des méthodes d'accès à ces variables stockées.
+
+ Une classe élémentaire d'algorithme doit implémenter la méthode "run".
+ """
+ def __init__(self):
+ """
+ L'initialisation présente permet de fabriquer des variables de stockage
+ disponibles de manière générique dans les algorithmes élémentaires. Ces
+ variables de stockage sont ensuite conservées dans un dictionnaire
+ interne à l'objet, mais auquel on accède par la méthode "get".
+
+ Les variables prévues sont :
+ - Analysis : l'analyse
+ - Innovation : l'innovation : d = Y - H Xb
+ - SigmaObs2 : correction optimale des erreurs d'observation
+ - SigmaBck2 : correction optimale des erreurs d'ébauche
+ - OMA : Observation moins Analysis : Y - Xa
+ - OMB : Observation moins Background : Y - Xb
+ - AMB : Analysis moins Background : Xa - Xb
+ - CovarianceAPosteriori : matrice A
+ On peut rajouter des variables à stocker dans l'initialisation de
+ l'algorithme élémentaire qui va hériter de cette classe
+ """
+ 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["Analysis"] = Persistence.OneVector(name = "Analysis")
+ self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
+ self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
+ self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
+ self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
+ 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):
+ """
+ Renvoie l'une des variables stockées identifiée par la clé, ou le
+ dictionnaire de l'ensemble des variables disponibles en l'absence de
+ clé. Ce sont directement les variables sous forme objet qui sont
+ renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
+ des classes de persistance.
+ """
+ if key is not None:
+ return self.StoredVariables[key]
+ else:
+ return self.StoredVariables
+
+ def has_key(self, key=None):
+ """
+ Vérifie si l'une des variables stockées est identifiée par la clé.
+ """
+ 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):
+ """
+ Doit implémenter l'opération élémentaire de calcul d'assimilation sous
+ sa forme mathématique la plus naturelle possible.
+ """
+ raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
+
+# ==============================================================================
+class Diagnostic:
+ """
+ Classe générale d'interface de type diagnostic
+
+ Ce template s'utilise de la manière suivante : il sert de classe "patron" en
+ même temps que l'une des classes de persistance, comme "OneScalar" par
+ exemple.
+
+ Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
+ méthode "_formula" pour écrire explicitement et proprement la formule pour
+ l'écriture mathématique du calcul du diagnostic (méthode interne non
+ publique), et "calculate" pour activer la précédente tout en ayant vérifié
+ et préparé les données, et pour stocker les résultats à chaque pas (méthode
+ externe d'activation).
+ """
+ def __init__(self, name = "", parameters = {}):
+ self.name = str(name)
+ self.parameters = dict( parameters )
+
+ def _formula(self, *args):
+ """
+ Doit implémenter l'opération élémentaire de diagnostic sous sa forme
+ mathématique la plus naturelle possible.
+ """
+ raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
+
+ def calculate(self, *args):
+ """
+ Active la formule de calcul avec les arguments correctement rangés
+ """
+ raise NotImplementedError("Diagnostic activation method has not been implemented!")
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Ce module permet de mettre en place un logging utilisable partout dans
+ l'application, par défaut à la console, et si nécessaire dans un fichier.
+
+ Il doit être appelé en premier dans AssimilationStudy (mais pas directement
+ dans les applications utilisateurs), en l'important et en instanciant un
+ objet :
+ import Logging ; Logging.Logging()
+
+ Par défaut, seuls les messages du niveau WARNING ou au-delà sont disponibles
+ (donc les simples messages d'info ne sont pas disponibles), ce que l'on peut
+ changer à l'instanciation avec le mot-clé "level" :
+ import Logging ; Logging.Logging(level=20)
+
+ On peut éventuellement demander à l'objet de sortir aussi les messages dans
+ un fichier (noms par défaut : AssimilationStudy.log, niveau NOTSET) :
+ import Logging ; Logging.Logging().setLogfile()
+
+ Si on veut changer le nom du fichier ou le niveau global de message, il faut
+ récupérer l'instance et appliquer les méthodes :
+ import Logging
+ log = Logging.Logging()
+ import logging
+ log.setLevel(logging.DEBUG)
+ log.setLogfile(filename="toto.log", filemode="a", level=logging.WARNING)
+ et on change éventuellement le niveau avec :
+ log.setLogfileLevel(logging.INFO)
+
+ Ensuite, n'importe où dans les applications, il suffit d'utiliser le module
+ "logging" (avec un petit "l") :
+ import logging
+ log = logging.getLogger(NAME) # Avec rien (recommandé) ou un nom NAME
+ log.critical("...")
+ log.error("...")
+ log.warning("...")
+ log.info("...")
+ log.debug("...")
+ ou encore plus simplement :
+ import logging
+ logging.info("...")
+
+ Dans une application, à n'importe quel endroit et autant de fois qu'on veut,
+ on peut changer le niveau global de message en utilisant par exemple :
+ import logging
+ logging.setLevel(logging.DEBUG)
+
+ On rappelle les niveaux (attributs de "logging") et leur ordre :
+ NOTSET=0 < DEBUG=10 < INFO=20 < WARNING=30 < ERROR=40 < CRITICAL=50
+"""
+__author__ = "Jean-Philippe ARGAUD - Octobre 2008"
+
+import os
+import sys
+import logging
+from PlatformInfo import PlatformInfo
+
+LOGFILE = os.path.join(os.path.abspath(os.curdir),"AssimilationStudy.log")
+
+# ==============================================================================
+class Logging:
+ def __init__(self, level=logging.WARNING):
+ """
+ Initialise un logging à la console pour TOUS les niveaux de messages.
+ """
+ logging.basicConfig(
+ format = '%(levelname)-8s %(message)s',
+ level = level,
+ stream = sys.stdout,
+ )
+ self.__logfile = None
+ #
+ # Initialise l'affichage de logging
+ # ---------------------------------
+ p = PlatformInfo()
+ #
+ logging.info( "--------------------------------------------------" )
+ logging.info( "Lancement de "+p.getName()+" "+p.getVersion() )
+ logging.info( "--------------------------------------------------" )
+ logging.info( "Versions logicielles :" )
+ logging.info( "- Python "+p.getPythonVersion() )
+ logging.info( "- Numpy "+p.getNumpyVersion() )
+ logging.info( "- Scipy "+p.getScipyVersion() )
+ logging.info( "" )
+
+ def setLogfileLevel(self, level=logging.NOTSET ):
+ """
+ Permet de changer globalement le niveau des messages disponibles.
+ """
+ logging.getLogger().setLevel(level)
+
+ def setLogfile(self, filename=LOGFILE, filemode="w", level=logging.NOTSET):
+ """
+ Permet de disposer des messages dans un fichier EN PLUS de la console.
+ """
+ if self.__logfile is not None:
+ # Supprime le précédent mode de stockage fichier s'il exsitait
+ logging.getLogger().removeHandler(self.__logfile)
+ self.__logfile = logging.FileHandler(filename, filemode)
+ self.__logfile.setLevel(level)
+ self.__logfile.setFormatter(
+ logging.Formatter('%(asctime)s %(levelname)-8s %(message)s',
+ '%d %b %Y %H:%M:%S'))
+ logging.getLogger().addHandler(self.__logfile)
+
+ def setLogfileLevel(self, level=logging.NOTSET ):
+ """
+ Permet de changer le niveau des messages stockés en fichier. Il ne sera
+ pris en compte que s'il est supérieur au niveau global.
+ """
+ self.__logfile.setLevel(level)
+
+ def getLevel(self):
+ """
+ Renvoie le niveau de Logging sous forme texte
+ """
+ return logging.getLevelName( logging.getLogger().getEffectiveLevel() )
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+ import os.path
+
+ l = Logging(level = logging.NOTSET)
+
+ logging.info("Message numéro 1 uniquement disponible sur console")
+
+ l.setLogfile(level = logging.WARNING)
+ if not os.path.isfile(LOGFILE):
+ raise ValueError("Le fichier de log \"%s\" n'a pas pu être créé."%LOGFILE)
+
+ logging.info("Message numéro 2 uniquement disponible sur console")
+ logging.warning("Message numéro 3 conjointement disponible sur console et fichier")
+
+ l.setLogfileLevel(logging.INFO)
+
+ logging.info("Message numéro 4 conjointement disponible sur console et fichier")
+
+ print
+ print " Le logging a été correctement initialisé. Le fichier suivant"
+ print " %s"%os.path.basename(LOGFILE)
+ print " a été correctement créé, et peut être effacé après vérification."
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Définit des outils de persistence et d'enregistrement de séries de valeurs
+ pour analyse ultérieure ou utilisation de calcul.
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+import numpy
+
+from PlatformInfo import PathManagement ; PathManagement()
+
+# ==============================================================================
+class Persistence:
+ """
+ Classe générale de persistence définissant les accesseurs nécessaires
+ (Template)
+ """
+ def __init__(self, name="", unit="", basetype=str):
+ """
+ name : nom courant
+ unit : unité
+ basetype : type de base de l'objet stocké à chaque pas
+
+ 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) :
+ __step : numérotation par défaut du pas courant
+ __basetype : le type de base de chaque valeur, sous la forme d'un type
+ permettant l'instanciation ou le casting Python
+ __steps : les pas de stockage. Par défaut, c'est __step
+ __values : les valeurs de stockage. Par défaut, c'est None
+ """
+ self.__name = str(name)
+ self.__unit = str(unit)
+ #
+ self.__step = -1
+ self.__basetype = basetype
+ #
+ self.__steps = []
+ self.__values = []
+
+ def basetype(self, basetype=None):
+ """
+ Renvoie ou met en place le type de base des objets stockés
+ """
+ if basetype is None:
+ return self.__basetype
+ else:
+ self.__basetype = basetype
+
+ def store(self, value=None, step=None):
+ """
+ Stocke une valeur à un pas. Une instanciation est faite avec le type de
+ base pour stocker l'objet. Si le pas n'est pas fournit, on utilise
+ l'étape de stockage comme valeur de pas.
+ """
+ if value is None: raise ValueError("Value argument required")
+ self.__step += 1
+ if step is not None:
+ self.__steps.append(step)
+ else:
+ self.__steps.append(self.__step)
+ #
+ self.__values.append(self.__basetype(value))
+
+ def shape(self):
+ """
+ Renvoie la taille sous forme numpy du dernier objet stocké. Si c'est un
+ objet numpy, renvoie le shape. Si c'est un entier, un flottant, un
+ complexe, renvoie 1. Si c'est une liste ou un dictionnaire, renvoie la
+ longueur. Par défaut, renvoie 1.
+ """
+ if len(self.__values) > 0:
+ if self.__basetype in [numpy.matrix, numpy.array]:
+ return self.__values[-1].shape
+ elif self.__basetype in [int, float]:
+ return (1,)
+ elif self.__basetype in [list, dict]:
+ return (len(self.__values[-1]),)
+ else:
+ return (1,)
+ else:
+ raise ValueError("Object has no shape before its first storage")
+
+ def __len__(self):
+ """
+ Renvoie le nombre d'éléments dans un séquence ou la plus grande
+ dimension d'une matrice
+ """
+ return max( self.shape() )
+
+ # ---------------------------------------------------------
+ def stepserie(self, item=None, step=None):
+ """
+ Renvoie par défaut toute la liste des pas de temps. Si l'argument "step"
+ existe dans la liste des pas de stockage effectués, renvoie ce pas
+ "step". Si l'argument "item" est correct, renvoie le pas stockée au
+ numéro "item".
+ """
+ if step is not None and step in self.__steps:
+ return step
+ elif item is not None and item < len(self.__steps):
+ return self.__steps[item]
+ else:
+ return self.__steps
+
+ def valueserie(self, item=None, step=None):
+ """
+ Renvoie par défaut toute la liste des valeurs/objets. Si l'argument
+ "step" existe dans la liste des pas de stockage effectués, renvoie la
+ valeur stockée à ce pas "step". Si l'argument "item" est correct,
+ renvoie la valeur stockée au numéro "item".
+ """
+ if step is not None and step in self.__steps:
+ index = self.__steps.index(step)
+ return self.__values[index]
+ elif item is not None and item < len(self.__values):
+ return self.__values[item]
+ else:
+ return self.__values
+
+ def stepnumber(self):
+ """
+ Renvoie le nombre de pas de stockage.
+ """
+ return len(self.__steps)
+
+ # ---------------------------------------------------------
+ def mean(self):
+ """
+ Renvoie la valeur moyenne des données à chaque pas. Il faut que le type
+ de base soit compatible avec les types élémentaires numpy.
+ """
+ try:
+ return [numpy.matrix(item).mean() for item in self.__values]
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def std(self, ddof=0):
+ """
+ Renvoie l'écart-type des données à chaque pas. Il faut que le type de
+ base soit compatible avec les types élémentaires numpy.
+
+ ddof : c'est le nombre de degrés de liberté pour le calcul de
+ l'écart-type, qui est dans le diviseur. Inutile avant Numpy 1.1
+ """
+ try:
+ if numpy.version.version >= '1.1.0':
+ return [numpy.matrix(item).std(ddof=ddof) for item in self.__values]
+ else:
+ return [numpy.matrix(item).std() for item in self.__values]
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def sum(self):
+ """
+ Renvoie la somme des données à chaque pas. Il faut que le type de
+ base soit compatible avec les types élémentaires numpy.
+ """
+ try:
+ return [numpy.matrix(item).sum() for item in self.__values]
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def min(self):
+ """
+ Renvoie le minimum des données à chaque pas. Il faut que le type de
+ base soit compatible avec les types élémentaires numpy.
+ """
+ try:
+ return [numpy.matrix(item).min() for item in self.__values]
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def max(self):
+ """
+ Renvoie le maximum des données à chaque pas. Il faut que le type de
+ base soit compatible avec les types élémentaires numpy.
+ """
+ try:
+ return [numpy.matrix(item).max() for item in self.__values]
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def plot(self, item=None, step=None,
+ steps = None,
+ title = "",
+ xlabel = "",
+ ylabel = "",
+ ltitle = None,
+ geometry = "600x400",
+ filename = "",
+ persist = False,
+ pause = True,
+ ):
+ """
+ Renvoie un affichage de la valeur à chaque pas, si elle est compatible
+ avec un affichage Gnuplot (donc essentiellement un vecteur). Si
+ l'argument "step" existe dans la liste des pas de stockage effectués,
+ renvoie l'affichage de la valeur stockée à ce pas "step". Si l'argument
+ "item" est correct, renvoie l'affichage de la valeur stockée au numéro
+ "item". Par défaut ou en l'absence de "step" ou "item", renvoie un
+ affichage successif de tous les pas.
+
+ Arguments :
+ - step : valeur du pas à afficher
+ - item : index de la valeur à afficher
+ - steps : liste unique des pas de l'axe des X, ou None si c'est
+ la numérotation par défaut
+ - title : base du titre général, qui sera automatiquement
+ complétée par la mention du pas
+ - xlabel : label de l'axe des X
+ - ylabel : label de l'axe des Y
+ - ltitle : titre associé au vecteur tracé
+ - geometry : taille en pixels de la fenêtre et position du coin haut
+ gauche, au format X11 : LxH+X+Y (défaut : 600x400)
+ - filename : base de nom de fichier Postscript pour une sauvegarde,
+ qui est automatiquement complétée par le numéro du
+ fichier calculé par incrément simple de compteur
+ - persist : booléen indiquant que la fenêtre affichée sera
+ conservée lors du passage au dessin suivant
+ Par défaut, persist = False
+ - pause : booléen indiquant une pause après chaque tracé, et
+ attendant un Return
+ Par défaut, pause = True
+ """
+ import os
+ #
+ # Vérification de la disponibilité du module Gnuplot
+ try:
+ import Gnuplot
+ self.__gnuplot = Gnuplot
+ except:
+ raise ImportError("The Gnuplot module is required to plot the object.")
+ #
+ # Vérification et compléments sur les paramètres d'entrée
+ if persist:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -persist -geometry '+geometry
+ else:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -geometry '+geometry
+ if ltitle is None:
+ ltitle = ""
+ self.__g = self.__gnuplot.Gnuplot() # persist=1
+ self.__g('set terminal '+self.__gnuplot.GnuplotOpts.default_term)
+ self.__g('set style data lines')
+ self.__g('set grid')
+ self.__g('set autoscale')
+ self.__g('set xlabel "'+str(xlabel).encode('ascii','replace')+'"')
+ self.__g('set ylabel "'+str(ylabel).encode('ascii','replace')+'"')
+ #
+ # Tracé du ou des vecteurs demandés
+ indexes = []
+ if step is not None and step in self.__steps:
+ indexes.append(self.__steps.index(step))
+ elif item is not None and item < len(self.__values):
+ indexes.append(item)
+ else:
+ indexes = indexes + range(len(self.__values))
+ #
+ i = -1
+ for index in indexes:
+ self.__g('set title "'+str(title).encode('ascii','replace')+' (pas '+str(index)+')"')
+ if ( type(steps) is type([]) ) or ( type(steps) is type(numpy.array([])) ):
+ Steps = list(steps)
+ else:
+ Steps = range(len(self.__values[index]))
+ #
+ self.__g.plot( self.__gnuplot.Data( Steps, self.__values[index], title=ltitle ) )
+ #
+ if filename != "":
+ i += 1
+ stepfilename = "%s_%03i.ps"%(filename,i)
+ if os.path.isfile(stepfilename):
+ raise ValueError("Error: a file with this name \"%s\" already exists."%stepfilename)
+ self.__g.hardcopy(filename=stepfilename, color=1)
+ if pause:
+ raw_input('Please press return to continue...\n')
+
+ # ---------------------------------------------------------
+ def stepmean(self):
+ """
+ Renvoie la moyenne sur toutes les valeurs sans tenir compte de la
+ longueur des pas. Il faut que le type de base soit compatible avec
+ les types élémentaires numpy.
+ """
+ try:
+ return numpy.matrix(self.__values).mean()
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def stepstd(self, ddof=0):
+ """
+ Renvoie l'écart-type de toutes les valeurs sans tenir compte de la
+ longueur des pas. Il faut que le type de base soit compatible avec
+ les types élémentaires numpy.
+
+ ddof : c'est le nombre de degrés de liberté pour le calcul de
+ l'écart-type, qui est dans le diviseur. Inutile avant Numpy 1.1
+ """
+ try:
+ if numpy.version.version >= '1.1.0':
+ return numpy.matrix(self.__values).std(ddof=ddof)
+ else:
+ return numpy.matrix(self.__values).std()
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def stepsum(self):
+ """
+ Renvoie la somme de toutes les valeurs sans tenir compte de la
+ longueur des pas. Il faut que le type de base soit compatible avec
+ les types élémentaires numpy.
+ """
+ try:
+ return numpy.matrix(self.__values).sum()
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def stepmin(self):
+ """
+ Renvoie le minimum de toutes les valeurs sans tenir compte de la
+ longueur des pas. Il faut que le type de base soit compatible avec
+ les types élémentaires numpy.
+ """
+ try:
+ return numpy.matrix(self.__values).min()
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def stepmax(self):
+ """
+ Renvoie le maximum de toutes les valeurs sans tenir compte de la
+ longueur des pas. Il faut que le type de base soit compatible avec
+ les types élémentaires numpy.
+ """
+ try:
+ return numpy.matrix(self.__values).max()
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ def cumsum(self):
+ """
+ Renvoie la somme cumulée de toutes les valeurs sans tenir compte de la
+ longueur des pas. Il faut que le type de base soit compatible avec
+ les types élémentaires numpy.
+ """
+ try:
+ return numpy.matrix(self.__values).cumsum(axis=0)
+ except:
+ raise TypeError("Base type is incompatible with numpy")
+
+ # On pourrait aussi utiliser les autres attributs d'une "matrix", comme
+ # "tofile", "min"...
+
+ def stepplot(self,
+ steps = None,
+ title = "",
+ xlabel = "",
+ ylabel = "",
+ ltitle = None,
+ geometry = "600x400",
+ filename = "",
+ persist = False,
+ pause = True,
+ ):
+ """
+ Renvoie un affichage unique pour l'ensemble des valeurs à chaque pas, si
+ elles sont compatibles avec un affichage Gnuplot (donc essentiellement
+ un vecteur). Si l'argument "step" existe dans la liste des pas de
+ stockage effectués, renvoie l'affichage de la valeur stockée à ce pas
+ "step". Si l'argument "item" est correct, renvoie l'affichage de la
+ valeur stockée au numéro "item".
+
+ Arguments :
+ - steps : liste unique des pas de l'axe des X, ou None si c'est
+ la numérotation par défaut
+ - title : base du titre général, qui sera automatiquement
+ complétée par la mention du pas
+ - xlabel : label de l'axe des X
+ - ylabel : label de l'axe des Y
+ - ltitle : titre associé au vecteur tracé
+ - geometry : taille en pixels de la fenêtre et position du coin haut
+ gauche, au format X11 : LxH+X+Y (défaut : 600x400)
+ - filename : nom de fichier Postscript pour une sauvegarde,
+ - persist : booléen indiquant que la fenêtre affichée sera
+ conservée lors du passage au dessin suivant
+ Par défaut, persist = False
+ - pause : booléen indiquant une pause après chaque tracé, et
+ attendant un Return
+ Par défaut, pause = True
+ """
+ import os
+ #
+ # Vérification de la disponibilité du module Gnuplot
+ try:
+ import Gnuplot
+ self.__gnuplot = Gnuplot
+ except:
+ raise ImportError("The Gnuplot module is required to plot the object.")
+ #
+ # Vérification et compléments sur les paramètres d'entrée
+ if persist:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -persist -geometry '+geometry
+ else:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -geometry '+geometry
+ if ltitle is None:
+ ltitle = ""
+ if ( type(steps) is type([]) ) or ( type(steps) is type(numpy.array([])) ):
+ Steps = list(steps)
+ else:
+ Steps = range(len(self.__values[0]))
+ self.__g = self.__gnuplot.Gnuplot() # persist=1
+ self.__g('set terminal '+self.__gnuplot.GnuplotOpts.default_term)
+ self.__g('set style data lines')
+ self.__g('set grid')
+ self.__g('set autoscale')
+ self.__g('set title "'+str(title).encode('ascii','replace') +'"')
+ self.__g('set xlabel "'+str(xlabel).encode('ascii','replace')+'"')
+ self.__g('set ylabel "'+str(ylabel).encode('ascii','replace')+'"')
+ #
+ # Tracé du ou des vecteurs demandés
+ indexes = range(len(self.__values))
+ self.__g.plot( self.__gnuplot.Data( Steps, self.__values[indexes.pop(0)], title=ltitle+" (pas 0)" ) )
+ for index in indexes:
+ self.__g.replot( self.__gnuplot.Data( Steps, self.__values[index], title=ltitle+" (pas %i)"%index ) )
+ #
+ if filename != "":
+ self.__g.hardcopy(filename=filename, color=1)
+ if pause:
+ raw_input('Please press return to continue...\n')
+
+# ==============================================================================
+class OneScalar(Persistence):
+ """
+ Classe définissant le stockage d'une valeur unique réelle (float) par pas
+
+ Le type de base peut être changé par la méthode "basetype", mais il faut que
+ le nouveau type de base soit compatible avec les types par éléments de
+ numpy. On peut même utiliser cette classe pour stocker des vecteurs/listes
+ ou des matrices comme dans les classes suivantes, mais c'est déconseillé
+ pour conserver une signification claire des noms.
+ """
+ def __init__(self, name="", unit="", basetype = float):
+ Persistence.__init__(self, name, unit, basetype)
+
+class OneVector(Persistence):
+ """
+ Classe définissant le stockage d'une liste (list) de valeurs homogènes par
+ hypothèse par pas. Pour éviter les confusions, ne pas utiliser la classe
+ "OneVector" pour des données hétérogènes, mais bien "OneList".
+ """
+ def __init__(self, name="", unit="", basetype = list):
+ Persistence.__init__(self, name, unit, basetype)
+
+class OneMatrix(Persistence):
+ """
+ Classe définissant le stockage d'une matrice de valeurs (numpy.matrix) par
+ pas
+ """
+ def __init__(self, name="", unit="", basetype = numpy.matrix):
+ Persistence.__init__(self, name, unit, basetype)
+
+class OneList(Persistence):
+ """
+ Classe définissant le stockage d'une liste de valeurs potentiellement
+ hétérogènes (list) par pas. Pour éviter les confusions, ne pas utiliser la
+ classe "OneVector" pour des données hétérogènes, mais bien "OneList".
+ """
+ def __init__(self, name="", unit="", basetype = list):
+ Persistence.__init__(self, name, unit, basetype)
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ print "======> Un flottant"
+ OBJET_DE_TEST = OneScalar("My float", unit="cm")
+ OBJET_DE_TEST.store( 5.)
+ OBJET_DE_TEST.store(-5.)
+ OBJET_DE_TEST.store( 1.)
+ print "Les pas de stockage :", OBJET_DE_TEST.stepserie()
+ print "Les valeurs :", OBJET_DE_TEST.valueserie()
+ print "La 2ème valeur :", OBJET_DE_TEST.valueserie(1)
+ print "La dernière valeur :", OBJET_DE_TEST.valueserie(-1)
+ print "Valeurs par pas :"
+ print " La moyenne :", OBJET_DE_TEST.mean()
+ print " L'écart-type :", OBJET_DE_TEST.std()
+ print " La somme :", OBJET_DE_TEST.sum()
+ print " Le minimum :", OBJET_DE_TEST.min()
+ print " Le maximum :", OBJET_DE_TEST.max()
+ print "Valeurs globales :"
+ print " La moyenne :", OBJET_DE_TEST.stepmean()
+ print " L'écart-type :", OBJET_DE_TEST.stepstd()
+ print " La somme :", OBJET_DE_TEST.stepsum()
+ print " Le minimum :", OBJET_DE_TEST.stepmin()
+ print " Le maximum :", OBJET_DE_TEST.stepmax()
+ print " La somme cumulée :", OBJET_DE_TEST.cumsum()
+ print "Taille \"shape\" :", OBJET_DE_TEST.shape()
+ print "Taille \"len\" :", len(OBJET_DE_TEST)
+ del OBJET_DE_TEST
+ print
+
+ print "======> Un entier"
+ OBJET_DE_TEST = OneScalar("My int", unit="cm", basetype=int)
+ OBJET_DE_TEST.store( 5 )
+ OBJET_DE_TEST.store(-5 )
+ OBJET_DE_TEST.store( 1.)
+ print "Les pas de stockage :", OBJET_DE_TEST.stepserie()
+ print "Les valeurs :", OBJET_DE_TEST.valueserie()
+ print "La 2ème valeur :", OBJET_DE_TEST.valueserie(1)
+ print "La dernière valeur :", OBJET_DE_TEST.valueserie(-1)
+ print "Valeurs par pas :"
+ print " La moyenne :", OBJET_DE_TEST.mean()
+ print " L'écart-type :", OBJET_DE_TEST.std()
+ print " La somme :", OBJET_DE_TEST.sum()
+ print " Le minimum :", OBJET_DE_TEST.min()
+ print " Le maximum :", OBJET_DE_TEST.max()
+ print "Valeurs globales :"
+ print " La moyenne :", OBJET_DE_TEST.stepmean()
+ print " L'écart-type :", OBJET_DE_TEST.stepstd()
+ print " La somme :", OBJET_DE_TEST.stepsum()
+ print " Le minimum :", OBJET_DE_TEST.stepmin()
+ print " Le maximum :", OBJET_DE_TEST.stepmax()
+ print " La somme cumulée :", OBJET_DE_TEST.cumsum()
+ print "Taille \"shape\" :", OBJET_DE_TEST.shape()
+ print "Taille \"len\" :", len(OBJET_DE_TEST)
+ del OBJET_DE_TEST
+ print
+
+ print "======> Un booléen"
+ OBJET_DE_TEST = OneScalar("My bool", unit="", basetype=bool)
+ OBJET_DE_TEST.store( True )
+ OBJET_DE_TEST.store( False )
+ OBJET_DE_TEST.store( True )
+ print "Les pas de stockage :", OBJET_DE_TEST.stepserie()
+ print "Les valeurs :", OBJET_DE_TEST.valueserie()
+ print "La 2ème valeur :", OBJET_DE_TEST.valueserie(1)
+ print "La dernière valeur :", OBJET_DE_TEST.valueserie(-1)
+ print "Taille \"shape\" :", OBJET_DE_TEST.shape()
+ print "Taille \"len\" :", len(OBJET_DE_TEST)
+ del OBJET_DE_TEST
+ print
+
+ print "======> Un vecteur de flottants"
+ OBJET_DE_TEST = OneVector("My float vector", unit="cm")
+ OBJET_DE_TEST.store( (5 , -5) )
+ OBJET_DE_TEST.store( (-5, 5 ) )
+ OBJET_DE_TEST.store( (1., 1.) )
+ print "Les pas de stockage :", OBJET_DE_TEST.stepserie()
+ print "Les valeurs :", OBJET_DE_TEST.valueserie()
+ print "La 2ème valeur :", OBJET_DE_TEST.valueserie(1)
+ print "La dernière valeur :", OBJET_DE_TEST.valueserie(-1)
+ print "Valeurs par pas :"
+ print " La moyenne :", OBJET_DE_TEST.mean()
+ print " L'écart-type :", OBJET_DE_TEST.std()
+ print " La somme :", OBJET_DE_TEST.sum()
+ print " Le minimum :", OBJET_DE_TEST.min()
+ print " Le maximum :", OBJET_DE_TEST.max()
+ print "Valeurs globales :"
+ print " La moyenne :", OBJET_DE_TEST.stepmean()
+ print " L'écart-type :", OBJET_DE_TEST.stepstd()
+ print " La somme :", OBJET_DE_TEST.stepsum()
+ print " Le minimum :", OBJET_DE_TEST.stepmin()
+ print " Le maximum :", OBJET_DE_TEST.stepmax()
+ print " La somme cumulée :", OBJET_DE_TEST.cumsum()
+ print "Taille \"shape\" :", OBJET_DE_TEST.shape()
+ print "Taille \"len\" :", len(OBJET_DE_TEST)
+ del OBJET_DE_TEST
+ print
+
+ print "======> Une liste hétérogène"
+ OBJET_DE_TEST = OneList("My list", unit="bool/cm")
+ OBJET_DE_TEST.store( (True , -5) )
+ OBJET_DE_TEST.store( (False, 5 ) )
+ OBJET_DE_TEST.store( (True , 1.) )
+ print "Les pas de stockage :", OBJET_DE_TEST.stepserie()
+ print "Les valeurs :", OBJET_DE_TEST.valueserie()
+ print "La 2ème valeur :", OBJET_DE_TEST.valueserie(1)
+ print "La dernière valeur :", OBJET_DE_TEST.valueserie(-1)
+ print "Valeurs par pas : attention, on peut les calculer car True=1, False=0, mais cela n'a pas de sens"
+ print " La moyenne :", OBJET_DE_TEST.mean()
+ print " L'écart-type :", OBJET_DE_TEST.std()
+ print " La somme :", OBJET_DE_TEST.sum()
+ print " Le minimum :", OBJET_DE_TEST.min()
+ print " Le maximum :", OBJET_DE_TEST.max()
+ print "Valeurs globales : attention, on peut les calculer car True=1, False=0, mais cela n'a pas de sens"
+ print " La moyenne :", OBJET_DE_TEST.stepmean()
+ print " L'écart-type :", OBJET_DE_TEST.stepstd()
+ print " La somme :", OBJET_DE_TEST.stepsum()
+ print " Le minimum :", OBJET_DE_TEST.stepmin()
+ print " Le maximum :", OBJET_DE_TEST.stepmax()
+ print " La somme cumulée :", OBJET_DE_TEST.cumsum()
+ print "Taille \"shape\" :", OBJET_DE_TEST.shape()
+ print "Taille \"len\" :", len(OBJET_DE_TEST)
+ del OBJET_DE_TEST
+ print
+
+ print "======> Utilisation directe de la classe Persistence"
+ OBJET_DE_TEST = Persistence("My object", unit="", basetype=int )
+ OBJET_DE_TEST.store( 1 )
+ OBJET_DE_TEST.store( 3 )
+ OBJET_DE_TEST.store( 7 )
+ print "Les pas de stockage :", OBJET_DE_TEST.stepserie()
+ print "Les valeurs :", OBJET_DE_TEST.valueserie()
+ print "La 2ème valeur :", OBJET_DE_TEST.valueserie(1)
+ print "La dernière valeur :", OBJET_DE_TEST.valueserie(-1)
+ print "Taille \"shape\" :", OBJET_DE_TEST.shape()
+ print "Taille \"len\" :", len(OBJET_DE_TEST)
+ del OBJET_DE_TEST
+ print
+
+ print "======> Affichage d'objets stockés"
+ OBJET_DE_TEST = Persistence("My object", unit="", basetype=numpy.array)
+ D = OBJET_DE_TEST
+ vect1 = [1, 2, 1, 2, 1]
+ vect2 = [-3, -3, 0, -3, -3]
+ vect3 = [-1, 1, -5, 1, -1]
+ vect4 = 100*[0.29, 0.97, 0.73, 0.01, 0.20]
+ print "Stockage de 3 vecteurs de longueur identique"
+ D.store(vect1)
+ D.store(vect2)
+ D.store(vect3)
+ print "Affichage de l'ensemble du stockage sur une même image"
+ D.stepplot(
+ title = "Tous les vecteurs",
+ filename="vecteurs.ps",
+ xlabel = "Axe X",
+ ylabel = "Axe Y",
+ pause = False )
+ print "Stockage d'un quatrième vecteur de longueur différente"
+ D.store(vect4)
+ print "Affichage séparé du dernier stockage"
+ D.plot(
+ item = 3,
+ title = "Vecteurs",
+ filename = "vecteur",
+ xlabel = "Axe X",
+ ylabel = "Axe Y",
+ pause = False )
+ print "Les images ont été stockées en fichiers Postscript"
+ print "Taille \"shape\" du dernier objet stocké",OBJET_DE_TEST.shape()
+ print "Taille \"len\" du dernier objet stocké",len(OBJET_DE_TEST)
+ del OBJET_DE_TEST
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Informations sur le code et la plateforme, et mise à jour des chemins
+
+ La classe "PlatformInfo" permet de récupérer les informations générales sur
+ le code et la plateforme sous forme de strings, ou d'afficher directement
+ les informations disponibles par les méthodes. L'impression directe d'un
+ objet de cette classe affiche les informations minimales. Par exemple :
+ print PlatformInfo()
+ print PlatformInfo().getVersion()
+ created = PlatformInfo().getDate()
+
+ La classe "PathManagement" permet de mettre à jour les chemins système pour
+ ajouter les outils numériques, matrices... On l'utilise en instanciant
+ simplement cette classe, sans meme récupérer d'objet :
+ PathManagement()
+"""
+__author__ = "Jean-Philippe ARGAUD - Mars 2008"
+
+import os
+
+# ==============================================================================
+class PlatformInfo:
+ """
+ Rassemblement des informations sur le code et la plateforme
+ """
+ def getName(self):
+ "Retourne le nom de l'application"
+ import version
+ return version.name
+
+ def getVersion(self):
+ "Retourne le numéro de la version"
+ import version
+ return version.version
+
+ def getDate(self):
+ "Retourne la date de création de la version"
+ import version
+ return version.date
+
+ def getPythonVersion(self):
+ "Retourne la version de python utilisée"
+ import sys
+ return ".".join(map(str,sys.version_info[0:3]))
+
+ def getNumpyVersion(self):
+ "Retourne la version de numpy utilisée"
+ import numpy.version
+ return numpy.version.version
+
+ def getScipyVersion(self):
+ "Retourne la version de scipy utilisée"
+ import scipy.version
+ return scipy.version.version
+
+ def getCurrentMemorySize(self):
+ "Retourne la taille mémoire courante utilisée"
+ return 1
+
+ def __str__(self):
+ import version
+ return "%s %s (%s)"%(version.name,version.version,version.date)
+
+# ==============================================================================
+class PathManagement:
+ """
+ Mise à jour du path système pour les répertoires d'outils
+ """
+ def __init__(self):
+ import os, sys
+ parent = os.path.abspath(os.path.join(os.path.dirname(__file__),".."))
+ self.__paths = {}
+ self.__paths["daExternals"] = os.path.join(parent,"daExternals")
+ self.__paths["daMatrices"] = os.path.join(parent,"daMatrices")
+ self.__paths["daNumerics"] = os.path.join(parent,"daNumerics")
+ #
+ for v in self.__paths.values():
+ sys.path.insert(0, v )
+ #
+ # Conserve en unique exemplaire chaque chemin
+ sys.path = list(set(sys.path))
+ del parent
+
+ def getpaths(self):
+ """
+ Renvoie le dictionnaire des chemins ajoutés
+ """
+ return self.__paths
+
+# ==============================================================================
+class SystemUsage:
+ """
+ Permet de récupérer les différentes tailles mémoires du process courant
+ """
+ #
+ # Le module resource renvoie 0 pour les tailles mémoire. On utilise donc
+ # plutôt : http://code.activestate.com/recipes/286222/ et les infos de
+ # http://www.redhat.com/docs/manuals/enterprise/RHEL-4-Manual/en-US/Reference_Guide/s2-proc-meminfo.html
+ #
+ _proc_status = '/proc/%d/status' % os.getpid()
+ _memo_status = '/proc/meminfo'
+ _scale = {
+ 'o': 1.0,
+ 'ko': 1024.0, 'mo': 1024.0*1024.0,
+ 'Ko': 1024.0, 'Mo': 1024.0*1024.0,
+ 'B': 1.0,
+ 'kB': 1024.0, 'mB': 1024.0*1024.0,
+ 'KB': 1024.0, 'MB': 1024.0*1024.0,
+ }
+ _max_mem = 0
+ _max_rss = 0
+ _max_sta = 0
+ #
+ def _VmA(self, VmKey, unit):
+ try:
+ t = open(self._memo_status)
+ v = t.read()
+ t.close()
+ except:
+ return 0.0 # non-Linux?
+ i = v.index(VmKey) # get VmKey line e.g. 'VmRSS: 9999 kB\n ...'
+ v = v[i:].split(None, 3) # whitespace
+ if len(v) < 3:
+ return 0.0 # invalid format?
+ # convert Vm value to bytes
+ mem = float(v[1]) * self._scale[v[2]]
+ return mem / self._scale[unit]
+ #
+ def getAvailablePhysicalMemory(self, unit="o"):
+ "Renvoie la mémoire physique utilisable en octets"
+ return self._VmA('MemTotal:', unit)
+ #
+ def getAvailableSwapMemory(self, unit="o"):
+ "Renvoie la mémoire swap utilisable en octets"
+ return self._VmA('SwapTotal:', unit)
+ #
+ def getAvailableMemory(self, unit="o"):
+ "Renvoie la mémoire totale (physique+swap) utilisable en octets"
+ return self._VmA('MemTotal:', unit) + self._VmA('SwapTotal:', unit)
+ #
+ def getUsableMemory(self, unit="o"):
+ """Renvoie la mémoire utilisable en octets
+ Rq : il n'est pas sûr que ce décompte soit juste...
+ """
+ return self._VmA('MemFree:', unit) + self._VmA('SwapFree:', unit) + \
+ self._VmA('Cached:', unit) + self._VmA('SwapCached:', unit)
+ #
+ def _VmB(self, VmKey, unit):
+ try:
+ t = open(self._proc_status)
+ v = t.read()
+ t.close()
+ except:
+ return 0.0 # non-Linux?
+ i = v.index(VmKey) # get VmKey line e.g. 'VmRSS: 9999 kB\n ...'
+ v = v[i:].split(None, 3) # whitespace
+ if len(v) < 3:
+ return 0.0 # invalid format?
+ # convert Vm value to bytes
+ mem = float(v[1]) * self._scale[v[2]]
+ return mem / self._scale[unit]
+ #
+ def getUsedMemory(self, unit="o"):
+ "Renvoie la mémoire totale utilisée en octets"
+ mem = self._VmB('VmSize:', unit)
+ self._max_mem = max(self._max_mem, mem)
+ return mem
+ #
+ def getUsedResident(self, unit="o"):
+ "Renvoie la mémoire résidente utilisée en octets"
+ mem = self._VmB('VmRSS:', unit)
+ self._max_rss = max(self._max_rss, mem)
+ return mem
+ #
+ def getUsedStacksize(self, unit="o"):
+ "Renvoie la taille du stack utilisé en octets"
+ mem = self._VmB('VmStk:', unit)
+ self._max_sta = max(self._max_sta, mem)
+ return mem
+ #
+ def getMaxUsedMemory(self):
+ "Renvoie la mémoire totale maximale mesurée"
+ return self._max_mem
+ #
+ def getMaxUsedResident(self):
+ "Renvoie la mémoire résidente maximale mesurée"
+ return self._max_rss
+ #
+ def getMaxUsedStacksize(self):
+ "Renvoie la mémoire du stack maximale mesurée"
+ return self._max_sta
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ print PlatformInfo()
+ print
+ p = PlatformInfo()
+ print "Les caractéristiques détaillées des applications et outils sont :"
+ print " - Application.......:",p.getName()
+ print " - Version...........:",p.getVersion()
+ print " - Date Application..:",p.getDate()
+ print " - Python............:",p.getPythonVersion()
+ print " - Numpy.............:",p.getNumpyVersion()
+ print " - Scipy.............:",p.getScipyVersion()
+ print
+
+ p = PathManagement()
+ print "Les chemins ajoutés au système pour des outils :"
+ for k,v in p.getpaths().items():
+ print " %12s : %s"%(k,os.path.basename(v))
+ print
+
+ m = SystemUsage()
+ print "La mémoire disponible est la suivante :"
+ print " - mémoire totale....: %4.1f Mo"%m.getAvailableMemory("Mo")
+ print " - mémoire physique..: %4.1f Mo"%m.getAvailablePhysicalMemory("Mo")
+ print " - mémoire swap......: %4.1f Mo"%m.getAvailableSwapMemory("Mo")
+ print " - utilisable........: %4.1f Mo"%m.getUsableMemory("Mo")
+ print "L'usage mémoire de cette exécution est le suivant :"
+ print " - mémoire totale....: %4.1f Mo"%m.getUsedMemory("Mo")
+ print " - mémoire résidente.: %4.1f Mo"%m.getUsedResident("Mo")
+ print " - taille de stack...: %4.1f Mo"%m.getUsedStacksize("Mo")
+ print "Création d'un objet range(1000000) et mesure mémoire"
+ x = range(1000000)
+ print " - mémoire totale....: %4.1f Mo"%m.getUsedMemory("Mo")
+ print "Destruction de l'objet et mesure mémoire"
+ del x
+ print " - mémoire totale....: %4.1f Mo"%m.getUsedMemory("Mo")
+ print "L'usage mémoire maximal de cette exécution est le suivant :"
+ print " - mémoire totale....: %4.1f Mo"%m.getMaxUsedMemory()
+ print " - mémoire résidente.: %4.1f Mo"%m.getMaxUsedResident()
+ print " - taille de stack...: %4.1f Mo"%m.getMaxUsedStacksize()
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+name = "Data Assimilation Package"
+version = "0.2.0"
+date = "lundi 23 septembre 2009, 11:11:11 (UTC+0200)"
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs
+ dependants au sens du test de Student.
+ Ce diagnostic utilise le calcul de la p-value pour le test de Student
+ pour 2 vecteurs dependants
+ En input : la tolerance
+ En output : le resultat du diagnostic est une reponse booleenne au test :
+ True si les moyennes sont egales au sens du Test de Student
+ False dans le cas contraire.
+"""
+__author__ = "Sophie RICCI - Octobre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from ComputeStudent import DependantVectors
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ """
+ Diagnostic qui effectueIndependantVectorsEqualVariance le test d egalite des moyennes de 2 vecteurs
+ dependants au sens du test de Student.
+ Ce diagnostic utilise le calcul de la p-value pour le test de Student
+ pour 2 vecteurs dependants
+ En input : la tolerance
+ En output : le resultat du diagnostic est une reponse booleenne au test :
+ True si les moyennes sont egales au sens du Test de Student
+ False dans le cas contraire.
+ """
+ def __init__(self, name="", unit="", basetype = None, parameters = {} ):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
+ if not self.parameters.has_key("tolerance"):
+ raise ValueError("A parameter named \"tolerance\" is required.")
+
+ def formula(self, V1, V2):
+ """
+ Effectue le calcul de la p-value de Student pour deux vecteurs.
+ """
+ [aire, Q, reponse, message] = DependantVectors(
+ vector1 = V1,
+ vector2 = V2,
+ tolerance = self.parameters["tolerance"] )
+ logging.info( message )
+ answerStudentTest = False
+ if (aire < (100.*self.parameters["tolerance"])) :
+ answerStudentTest = False
+ else:
+ answerStudentTest = True
+ return answerStudentTest
+
+ def calculate(self, vector1 = None, vector2 = None, step = None):
+ """
+ Active la formule de calcul
+ """
+ if (vector1 is None) or (vector2 is None) :
+ raise ValueError("Two vectors must be given to calculate the Student value")
+ V1 = numpy.array(vector1)
+ V2 = numpy.array(vector2)
+ if (V1.size < 1) or (V2.size < 1):
+ raise ValueError("The given vectors must not be empty")
+ if V1.size != V2.size:
+ raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
+ value = self.formula( V1, V2 )
+ self.store( value = value, step = step)
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ print " Test d'égalite des moyennes au sens de Student pour deux vecteurs"
+ print " dépendants."
+ print
+ #
+ # Initialisation des inputs et appel du diagnostic
+ # --------------------------------------------------------------------
+ tolerance = 0.05
+ D = ElementaryDiagnostic("ComputeMeanStudent_DependVect", parameters = {
+ "tolerance":tolerance,
+ })
+ #
+ # Tirage de l'echantillon aleatoire
+ # --------------------------------------------------------------------
+ x1 = numpy.array(([-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]))
+ x2 = numpy.array(([-0.23, 1.36, 0.32, 0.24, -0.66, -0.19, -0.31, 0.56, 1.21, 0.99]))
+ #
+ # Calcul
+ # --------------------------------------------------------------------
+ D.calculate(x1, x2)
+ #
+ if D.valueserie(0) :
+ print " L'hypothèse d'égalité des moyennes est valide."
+ print
+ else :
+ raise ValueError("The egality of the means is NOT valid")
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs
+ independants supposes de variances differentes au sens du test de Student.
+ En input : la tolerance
+ En output : le resultat du diagnostic est une reponse booleenne au test :
+ True si les moyennes sont egales au sens du Test de Student
+ False dans le cas contraire.
+"""
+__author__ = "Sophie RICCI - Octobre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from ComputeStudent import IndependantVectorsDifferentVariance
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ """
+ Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs
+ independants supposes de variances differentes au sens du test de Student.
+ En input : la tolerance
+ En output : le resultat du diagnostic est une reponse booleenne au test :
+ True si les moyennes sont egales au sens du Test de Student
+ False dans le cas contraire.
+ """
+ def __init__(self, name="", unit="", basetype = None, parameters = {} ):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
+ if not self.parameters.has_key("tolerance"):
+ raise ValueError("A parameter named \"tolerance\" is required.")
+
+ def formula(self, V1, V2):
+ """
+ Effectue le calcul de la p-value de Student pour deux vecteurs
+ independants supposes de variances differentes.
+ """
+ [aire, Q, reponse, message] = IndependantVectorsDifferentVariance(
+ vector1 = V1,
+ vector2 = V2,
+ tolerance = self.parameters["tolerance"],
+ )
+ logging.info( message )
+ answerStudentTest = False
+ if (aire < (100.*self.parameters["tolerance"])) :
+ answerStudentTest = False
+ else:
+ answerStudentTest = True
+ return answerStudentTest
+
+ def calculate(self, vector1 = None, vector2 = None, step = None):
+ """
+ Active la formule de calcul
+ """
+ if (vector1 is None) or (vector2 is None) :
+ raise ValueError("Two vectors must be given to calculate the Student value")
+ V1 = numpy.array(vector1)
+ V2 = numpy.array(vector2)
+ if (V1.size < 1) or (V2.size < 1):
+ raise ValueError("The given vectors must not be empty")
+ if V1.size != V2.size:
+ raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
+ value = self.formula( V1, V2 )
+ self.store( value = value, step = step)
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ print " Test d'égalite des moyennes au sens de Student pour deux vecteurs"
+ print " indépendants supposés de variances différentes."
+ print
+ #
+ # Initialisation des inputs et appel du diagnostic
+ # --------------------------------------------------------------------
+ tolerance = 0.05
+ D = ElementaryDiagnostic("IndependantVectorsDifferentVariance", parameters = {
+ "tolerance":tolerance,
+ })
+ #
+ # Tirage de l'echantillon aleatoire
+ # --------------------------------------------------------------------
+ x1 = numpy.array(([-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]))
+ x2 = numpy.array(([-0.23, 1.36, 0.32, 0.24, -0.66, -0.19, -0.31, 0.56, 1.21, 0.99]))
+ #
+ # Calcul
+ # --------------------------------------------------------------------
+ D.calculate(x1, x2)
+#
+ if D.valueserie(0) :
+ print " L'hypothèse d'égalité des moyennes est valide."
+ print
+ else :
+ raise ValueError("The egality of the means is NOT valid")
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs
+ independants supposes de variances egales au sens du test de Student.
+ En input : la tolerance
+ En output : le resultat du diagnostic est une reponse booleenne au test :
+ True si les moyennes sont egales au sens du Test de Student
+ False dans le cas contraire.
+"""
+__author__ = "Sophie RICCI - Octobre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from ComputeStudent import IndependantVectorsEqualVariance
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ """
+ Diagnostic qui effectue le test d egalite des moyennes de 2 vecteurs independants supposes de variances egales au sens du test de Student.
+ En input : la tolerance
+ En output : le resultat du diagnostic est une reponse booleenne au test :
+ True si les moyennes sont egales au sens du Test de Student
+ False dans le cas contraire.
+ """
+ def __init__(self, name="", unit="", basetype = None, parameters = {} ):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
+ if not self.parameters.has_key("tolerance"):
+ raise ValueError("A parameter named \"tolerance\" is required.")
+
+ def formula(self, V1, V2):
+ """
+ Effectue le calcul de la p-value de Student pour deux vecteurs
+ independants supposes de variances egales.
+ """
+ [aire, Q, reponse, message] = IndependantVectorsEqualVariance(
+ vector1 = V1,
+ vector2 = V2,
+ tolerance = self.parameters["tolerance"],
+ )
+ logging.info( message )
+ answerStudentTest = False
+ if (aire < (100.*self.parameters["tolerance"])) :
+ answerStudentTest = False
+ else:
+ answerStudentTest = True
+ return answerStudentTest
+
+ def calculate(self, vector1 = None, vector2 = None, step = None):
+ """
+ Active la formule de calcul
+ """
+ if (vector1 is None) or (vector2 is None) :
+ raise ValueError("Two vectors must be given to calculate the Student value")
+ V1 = numpy.array(vector1)
+ V2 = numpy.array(vector2)
+ if (V1.size < 1) or (V2.size < 1):
+ raise ValueError("The given vectors must not be empty")
+ if V1.size != V2.size:
+ raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
+ value = self.formula( V1, V2 )
+ self.store( value = value, step = step)
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ print " Test d'égalite des moyennes au sens de Student pour deux vecteurs"
+ print " indépendants supposés de variances égales"
+ print
+ #
+ # Initialisation des inputs et appel du diagnostic
+ # --------------------------------------------------------------------
+ tolerance = 0.05
+ D = ElementaryDiagnostic("ComputeMeanStudent_IndepVect_EgalVar", parameters = {
+ "tolerance":tolerance,
+ })
+ #
+ # Tirage de l'echantillon aleatoire
+ # --------------------------------------------------------------------
+ x1 = numpy.array(([-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]))
+ x2 = numpy.array(([-0.23, 1.36, 0.32, 0.24, -0.66, -0.19, -0.31, 0.56, 1.21, 0.99]))
+ #
+ # Calcul
+ # --------------------------------------------------------------------
+ D.calculate(x1, x2)
+ #
+ if D.valueserie(0) :
+ print " L'hypothèse d'égalité des moyennes est valide."
+ print
+ else :
+ raise ValueError("The egality of the means is NOT valid")
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic qui compare les variances de 2 vecteurs au sens de Fisher à
+ l'aide du calcul de la p-value pour le test de Fisher.
+ - entrée : la tolérance (tolerance) sous forme de paramètres dans le
+ dictionnaire Par, et les deux vecteurs d'échantillons.
+ - sortie : le résultat du diagnostic est une réponse booléenne au test :
+ True si l'égalite des variances est valide au sens du test de Fisher,
+ False dans le cas contraire
+"""
+__author__ = "Sophie RICCI - Juillet 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from ComputeFisher import ComputeFisher
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ """
+ Diagnostic qui compare les variances de 2 vecteurs au sens de Fisher à
+ l'aide du calcul de la p-value pour le test de Fisher.
+ - entrée : la tolérance (tolerance) sous forme de paramètres dans le
+ dictionnaire Par, et les deux vecteurs d'échantillons.
+ - sortie : le résultat du diagnostic est une réponse booléenne au test :
+ True si l'égalite des variances est valide au sens du test de Fisher,
+ False dans le cas contraire
+ """
+ def __init__(self, name="", unit="", basetype = None, parameters = {} ):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
+ if not self.parameters.has_key("tolerance"):
+ raise ValueError("A parameter named \"tolerance\" is required.")
+
+ def formula(self, V1, V2):
+ """
+ Effectue le test de Fisher avec la p-value pour 2 vecteurs
+ """
+ [aire, f, reponse, message] = ComputeFisher(
+ vector1 = V1,
+ vector2 = V2,
+ tolerance = self.parameters["tolerance"],
+ )
+ answerKhisquareTest = False
+ if (aire < (100.*self.parameters["tolerance"])) :
+ answerKhisquareTest = False
+ else:
+ answerKhisquareTest = True
+ logging.info( message )
+ #
+ return answerKhisquareTest
+
+ def calculate(self, vector1 = None, vector2 = None, step = None):
+ """
+ Active la formule de calcul
+ """
+ if (vector1 is None) or (vector2 is None) :
+ raise ValueError("Two vectors must be given to calculate the Fisher p-value")
+ V1 = numpy.array(vector1)
+ V2 = numpy.array(vector2)
+ if (V1.size < 1) or (V2.size < 1):
+ raise ValueError("The given vectors must not be empty")
+ if V1.size != V2.size:
+ raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
+ #
+ value = self.formula( V1, V2 )
+ #
+ self.store( value = value, step = step)
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ print " Test d'égalite des variances pour deux vecteurs de taille 10"
+ print
+ #
+ # Initialisation des inputs et appel du diagnostic
+ # --------------------------------------------------------------------
+ tolerance = 0.05
+ D = ElementaryDiagnostic("CompareVarianceFisher", parameters = {
+ "tolerance":tolerance,
+ })
+ #
+ # Tirage de l'echantillon aleatoire
+ # --------------------------------------------------------------------
+ x1 = numpy.array(([-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]))
+ x2 = numpy.array(([-0.23, 1.36, 0.32, 0.24, -0.66, -0.19, -0.31, 0.56, 1.21, 0.99]))
+ #
+ # Calcul
+ # --------------------------------------------------------------------
+ D.calculate(x1, x2)
+ #
+ if D.valueserie(0) :
+ print " L'hypothèse d'égalité des deux variances est correcte."
+ print
+ else :
+ raise ValueError("L'hypothèse d'égalité des deux variances est fausse.")
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Calcul 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 sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from AssimilationStudy import AssimilationStudy
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = 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
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Calcul de la fonction coût
+"""
+__author__ = "Sophie RICCI - Octobre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from AssimilationStudy import AssimilationStudy
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = float)
+
+ def _formula(self, X, HX, Xb, Y, R, B):
+ """
+ Calcul de la fonction cout
+ """
+ Jb = 1./2. * (X - Xb).T * B.I * (X - Xb)
+ logging.info( "Partial cost function : Jb = %s"%Jb )
+ #
+ Jo = 1./2. * (Y - HX).T * R.I * (Y - HX)
+ logging.info( "Partial cost function : Jo = %s"%Jo )
+ #
+ J = Jb + Jo
+ logging.info( "Total cost function : J = Jo + Jb = %s"%J )
+ return J
+
+ def calculate(self, x = None, 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 = 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
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Calcul de la fonction coût avec Hlin
+ HX = Hxb + Hlin dx
+"""
+__author__ = "Sophie RICCI - Octobre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from AssimilationStudy import AssimilationStudy
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = float)
+ self.__name = str( name )
+
+ def _formula(self, X = None, dX = None, Hlin = None, Xb=None, HXb = None, Y=None, R=None, B=None):
+
+ """
+ Calcul de la fonction cout
+ """
+ HX = HXb + Hlin.T * dX
+ if hasattr(HX, 'A1') :
+ HX = HX.A1
+ #
+ Jb = 1./2. * (X - Xb).T * B.I * (X - Xb)
+ logging.info( "Partial cost function : Jb = %s"%Jb )
+ #
+ Jo = 1./2. * (Y - HX).T * R.I * (Y - HX)
+ logging.info( "Partial cost function : Jo = %s"%Jo )
+ #
+ J = Jb + Jo
+ logging.info( "Total cost function : J = Jo + Jb = %s"%J )
+ return J
+
+ def calculate(self, x = None, dx = None, Hlin = None, xb = None, Hxb = None, yo = None, R = None, B = None , step = None):
+ """
+ Teste les arguments, active la formule de calcul et stocke le résultat
+ """
+ if (x is None) or (xb is None) or (yo is None) or (dx is None):
+ raise ValueError("Vectors x, dx, xb and yo must be given to compute J")
+ dX = dx
+ if hasattr(numpy.matrix(x), 'A1') :
+ X = numpy.matrix(x).A1
+ if hasattr(numpy.matrix(xb), 'A1') :
+ Xb = numpy.matrix(xb).A1
+ if hasattr(numpy.matrix(yo), 'A1') :
+ Y = numpy.matrix(yo).A1
+ B = numpy.matrix(B)
+ R = numpy.matrix(R)
+ if (Hlin is None ) :
+ raise ValueError("HlinT vector must be given")
+ if (Hxb is None ) :
+ raise ValueError("The given vector must be given")
+ HXb = Hxb
+ if (B is None ) or (R is None ):
+ raise ValueError("The matrices B and R must be given")
+ #
+ value = self._formula(X, dX, Hlin, Xb, HXb, Y, R, B)
+ #
+ self.store( value = value, step = step )
+
+#===============================================================================
+if __name__ == "__main__":
+ print "\nAUTOTEST\n"
+ #
+ D = ElementaryDiagnostic("Ma fonction cout")
+ #
+ # Vecteur de type array
+ # ---------------------
+ x = numpy.array([1., 2.])
+ dx = numpy.array([0.1, 0.2])
+ xb = numpy.array([2., 2.])
+ yo = numpy.array([5., 6.])
+ Hlin = numpy.matrix(numpy.identity(2))
+ Hxb = Hlin *xb
+ Hxb = Hxb.T
+ Hxb = Hxb.A1
+ B = numpy.matrix(numpy.identity(2))
+ R = numpy.matrix(numpy.identity(2))
+ #
+ D.calculate( x = x, dx = dx, Hlin = Hlin, xb = xb, Hxb = Hxb, yo = yo, R = R, B = B)
+ print "Le vecteur x choisi est...:", x
+ print "L ebauche xb choisie est...:", xb
+ print "Le vecteur d observation est...:", yo
+ print "B = ", B
+ print "R = ", R
+ print "La fonction cout J vaut ...: %.2e"%D.valueserie(0)
+ #
+ if (abs(D.valueserie(0) - 11.925) > 1.e-6) :
+ raise ValueError("The computation of the cost function is NOT correct")
+ else :
+ print "The computation of the cost function is OK"
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Calcul de la 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 sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from AssimilationStudy import AssimilationStudy
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = 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
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic qui effectue le test 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 sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+from numpy import random
+import Persistence
+from BasicObjects import Diagnostic
+from ComputeKhi2 import ComputeKhi2_Gauss
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ """
+ """
+ def __init__(self, name="", unit="", basetype = None, parameters = {} ):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = 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 = 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 = 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.")
+
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ 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 sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from RMS import ElementaryDiagnostic as RMS
+from AssimilationStudy import AssimilationStudy
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name="", unit="", basetype = None, parameters = {} ):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = 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()
+ 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
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic qui effectue le test 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 sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+from numpy import random
+
+import Persistence
+from BasicObjects import Diagnostic
+from ComputeKhi2 import ComputeKhi2_Homogen
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name="", unit="", basetype = None, parameters = {} ):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = 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 = random.normal(50.,1.5,10000)
+ numpy.random.seed(2490)
+ x2 = 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
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Classe pour tracer simplement un vecteur à chaque pas
+"""
+__author__ = "Jean-Philippe ARGAUD - Juillet 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import os.path
+import numpy
+from BasicObjects import Diagnostic
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ try:
+ import Gnuplot
+ self.__gnuplot = Gnuplot
+ except:
+ raise ImportError("The Gnuplot module is required to plot the vector")
+
+ def _formula(self,
+ Vector, Steps,
+ title, xlabel, ylabel, ltitle,
+ geometry,
+ filename,
+ persist,
+ pause ):
+ """
+ Trace en gnuplot le vecteur Vector, avec une légende générale, en X et
+ en Y
+ """
+ if persist:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -persist -geometry '+geometry
+ else:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -geometry '+geometry
+ #
+ self.__g = self.__gnuplot.Gnuplot() # persist=1
+ self.__g('set terminal '+self.__gnuplot.GnuplotOpts.default_term)
+ self.__g('set style data lines')
+ self.__g('set grid')
+ self.__g('set autoscale')
+ self.__g('set title "'+title +'"')
+ self.__g('set xlabel "'+xlabel+'"')
+ self.__g('set ylabel "'+ylabel+'"')
+ self.__g.plot( self.__gnuplot.Data( Steps, Vector, title=ltitle ) )
+ if filename != "":
+ self.__g.hardcopy(filename=filename, color=1)
+ if pause:
+ raw_input('Please press return to continue...\n')
+ #
+ return 1
+
+ def calculate(self, vector = None, steps = None,
+ title = "", xlabel = "", ylabel = "", ltitle = None,
+ geometry = "600x400",
+ filename = "",
+ persist = False,
+ pause = True ):
+ """
+ Arguments :
+ - vector : le vecteur à tracer, en liste ou en numpy.array
+ - steps : liste unique des pas de l'axe des X, ou None si c'est
+ la numérotation par défaut
+ - title : titre général du dessin
+ - xlabel : label de l'axe des X
+ - ylabel : label de l'axe des Y
+ - ltitle : titre associé au vecteur tracé
+ - geometry : taille en pixels de la fenêtre et position du coin haut
+ gauche, au format X11 : LxH+X+Y (défaut : 600x400)
+ - filename : nom de fichier Postscript pour une sauvegarde à 1 pas
+ Attention, il faut changer le nom à l'appel pour
+ plusieurs pas de sauvegarde
+ - persist : booléen indiquant que la fenêtre affichée sera
+ conservée lors du passage au dessin suivant
+ Par défaut, persist = False
+ - pause : booléen indiquant une pause après chaque tracé, et
+ attendant un Return
+ Par défaut, pause = True
+ """
+ if vector is None:
+ raise ValueError("One vector must be given to plot it.")
+ if ltitle is None:
+ ltitle = ""
+ Vector = numpy.array(vector)
+ if Vector.size < 1:
+ raise ValueError("The given vector must not be empty")
+ if steps is None:
+ Steps = range(len( vector ))
+ elif not ( type(steps) is type([]) or type(steps) is not type(numpy.array([])) ):
+ raise ValueError("The steps must be given as a list/tuple.")
+ else:
+ Steps = list(steps)
+ if os.path.isfile(filename):
+ raise ValueError("Error: a file with this name \"%s\" already exists."%filename)
+ #
+ value = self._formula(
+ Vector = Vector,
+ Steps = Steps,
+ title = str(title).encode('ascii','replace'),
+ xlabel = str(xlabel).encode('ascii','replace'),
+ ylabel = str(ylabel).encode('ascii','replace'),
+ ltitle = str(ltitle),
+ geometry = str(geometry),
+ filename = str(filename),
+ persist = bool(persist),
+ pause = bool(pause) )
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ D = ElementaryDiagnostic("Mon Plot")
+
+ vect = [1, 2, 1, 2, 1]
+ D.calculate(vect, title = "Vecteur 1", xlabel = "Axe X", ylabel = "Axe Y" )
+ vect = [1, 3, 1, 3, 1]
+ D.calculate(vect, title = "Vecteur 2", filename = "vecteur.ps")
+ vect = [1, 1, 1, 1, 1]
+ D.calculate(vect, title = "Vecteur 3")
+ vect = [0.29, 0.97, 0.73, 0.01, 0.20]
+ D.calculate(vect, title = "Vecteur 4")
+ vect = [-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]
+ D.calculate(vect, title = "Vecteur 5")
+ vect = [0.29, 0.97, 0.73, 0.01, 0.20]
+ D.calculate(vect, title = "Vecteur 6 affiche avec une autre geometrie et position", geometry="800x200+50+50")
+ vect = 100*[0.29, 0.97, 0.73, 0.01, 0.20]
+ D.calculate(vect, title = "Vecteur 7 : long construit par repetition")
+ vect = [0.29, 0.97, 0.73, 0.01, 0.20]
+ D.calculate(vect, title = "Vecteur 8", ltitle = "Vecteur 8")
+ temps = [0.1,0.2,0.3,0.4,0.5]
+ D.calculate(vect, temps, title = "Vecteur 8 avec axe du temps modifie")
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Classe pour tracer simplement une liste de vecteurs à chaque pas
+"""
+__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import os.path
+import numpy
+from BasicObjects import Diagnostic
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ try:
+ import Gnuplot
+ self.__gnuplot = Gnuplot
+ except:
+ raise ImportError("The Gnuplot module is required to plot the vector")
+
+ def _formula(self,
+ Vector, Steps,
+ title, xlabel, ylabel, ltitle,
+ geometry,
+ filename,
+ persist,
+ pause ):
+ """
+ Trace en gnuplot chaque vecteur de la liste Vector, avec une légende
+ générale, en X et en Y
+ """
+ if persist:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -persist -geometry '+geometry
+ else:
+ self.__gnuplot.GnuplotOpts.gnuplot_command = 'gnuplot -geometry '+geometry
+ #
+ self.__g = self.__gnuplot.Gnuplot() # persist=1
+ self.__g('set terminal '+self.__gnuplot.GnuplotOpts.default_term)
+ self.__g('set style data lines')
+ self.__g('set grid')
+ self.__g('set autoscale')
+ self.__g('set title "'+title +'"')
+ self.__g('set xlabel "'+xlabel+'"')
+ self.__g('set ylabel "'+ylabel+'"')
+ self.__g.plot( self.__gnuplot.Data( Steps, Vector.pop(0), title=ltitle.pop(0) ) )
+ for vector in Vector:
+ self.__g.replot( self.__gnuplot.Data( Steps, vector, title=ltitle.pop(0) ) )
+ if filename != "":
+ self.__g.hardcopy(filename=filename, color=1)
+ if pause:
+ raw_input('Please press return to continue...\n')
+ #
+ return 1
+
+ def calculate(self, vector = None, steps = None,
+ title = "", xlabel = "", ylabel = "", ltitle = None,
+ geometry = "600x400",
+ filename = "",
+ persist = False,
+ pause = True ):
+ """
+ Arguments :
+ - vector : liste des vecteurs à tracer, chacun étant en liste ou
+ en numpy.array
+ - steps : liste unique des pas, ou None si c'est la numérotation
+ par défaut
+ - title : titre général du dessin
+ - xlabel : label de l'axe des X
+ - ylabel : label de l'axe des Y
+ - ltitle : liste des titres associés à chaque vecteur, dans le
+ même ordre que les vecteurs eux-mêmes
+ - geometry : taille en pixels de la fenêtre et position du coin haut
+ gauche, au format X11 : LxH+X+Y (défaut : 600x400)
+ - filename : nom de fichier Postscript pour une sauvegarde à 1 pas
+ Attention, il faut changer le nom à l'appel pour
+ plusieurs pas de sauvegarde
+ - persist : booléen indiquant que la fenêtre affichée sera
+ conservée lors du passage au dessin suivant
+ Par défaut, persist = False
+ - pause : booléen indiquant une pause après chaque tracé, et
+ attendant un Return
+ Par défaut, pause = True
+ """
+ if vector is None:
+ raise ValueError("One vector must be given to plot it.")
+ if type(vector) is not type([]) and type(vector) is not type(()):
+ raise ValueError("The vector(s) must be given as a list/tuple.")
+ if ltitle is None or len(ltitle) != len(vector):
+ ltitle = ["" for i in range(len(vector))]
+ VectorList = []
+ for onevector in vector:
+ VectorList.append( numpy.array( onevector ) )
+ if VectorList[-1].size < 1:
+ raise ValueError("Each given vector must not be empty.")
+ if steps is None:
+ Steps = range(len(vector[0]))
+ elif not ( type(steps) is type([]) or type(steps) is not type(numpy.array([])) ):
+ raise ValueError("The steps must be given as a list/tuple.")
+ else:
+ Steps = list(steps)
+ if os.path.isfile(filename):
+ raise ValueError("Error: a file with this name \"%s\" already exists."%filename)
+ #
+ value = self._formula(
+ Vector = VectorList,
+ Steps = Steps,
+ title = str(title).encode('ascii','replace'),
+ xlabel = str(xlabel).encode('ascii','replace'),
+ ylabel = str(ylabel).encode('ascii','replace'),
+ ltitle = [str(lt) for lt in ltitle],
+ geometry = str(geometry),
+ filename = str(filename),
+ persist = bool(persist),
+ pause = bool(pause),
+ )
+
+# ==============================================================================
+if __name__ == "__main__":
+ print '\n AUTODIAGNOSTIC \n'
+
+ D = ElementaryDiagnostic("Mon Plot")
+
+ vect1 = [1, 2, 1, 2, 1]
+ D.calculate([vect1,], title = "Vecteur 1", xlabel = "Axe X", ylabel = "Axe Y" )
+ vect2 = [1, 3, 1, 3, 1]
+ D.calculate([vect1,vect2], title = "Vecteurs 1 et 2", filename = "liste_de_vecteurs.ps")
+ vect3 = [-1, 1, -1, 1, -1]
+ D.calculate((vect1,vect2,vect3), title = "Vecteurs 1 a 3")
+ vect4 = 100*[0.29, 0.97, 0.73, 0.01, 0.20]
+ D.calculate([vect4,], title = "Vecteur 4 : long construit par repetition")
+ D.calculate(
+ (vect1,vect2,vect3),
+ [0.1,0.2,0.3,0.4,0.5],
+ title = "Vecteurs 1 a 3, temps modifie",
+ ltitle = ["Vecteur 1","Vecteur 2","Vecteur 3"])
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Calcul d'une RMS
+"""
+__author__ = "Jean-Philippe ARGAUD - Juillet 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import math
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = float)
+
+ def _formula(self, V1, V2):
+ """
+ Fait un écart RMS entre deux vecteurs V1 et V2
+ """
+ rms = math.sqrt( ((V2 - V1)**2).sum() / float(V1.size) )
+ #
+ return rms
+
+ def calculate(self, vector1 = None, vector2 = None, step = None):
+ """
+ Teste les arguments, active la formule de calcul et stocke le résultat
+ """
+ if vector1 is None or vector2 is None:
+ raise ValueError("Two vectors must be given to calculate their RMS")
+ 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'
+
+ D = ElementaryDiagnostic("Ma RMS")
+
+ vect1 = [1, 2, 1, 2, 1]
+ vect2 = [2, 1, 2, 1, 2]
+ D.calculate(vect1,vect2)
+ vect1 = [1, 3, 1, 3, 1]
+ vect2 = [2, 2, 2, 2, 2]
+ D.calculate(vect1,vect2)
+ vect1 = [1, 1, 1, 1, 1]
+ vect2 = [2, 2, 2, 2, 2]
+ D.calculate(vect1,vect2)
+ vect1 = [1, 1, 1, 1, 1]
+ vect2 = [4, -2, 4, -2, -2]
+ D.calculate(vect1,vect2)
+ vect1 = [0.29, 0.97, 0.73, 0.01, 0.20]
+ vect2 = [0.92, 0.86, 0.11, 0.72, 0.54]
+ D.calculate(vect1,vect2)
+ vect1 = [-0.23262176, 1.36065207, 0.32988102, 0.24400551, -0.66765848, -0.19088483, -0.31082575, 0.56849814, 1.21453443, 0.99657516]
+ vect2 = [0,0,0,0,0,0,0,0,0,0]
+ D.calculate(vect1,vect2)
+ print " Les valeurs de RMS attendues sont les suivantes : [1.0, 1.0, 1.0, 3.0, 0.53162016515553656, 0.73784217096601323]"
+ print " Les RMS obtenues................................:", D.valueserie()
+ print " La moyenne......................................:", D.stepmean()
+ print
+
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic sur la reduction du biais lors de l'analyse
+"""
+__author__ = "Sophie RICCI - Aout 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from AssimilationStudy import AssimilationStudy
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = 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
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic sur la reduction de la variance lors de l'analyse
+"""
+__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from AssimilationStudy import AssimilationStudy
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = bool )
+
+ def _formula(self, V1, V2):
+ """
+ Vérification de la reduction de variance sur les écarts entre OMB et OMA
+ lors de l'analyse
+ """
+ varianceOMB = V1.var()
+ varianceOMA = V2.var()
+ #
+ if varianceOMA > varianceOMB:
+ reducevariance = False
+ else :
+ reducevariance = True
+ #
+ return reducevariance
+
+ 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 variance 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 ReduceVariance")
+ #
+ # Vecteur de type matrix
+ # ----------------------
+ 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 " La variance de OMB est de...........................:", x1.var()
+ 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()
+ print " La variance de OMA est de...........................:", x2.var()
+ #
+ D.calculate( vectorOMB = x1, vectorOMA = x2)
+ if not D.valueserie(0) :
+ print " Résultat : l'analyse NE RÉDUIT PAS la variance"
+ else :
+ print " Résultat : l'analyse RÉDUIT la variance"
+ print
+ #
+ # Vecteur de type array
+ # ---------------------
+ x1 = numpy.array(range(11))
+ x2 = numpy.matrix(range(-10,12,2))
+ 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 " La variance de OMB est de...........................:", x1.var()
+ 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()
+ print " La variance de OMA est de...........................:", x2.var()
+ #
+ D.calculate( vectorOMB = x1, vectorOMA = x2)
+ if not D.valueserie(1) :
+ print " Résultat : l'analyse NE RÉDUIT PAS la variance"
+ else :
+ print " Résultat : l'analyse RÉDUIT la variance"
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic sur l'arrêt (ou le ralentissement) de la réduction de la variance
+ au fil des pas (ou itérations) de l'analyse.
+ Ce diagnostic s'applique typiquement au vecteur de différence entre la
+ variance de OMB et la variance de OMA au fil du temps ou des itérations:
+ V[i] = vecteur des VAR(OMB)[i] - VAR(OMA)[i] au temps ou itération i.
+"""
+__author__ = "Sophie Ricci - Septembre 2008"
+
+import sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from AssimilationStudy import AssimilationStudy
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = int )
+
+ def _formula(self, V, CutOffSlope, MultiSlope0):
+ """
+ Recherche du pas de temps ou iteration pour laquelle la reduction
+ de la variance est
+ - inferieure a la valeur seuil CutOffSlope
+ (si une valeure est donnee a CutOffSlope)
+ - inferieure a MultiSlope0 * la pente a la premiere iteration
+ (si une valeure est donnee a MultiSlope0)
+ V[i] = vecteur des VAR(OMB)[i] - VAR(OMA)[i] au temps ou iteration i.
+ """
+ N = V.size
+ pente = numpy.matrix(numpy.zeros((N,))).T
+ iterstopreduction = 0.
+ for i in range (1, N) :
+ pente[i] = V[i]- V[i-1]
+ if pente[i] > 0.0 :
+ raise ValueError("The analysis is INCREASING the variance a l iteration ", i)
+ if CutOffSlope is not None:
+ if numpy.abs(pente[i]) < CutOffSlope :
+ iterstopreduction = i
+ break
+ if MultiSlope0 is not None:
+ if numpy.abs(pente[i]) < MultiSlope0 * numpy.abs(pente[1]) :
+ iterstopreduction = i
+ break
+ #
+ return iterstopreduction
+
+ def calculate(self, vector = None, CutOffSlope = None, MultiSlope0 = None, step = None) :
+ """
+ Teste les arguments, active la formule de calcul et stocke le resultat
+ Arguments :
+ - vector : vecteur des VAR(OMB) - VAR(OMA) au fil des iterations
+ - CutOffSlope : valeur minimale de la pente
+ - MultiSlope0 : Facteur multiplicatif de la pente initiale pour comparaison
+ """
+ if (vector is None) :
+ raise ValueError("One vector must be given to test the convergence of the variance after analysis")
+ V = numpy.array(vector)
+ if V.size < 1 :
+ raise ValueError("The given vector must not be empty")
+ if (MultiSlope0 is None) and (CutOffSlope is None) :
+ raise ValueError("You must set the value of ONE of the CutOffSlope of MultiSlope0 key word")
+ #
+ value = self._formula( V, CutOffSlope, MultiSlope0 )
+ #
+ self.store( value = value, step = step )
+
+#===============================================================================
+if __name__ == "__main__":
+ print "\n AUTODIAGNOSTIC \n"
+
+ # Instanciation de l'objet diagnostic
+ # ------------------------------------------------
+ D = ElementaryDiagnostic("Mon StopReductionVariance")
+
+ # Vecteur de reduction VAR(OMB)-VAR(OMA)
+ # ------------------------------------------------
+ x = numpy.array(([0.60898111, 0.30449056, 0.15224528, 0.07612264, 0.03806132, 0.01903066, 0.00951533, 0.00475766, 0.00237883, 0.00118942]))
+ print " Le vecteur choisi est :", x
+ print " Sur ce vecteur, la reduction a l iteration N = 7 est inferieure a 0.005"
+ print " Sur ce vecteur, la reduction a l iteration N = 8 est inferieure a 0.01 * la reduction a l iteration 1"
+
+ # Comparaison a la valeur seuil de la reduction
+ # ------------------------------------------------
+ D.calculate( vector = x, CutOffSlope = 0.005, MultiSlope0 = None)
+ if (D.valueserie(0) - 7.) < 1.e-15 :
+ print " Test : La comparaison a la valeur seuil de la reduction est juste"
+ else :
+ print " Test : La comparaison a la valeur seuil de la reduction est fausse"
+
+ # Comparaison a alpha* la reduction a la premiere iteration
+ # ------------------------------------------------
+ D.calculate( vector = x, CutOffSlope = None, MultiSlope0 = 0.01)
+ if (D.valueserie(1) - 8.) < 1.e-15 :
+ print " Test : La comparaison a la reduction a la premiere iteration est juste"
+ else :
+ print " Test : La comparaison a la reduction a la premiere iteration est fausse"
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Diagnostic sur 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 sys ; sys.path.insert(0, "../daCore")
+
+import numpy
+import Persistence
+from BasicObjects import Diagnostic
+from scipy.linalg import eig
+import logging
+
+# ==============================================================================
+class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
+ def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
+ Diagnostic.__init__(self, name, parameters)
+ Persistence.OneScalar.__init__( self, name, unit, basetype = 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
+
+
--- /dev/null
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
--- /dev/null
+<?xml version='1.0' encoding='iso-8859-1' ?>
+<proc>
+ <property name="DefaultStudyID" value="1"/>
+ <type name="bool" kind="bool"/>
+ <sequence name="boolvec" content="bool"/>
+ <type name="double" kind="double"/>
+ <sequence name="dblevec" content="double"/>
+ <objref name="file" id="file"/>
+ <type name="int" kind="int"/>
+ <sequence name="intvec" content="int"/>
+ <objref name="pyobj" id="python:obj:1.0"/>
+ <sequence name="seqdblevec" content="dblevec"/>
+ <type name="string" kind="string"/>
+ <sequence name="stringvec" content="string"/>
+ <container name="DefaultContainer">
+ <property name="container_name" value="FactoryServer"/>
+ <property name="hostname" value="localhost"/>
+ </container>
+ <inline name="Building_B">
+ <script><code><![CDATA[# Construction de la matrice de covariances d'erreurs d'ebauche
+import numpy
+
+dimension = len( parametres )
+
+B = numpy.matrix(numpy.core.identity(dimension))
+B = B.A1
+]]></code></script>
+ <inport name="parametres" type="pyobj"/>
+ <outport name="B" type="dblevec"/>
+ </inline>
+ <inline name="Building_R">
+ <script><code><![CDATA[# Construction de la matrice de covariances d'erreurs de mesure
+import numpy
+
+dimension = len( experiences ) * len( experiences[0] )
+
+R = numpy.matrix(numpy.core.identity(dimension))
+R = R.A1
+]]></code></script>
+ <inport name="experiences" type="pyobj"/>
+ <outport name="R" type="dblevec"/>
+ </inline>
+ <inline name="Building_Xb">
+ <script><code><![CDATA[# Script pour extraire l'ebauche comme une liste
+# et une liste de bornes [min,max]
+
+dx = 1.e-2 # Increment en pourcent de Xb
+
+Xb = []
+dXb = []
+Bornes = []
+
+for parametre in parametres:
+ Xb.append( parametre[1] )
+ dXb.append( dx * parametre[1] )
+ Bornes.append( parametre[2:4] )
+]]></code></script>
+ <inport name="parametres" type="pyobj"/>
+ <outport name="Xb" type="dblevec"/>
+ <outport name="Bornes" type="seqdblevec"/>
+ <outport name="dXb" type="dblevec"/>
+ </inline>
+ <inline name="Building_Yo">
+ <script><code><![CDATA[# Script pour extraire l'es mesures/observations et leur nom
+# comme une liste
+
+if len(calcul) != len(experiences):
+ raise ValueError("Les nombres de variables calculees et observees doient etre les memes")
+
+if len(experiences) >=1:
+ nb_observations_par_experience = len( experiences[0] )
+
+Yo = []
+Nom_Yo = []
+
+for i in range( len(experiences) ):
+ Nom_Yo.append( calcul[i][2] )
+ if len( experiences[i] ) != nb_observations_par_experience:
+ raise ValueError("Le nombre de mesures par serie experimentale doit etre constant egal a %i mais la serie %i est longue de %i"%(nb_observations_par_experience,i,len(experiences[i])))
+ for observation in experiences[i]:
+ Yo.append( observation[1] )
+]]></code></script>
+ <inport name="calcul" type="pyobj"/>
+ <inport name="experiences" type="pyobj"/>
+ <outport name="Yo" type="dblevec"/>
+ <outport name="Nom_Yo" type="stringvec"/>
+ </inline>
+ <inline name="Entrees du calcul AD">
+ <script><code><![CDATA[import numpy
+
+print
+print "### ============================================================="
+print "### Verification des arguments fabriques pour passer d'Aster a AD"
+print
+print " Tailles et types :"
+print " de Xb",len(Xb),type(Xb)
+print " de Yo",len(Yo),type(Yo)
+print " de B ",len(B),type(B)
+print " de R ",len(R),type(R)
+print " de H ",len(H),type(H)
+print
+print " Verification de remise a l'echelle :"
+dimensionXb = len( Xb )
+dimensionYo = len( Yo )
+print " de B"
+numpyB = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
+print " de R"
+numpyR = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
+print " de H"
+numpyH = numpy.matrix( H, numpy.float ).reshape((dimensionYo,dimensionXb))
+print "### ============================================================="
+print
+]]></code></script>
+ <inport name="Xb" type="dblevec"/>
+ <inport name="Yo" type="dblevec"/>
+ <inport name="B" type="dblevec"/>
+ <inport name="R" type="dblevec"/>
+ <inport name="H" type="dblevec"/>
+ <outport name="Xb" type="dblevec"/>
+ <outport name="Yo" type="dblevec"/>
+ <outport name="B" type="dblevec"/>
+ <outport name="R" type="dblevec"/>
+ <outport name="H" type="dblevec"/>
+ </inline>
+ <inline name="Sorties du calcul AD">
+ <script><code><![CDATA[import numpy
+print "Diagnostics de sortie de test :"
+print " Remise en numpy des resultats"
+Xa = numpy.array(xa)
+Xb = numpy.array(xb)
+dimensionXb = len( Xb )
+dimensionYo = len( Yo )
+print "Verification de remise a l'echelle de B, R et H"
+B = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
+R = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
+H = numpy.matrix( H, numpy.float ).reshape((dimensionYo,dimensionXb))
+I = numpy.matrix(numpy.core.identity(dimensionYo))
+deltaB = Xa - Xb
+deltaA = (Yo - numpy.dot(H,Xa)).A1
+deltaI = (Yo - numpy.dot(H,Xb)).A1
+print "Calcul de LS"
+print "deltaA",deltaA
+LS = float(numpy.dot(deltaA, deltaA))
+print "Calcul de LSI"
+LSI = float(numpy.dot(deltaI, deltaI))
+print "Calcul de J"
+J = float( numpy.dot(deltaB, numpy.dot(B,deltaB).A1) + numpy.dot(deltaA, numpy.dot(R,deltaA).A1) )
+print "Calcul de JI"
+JI = float( numpy.dot(deltaI, numpy.dot(R,deltaI).A1) )
+
+print
+print "Sortie du test :"
+print " Xb =",Xb
+print " Xa =",Xa
+print "Difference Xa-Xb :"
+print " Xa-Xb =",deltaB
+print " max(Xa-Xb) =",max(deltaB)
+print " min(Xa-Xb) =",min(deltaB)
+print "Fonctionnelles d'ecarts :"
+print " Initialement : J =",JI
+print " LS =",LSI
+print " Analyse : J =",J,"(Baise de %i%s)"%(100.*(JI-J)/J,"%")
+print " LS =",LS,"(Baise de %i%s)"%(100.*(LSI-LS)/LSI,"%")
+print "Autres informations :"
+print " d =",Innovation
+print " A =",A
+print
+]]></code></script>
+ <inport name="xa" type="dblevec"/>
+ <inport name="A" type="dblevec"/>
+ <inport name="Innovation" type="dblevec"/>
+ <inport name="xb" type="dblevec"/>
+ <inport name="Yo" type="dblevec"/>
+ <inport name="B" type="dblevec"/>
+ <inport name="R" type="dblevec"/>
+ <inport name="H" type="dblevec"/>
+ </inline>
+ <parameter>
+ <tonode>Building_B</tonode><toport>parametres</toport>
+ <value><objref>(lp1
+.</objref></value>
+ </parameter>
+ <parameter>
+ <tonode>Building_Xb</tonode><toport>parametres</toport>
+ <value><objref>(lp1
+.</objref></value>
+ </parameter>
+ <parameter>
+ <tonode>Building_Yo</tonode><toport>calcul</toport>
+ <value><objref>(lp1
+.</objref></value>
+ </parameter>
+ <parameter>
+ <tonode>Building_Yo</tonode><toport>experiences</toport>
+ <value><objref>(lp1
+.</objref></value>
+ </parameter>
+ <parameter>
+ <tonode>Sorties du calcul AD</tonode><toport>xa</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Sorties du calcul AD</tonode><toport>A</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Sorties du calcul AD</tonode><toport>Innovation</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Sorties du calcul AD</tonode><toport>xb</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Sorties du calcul AD</tonode><toport>Yo</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Sorties du calcul AD</tonode><toport>B</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Sorties du calcul AD</tonode><toport>R</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Sorties du calcul AD</tonode><toport>H</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Building_R</tonode><toport>experiences</toport>
+ <value><objref>(lp1
+.</objref></value>
+ </parameter>
+ <parameter>
+ <tonode>Entrees du calcul AD</tonode><toport>Xb</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Entrees du calcul AD</tonode><toport>Yo</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Entrees du calcul AD</tonode><toport>B</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Entrees du calcul AD</tonode><toport>R</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>Entrees du calcul AD</tonode><toport>H</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <presentation name="Building_B" x="9" y="81.5" width="227.5" height="110"/>
+ <presentation name="Building_Xb" x="9" y="233" width="227.5" height="168"/>
+ <presentation name="Building_Yo" x="249" y="231" width="227.5" height="139"/>
+ <presentation name="Sorties du calcul AD" x="246" y="432" width="223" height="313"/>
+ <presentation name="Building_R" x="250" y="84.5" width="227.5" height="110"/>
+ <presentation name="Entrees du calcul AD" x="9" y="429" width="227.5" height="226"/>
+ <presentation name="__ROOT__" x="0" y="0" width="487" height="754"/>
+</proc>
--- /dev/null
+<?xml version='1.0' encoding='iso-8859-1' ?>
+<proc>
+ <property name="DefaultStudyID" value="1"/>
+ <type name="bool" kind="bool"/>
+ <sequence name="boolvec" content="bool"/>
+ <type name="double" kind="double"/>
+ <sequence name="dblevec" content="double"/>
+ <objref name="file" id="file"/>
+ <type name="int" kind="int"/>
+ <sequence name="intvec" content="int"/>
+ <objref name="pyobj" id="python:obj:1.0"/>
+ <sequence name="seqdblevec" content="dblevec"/>
+ <sequence name="seqint" content="int"/>
+ <sequence name="seqintvec" content="intvec"/>
+ <type name="string" kind="string"/>
+ <sequence name="stringvec" content="string"/>
+ <container name="DefaultContainer">
+ <property name="container_name" value="FactoryServer"/>
+ <property name="hostname" value="localhost"/>
+ </container>
+ <bloc name="H_linearization">
+ <foreach name="Finite_differences_derivation" type="int">
+ <bloc name="Elementary_calculation">
+ <inline name="Perturbated_point_X">
+ <script><code><![CDATA[print iter
+print seq_X[iter]
+X = seq_X[iter]
+]]></code></script>
+ <inport name="iter" type="int"/>
+ <inport name="seq_X" type="seqdblevec"/>
+ <outport name="X" type="dblevec"/>
+ <outport name="iter" type="int"/>
+ </inline>
+ <inline name="ASTER">
+ <script><code><![CDATA[print 'Debut ASTER_recal'
+
+import os
+execfile( os.path.join(SOURCES_ROOT, 'N_Code_Aster.py') )
+
+print 'RESU_CALC (%s): %s ' % (iter, RESU_CALC)
+print 'DIAG (%s): %s' % (iter, DIAG)
+]]></code></script>
+ <inport name="X" type="dblevec"/>
+ <inport name="iter" type="int"/>
+ <inport name="ASTER_ROOT" type="string"/>
+ <inport name="rcdir" type="string"/>
+ <inport name="debug" type="bool"/>
+ <inport name="DISPLAY" type="string"/>
+ <inport name="SOURCES_ROOT" type="string"/>
+ <inport name="export" type="string"/>
+ <inport name="parametres" type="pyobj"/>
+ <inport name="calcul" type="pyobj"/>
+ <inport name="experience" type="pyobj"/>
+ <inport name="fileparameters" type="string"/>
+ <outport name="FX" type="dblevec"/>
+ <outport name="FY" type="dblevec"/>
+ <outport name="DIMS" type="intvec"/>
+ <outport name="DIAG" type="string"/>
+ <outport name="iter" type="int"/>
+ </inline>
+ <control> <fromnode>Perturbated_point_X</fromnode> <tonode>ASTER</tonode> </control>
+ <datalink control="false">
+ <fromnode>Perturbated_point_X</fromnode> <fromport>X</fromport>
+ <tonode>ASTER</tonode> <toport>X</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Perturbated_point_X</fromnode> <fromport>iter</fromport>
+ <tonode>ASTER</tonode> <toport>iter</toport>
+ </datalink>
+ </bloc>
+ </foreach>
+ <inline name="Gradient">
+ <script><code><![CDATA[
+print "seq_FX=", seq_FX
+print "seq_FY=", seq_FY
+print "seq_DIMS=", seq_DIMS
+print "lst_DIAG=", lst_DIAG
+print "lst_iter=", lst_iter
+print "dX=", dX
+
+import os
+execfile( os.path.join(SOURCES_ROOT, 'N_Gradient.py') )
+
+print "fonctionnelle=", fonctionnelle
+print "gradient=", gradient
+]]></code></script>
+ <inport name="seq_FX" type="seqdblevec"/>
+ <inport name="seq_FY" type="seqdblevec"/>
+ <inport name="seq_DIMS" type="seqintvec"/>
+ <inport name="lst_DIAG" type="stringvec"/>
+ <inport name="lst_iter" type="intvec"/>
+ <inport name="dX" type="dblevec"/>
+ <inport name="SOURCES_ROOT" type="string"/>
+ <outport name="fonctionnelle" type="dblevec"/>
+ <outport name="gradient" type="seqdblevec"/>
+ </inline>
+ <inline name="Input">
+ <function name="inputctrl">
+ <code><![CDATA[# debug : nb permet de limiter le nb de directions pour les differences finies
+nb = None # on calcule pour tous les parametres
+#nb = 1
+
+print
+print "XXXXXXXXXXXXX======== Demarrage de la procedure"
+print
+import copy
+
+def inputctrl( X, dX ):
+
+ if len( X ) != len( dX ):
+ raise ValueError("X and dX has to present the same lenght")
+
+ nbparam = len(X)
+ if nb: nbparam = min(nb, nbparam)
+ seq_X = [ X ]
+ for i in range(nbparam):
+ Xplus = copy.copy(X)
+ Xplus[i] = X[i] + dX[i]
+ seq_X.append( Xplus )
+
+ nb_core = 2
+ itervect = range( len(seq_X) )
+ nbBranches = min(nb_core, len( itervect ))
+
+ print " Taille du vecteur :", nbparam
+ print "seq_X:", seq_X
+ print "nbBranches:", nbBranches
+ print "itervect:", itervect
+
+ # debug
+ #nbBranches = 1
+ #if nbBranches==1:
+ # itervect = [0]
+ # seq_X = [ X ]
+
+ return nbBranches, itervect, seq_X, dX
+]]></code>
+ </function>
+ <inport name="X" type="dblevec"/>
+ <inport name="dX" type="dblevec"/>
+ <outport name="nbBranches" type="int"/>
+ <outport name="itervect" type="intvec"/>
+ <outport name="seq_X" type="seqdblevec"/>
+ <outport name="dX" type="dblevec"/>
+ </inline>
+ <inline name="Temporary_Parameters">
+ <script><code><![CDATA[# Boitier pour une interface simple de variables temporaires]]></code></script>
+ <inport name="ASTER_ROOT" type="string"/>
+ <inport name="rcdir" type="string"/>
+ <inport name="debug" type="bool"/>
+ <inport name="DISPLAY" type="string"/>
+ <inport name="SOURCES_ROOT" type="string"/>
+ <inport name="export" type="string"/>
+ <inport name="parametres" type="pyobj"/>
+ <inport name="calcul" type="pyobj"/>
+ <inport name="experience" type="pyobj"/>
+ <inport name="fileparameters" type="string"/>
+ <outport name="ASTER_ROOT" type="string"/>
+ <outport name="rcdir" type="string"/>
+ <outport name="debug" type="bool"/>
+ <outport name="DISPLAY" type="string"/>
+ <outport name="SOURCES_ROOT" type="string"/>
+ <outport name="export" type="string"/>
+ <outport name="parametres" type="pyobj"/>
+ <outport name="calcul" type="pyobj"/>
+ <outport name="experience" type="pyobj"/>
+ <outport name="fileparameters" type="string"/>
+ </inline>
+ <control> <fromnode>Finite_differences_derivation</fromnode> <tonode>Gradient</tonode> </control>
+ <control> <fromnode>Input</fromnode> <tonode>Finite_differences_derivation</tonode> </control>
+ <control> <fromnode>Input</fromnode> <tonode>Gradient</tonode> </control>
+ <control> <fromnode>Temporary_Parameters</fromnode> <tonode>Finite_differences_derivation</tonode> </control>
+ <control> <fromnode>Temporary_Parameters</fromnode> <tonode>Gradient</tonode> </control>
+ <datalink control="false">
+ <fromnode>Finite_differences_derivation</fromnode> <fromport>SmplPrt</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.Perturbated_point_X</tonode> <toport>iter</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Input</fromnode> <fromport>nbBranches</fromport>
+ <tonode>Finite_differences_derivation</tonode> <toport>nbBranches</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Input</fromnode> <fromport>itervect</fromport>
+ <tonode>Finite_differences_derivation</tonode> <toport>SmplsCollection</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Input</fromnode> <fromport>seq_X</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.Perturbated_point_X</tonode> <toport>seq_X</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Input</fromnode> <fromport>dX</fromport>
+ <tonode>Gradient</tonode> <toport>dX</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>ASTER_ROOT</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>ASTER_ROOT</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>rcdir</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>rcdir</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>debug</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>debug</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>DISPLAY</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>DISPLAY</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>SOURCES_ROOT</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>SOURCES_ROOT</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>SOURCES_ROOT</fromport>
+ <tonode>Gradient</tonode> <toport>SOURCES_ROOT</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>export</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>export</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>parametres</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>parametres</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>calcul</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>calcul</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>experience</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>experience</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Temporary_Parameters</fromnode> <fromport>fileparameters</fromport>
+ <tonode>Finite_differences_derivation.Elementary_calculation.ASTER</tonode> <toport>fileparameters</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>FX</fromport>
+ <tonode>Gradient</tonode> <toport>seq_FX</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>FY</fromport>
+ <tonode>Gradient</tonode> <toport>seq_FY</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>DIMS</fromport>
+ <tonode>Gradient</tonode> <toport>seq_DIMS</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>DIAG</fromport>
+ <tonode>Gradient</tonode> <toport>lst_DIAG</toport>
+ </datalink>
+ <datalink control="false">
+ <fromnode>Finite_differences_derivation.Elementary_calculation.ASTER</fromnode> <fromport>iter</fromport>
+ <tonode>Gradient</tonode> <toport>lst_iter</toport>
+ </datalink>
+ </bloc>
+ <parameter>
+ <tonode>H_linearization.Finite_differences_derivation.Elementary_calculation.ASTER</tonode><toport>X</toport>
+ <value><array><data>
+<value><double>80000</double></value>
+<value><double>1000</double></value>
+<value><double>30</double></value>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>ASTER_ROOT</toport>
+ <value><string>''</string></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>rcdir</toport>
+ <value><string>''</string></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>debug</toport>
+ <value><boolean>0</boolean></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>DISPLAY</toport>
+ <value><string>:0.0</string></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>SOURCES_ROOT</toport>
+ <value><string>.</string></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>export</toport>
+ <value><string>''</string></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>parametres</toport>
+ <value><objref><![CDATA[(lp1
+.]]></objref></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>calcul</toport>
+ <value><objref><![CDATA[(lp1
+.]]></objref></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>experience</toport>
+ <value><objref><![CDATA[(lp1
+.]]></objref></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Temporary_Parameters</tonode><toport>fileparameters</toport>
+ <value><string>[]</string></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Input</tonode><toport>X</toport>
+ <value><array><data>
+<value><double>80000</double></value>
+<value><double>1000</double></value>
+<value><double>30</double></value>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>H_linearization.Input</tonode><toport>dX</toport>
+ <value><array><data>
+<value><double>0.001</double></value>
+<value><double>0.001</double></value>
+<value><double>0.0001</double></value>
+</data></array></value>
+ </parameter>
+ <presentation name="H_linearization.Finite_differences_derivation.Elementary_calculation.ASTER" x="396" y="80.5" width="227.5" height="429"/>
+ <presentation name="H_linearization.Finite_differences_derivation" x="368.856" y="80.5" width="651.5" height="657.5"/>
+ <presentation name="H_linearization.Finite_differences_derivation.Elementary_calculation.Perturbated_point_X" x="9" y="81" width="227.5" height="139"/>
+ <presentation name="H_linearization.Gradient" x="1064.38" y="338.546" width="227.5" height="284"/>
+ <presentation name="H_linearization.Finite_differences_derivation.Elementary_calculation" x="9" y="130" width="633" height="518.5"/>
+ <presentation name="H_linearization.Temporary_Parameters" x="10.3499" y="364.775" width="258.5" height="371"/>
+ <presentation name="H_linearization" x="10.06" y="80.75" width="1301.38" height="747"/>
+ <presentation name="H_linearization.Input" x="9" y="83.5" width="227.5" height="197"/>
+ <presentation name="__ROOT__" x="0" y="0" width="1320.44" height="836.75"/>
+</proc>
--- /dev/null
+<?xml version='1.0' encoding='iso-8859-1' ?>
+<proc>
+ <property name="DefaultStudyID" value="1"/>
+ <type name="bool" kind="bool"/>
+ <sequence name="boolvec" content="bool"/>
+ <type name="double" kind="double"/>
+ <sequence name="dblevec" content="double"/>
+ <objref name="file" id="file"/>
+ <type name="int" kind="int"/>
+ <sequence name="intvec" content="int"/>
+ <objref name="pyobj" id="python:obj:1.0"/>
+ <type name="string" kind="string"/>
+ <sequence name="stringvec" content="string"/>
+ <container name="DefaultContainer">
+ <property name="container_name" value="FactoryServer"/>
+ <property name="hostname" value="localhost"/>
+ </container>
+ <inline name="BLUE par matrices">
+ <function name="algorithm">
+ <code><![CDATA[import sys, os
+sys.path.insert(0, "../../Sources/daCore")
+sys.path.insert(0, "../../ComposantAD/daCore")
+#sys.path.insert(0, os.path.join(os.environ["HOME"],"SALOME5/supplements_JPA/ComposantAD/daCore"))
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+def algorithm(Yo, B, R, H, Xb):
+ #
+ # Remise en place des matrices
+ # -------------------
+ dimensionXb = len( Xb )
+ dimensionYo = len( Yo )
+ B = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
+ R = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
+ H = numpy.matrix( H, numpy.float ).reshape((dimensionYo,dimensionXb))
+ #
+ # Analyse
+ # -------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = Xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = Yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setAlgorithm(choice="Blue")
+ #
+ ADD.analyze()
+ #
+ Xa = ADD.get("Analysis").valueserie(0)
+ Innovation = ADD.get("Innovation").valueserie(0)
+ A = []
+ #
+ return Xa, A, Innovation
+]]></code>
+ </function>
+ <inport name="Yo" type="dblevec"/>
+ <inport name="B" type="dblevec"/>
+ <inport name="R" type="dblevec"/>
+ <inport name="H" type="dblevec"/>
+ <inport name="Xb" type="dblevec"/>
+ <outport name="Xa" type="dblevec"/>
+ <outport name="A" type="dblevec"/>
+ <outport name="Innovation" type="dblevec"/>
+ </inline>
+ <inline name="3D-VAR par matrices">
+ <function name="algorithm">
+ <code><![CDATA[import sys, os
+sys.path.insert(0, "../../Sources/daCore")
+sys.path.insert(0, "../../ComposantAD/daCore")
+#sys.path.insert(0, os.path.join(os.environ["HOME"],"SALOME5/supplements_JPA/ComposantAD/daCore"))
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+def algorithm(Yo, B, R, H, Xb):
+ #
+ # Remise en place des matrices
+ # -------------------
+ dimensionXb = len( Xb )
+ dimensionYo = len( Yo )
+ B = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
+ R = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
+ H = numpy.matrix( H, numpy.float ).reshape((dimensionYo,dimensionXb))
+ #
+ # Analyse
+ # -------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = Xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = Yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asMatrix = H )
+ #
+ ADD.setAlgorithm(choice="3DVAR")
+ #
+ ADD.analyze()
+ #
+ Xa = ADD.get("Analysis").valueserie(0)
+ Innovation = ADD.get("Innovation").valueserie(0)
+ A = []
+ #
+ return Xa, A, Innovation
+]]></code>
+ </function>
+ <inport name="Yo" type="dblevec"/>
+ <inport name="B" type="dblevec"/>
+ <inport name="R" type="dblevec"/>
+ <inport name="H" type="dblevec"/>
+ <inport name="Xb" type="dblevec"/>
+ <outport name="Xa" type="dblevec"/>
+ <outport name="A" type="dblevec"/>
+ <outport name="Innovation" type="dblevec"/>
+ </inline>
+ <inline name="3D-VAR par fonctions">
+ <function name="algorithm">
+ <code><![CDATA[import sys, os
+sys.path.insert(0, "../../Sources/daCore")
+sys.path.insert(0, "../../ComposantAD/daCore")
+#sys.path.insert(0, os.path.join(os.environ["HOME"],"SALOME5/supplements_JPA/ComposantAD/daCore"))
+
+import numpy
+from AssimilationStudy import AssimilationStudy
+
+def algorithm( Yo, B, R, FunctionH, TangentH, AdjointH, Xb, Bounds ):
+ #
+ # Remise en place des matrices
+ # -------------------
+ dimensionXb = len( Xb )
+ dimensionYo = len( Yo )
+ B = numpy.matrix( B, numpy.float ).reshape((dimensionXb,dimensionXb))
+ R = numpy.matrix( R, numpy.float ).reshape((dimensionYo,dimensionYo))
+ #
+ # Analyse
+ # -------
+ ADD = AssimilationStudy()
+ ADD.setBackground (asVector = Xb )
+ ADD.setBackgroundError (asCovariance = B )
+ ADD.setObservation (asVector = Yo )
+ ADD.setObservationError (asCovariance = R )
+ ADD.setObservationOperator(asFunction = {"Direct":FunctionH,
+ "Tangent":TangentH,
+ "Adjoint":AdjointH} )
+ #
+ ADD.setAlgorithm(choice="3DVAR")
+ ADD.setAlgorithmParameters(asDico={
+ "Minimizer":"LBFGSB",
+ "Bounds" :Bounds,
+ })
+ #
+ ADD.analyze()
+ #
+ Xa = ADD.get("Analysis").valueserie(0)
+ Innovation = ADD.get("Innovation").valueserie(0)
+ A = []
+ #
+ return Xa, A, Innovation
+]]></code>
+ </function>
+ <inport name="Yo" type="dblevec"/>
+ <inport name="B" type="dblevec"/>
+ <inport name="R" type="dblevec"/>
+ <inport name="FunctionH" type="pyobj"/>
+ <inport name="TangentH" type="pyobj"/>
+ <inport name="AdjointH" type="pyobj"/>
+ <inport name="Xb" type="dblevec"/>
+ <inport name="Bounds" type="pyobj"/>
+ <outport name="Xa" type="dblevec"/>
+ <outport name="A" type="dblevec"/>
+ <outport name="Innovation" type="dblevec"/>
+ </inline>
+ <parameter>
+ <tonode>BLUE par matrices</tonode><toport>Yo</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>BLUE par matrices</tonode><toport>B</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>BLUE par matrices</tonode><toport>R</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>BLUE par matrices</tonode><toport>H</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>BLUE par matrices</tonode><toport>Xb</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par matrices</tonode><toport>Yo</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par matrices</tonode><toport>B</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par matrices</tonode><toport>R</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par matrices</tonode><toport>H</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par matrices</tonode><toport>Xb</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par fonctions</tonode><toport>Yo</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par fonctions</tonode><toport>B</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par fonctions</tonode><toport>R</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par fonctions</tonode><toport>Xb</toport>
+ <value><array><data>
+</data></array></value>
+ </parameter>
+ <parameter>
+ <tonode>3D-VAR par fonctions</tonode><toport>Bounds</toport>
+ <value><objref>(lp1
+.</objref></value>
+ </parameter>
+ <presentation name="BLUE par matrices" x="9" y="80.5" width="227.5" height="226"/>
+ <presentation name="3D-VAR par matrices" x="237.5" y="80.5" width="227.5" height="226"/>
+ <presentation name="3D-VAR par fonctions" x="465.5" y="80.5" width="227.5" height="313"/>
+ <presentation name="__ROOT__" x="0" y="0" width="702" height="402.5"/>
+</proc>
--- /dev/null
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
--- /dev/null
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Outil numérique de calcul de la variable de Fisher pour comparer les
+ variances de 2 échantillons
+
+ Ce calcul nécessite :
+ - en input :
+ - les deux vecteurs (comme liste, array ou matrix) d'échantillons
+ dont on veut comparer la variance,
+ - la tolérance
+ - en output :
+ - la p-value,
+ - la valeur de la variable aléatoire,
+ - la réponse au test ainsi que
+ - le message qui interprete la reponse du test.
+"""
+__author__ = "Sophie RICCI - Juillet 2008"
+
+import numpy
+from scipy.stats import betai
+
+# ==============================================================================
+def ComputeFisher(vector1 = None, vector2 = None, tolerance = 0.05 ):
+ """
+ Outil numérique de calcul de la variable de Fisher pour comparer les
+ variances de 2 échantillons
+
+ Ce calcul nécessite :
+ - en input : les deux vecteurs (comme liste, array ou matrix)
+ d'échantillons dont on veut comparer la variance, la
+ tolérance
+ - en output : la p-value, la valeur de la variable aléatoire,
+ la réponse au test ainsi que le message qui interprete
+ la reponse du test.
+ """
+ if (vector1 is None) or (vector2 is None) :
+ raise ValueError("Two vectors must be given to calculate the Fisher value value")
+ V1 = numpy.array(vector1)
+ V2 = numpy.array(vector2)
+ if (V1.size < 1) or (V2.size < 1):
+ raise ValueError("The given vectors must not be empty")
+ #
+ # Calcul des variances des echantillons
+ # -------------------------------------
+ # où var est calculee comme : var = somme (xi -xmean)**2 /(n-1)
+ n1 = V1.size
+ n2 = V2.size
+ var1 = V1.std() * V1.std()
+ var2 = V2.std() * V2.std()
+ if (var1 > var2):
+ f = var1/var2
+ df1 = n1-1
+ df2 = n2-1
+ else:
+ f= var2/var1
+ df1 = n2-1
+ df2 = n1-1
+ prob1= betai(0.5*df2,0.5*df1,float(df2)/float(df2+df1*f))
+ prob2= (1. - betai(0.5*df1, 0.5*df2, float(df1)/float(df1+df2/f)))
+ prob = prob1 + prob2
+ #
+ # Calcul de la p-value
+ # --------------------
+ areafisher = 100 * prob
+ #
+ # Test
+ # ----
+ message = "Il y a %.2f%s de chance de se tromper en refusant l'hypothèse d'égalité des variances des 2 échantillons (si <%.2f%s, on refuse effectivement l'égalité)"%(areafisher,"%",100.*tolerance,"%")
+ if (areafisher < (100.*tolerance)) :
+ answerTestFisher = False
+ else:
+ answerTestFisher = True
+ # print "La reponse au test est", answerTestFisher
+
+ return areafisher, f, answerTestFisher, message
+
+# ==============================================================================
+if __name__ == "__main__":
+ print "\nAUTOTEST\n"
+ #
+ # Echantillons
+ # ------------
+ x1 = [-1., 0., 4., 2., -1., 3.]
+ x2 = [-1., 0., 4., 2., -1., 3.]
+ #
+ # Appel du calcul
+ # ---------------
+ [aire, f, reponse, message] = ComputeFisher(
+ vector1 = x1,
+ vector2 = x2,
+ tolerance = 0.05 )
+ #
+ print " aire.....:", aire
+ print " f........:", f
+ print " reponse..:", reponse
+ print " message..:", message
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Outil 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
+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 ComputeArea(self):
+ """
+ Calcul de la p-value
+ """
+ self.__areakhi2 = 100 * chisqprob(self.__Q, self.__ddl)
+ return self.__areakhi2
+
+ 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 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 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 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
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+__doc__ = """
+ Outil numérique de calcul 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
--- /dev/null
+#
+# Copyright (C) 2008-2009 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
