src/daComposant/daDiagnostics/HLinearity.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2010  EDF R&D
   4 #
   5 #  This library is free software; you can redistribute it and/or
   6 #  modify it under the terms of the GNU Lesser General Public
   7 #  License as published by the Free Software Foundation; either
   8 #  version 2.1 of the License.
   9 #
  10 #  This library is distributed in the hope that it will be useful,
  11 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 #  Lesser General Public License for more details.
  14 #
  15 #  You should have received a copy of the GNU Lesser General Public
  16 #  License along with this library; if not, write to the Free Software
  17 #  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 #  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 __doc__ = """
  22     Diagnotic de test sur la validité de l'hypothèse de linéarité de l'opérateur
  23     H entre xp et xm
  24
  25     Pour calculer Hlin on utilise un schéma différences finies centrées 2
  26     points. On définit un dxparam tel que :
  27         xp = xb + dxparam
  28     et
  29         xm = xb - dxparam
  30     On calcule Hxp et Hxm pour obtenir Hlin. Hlin est utilise dans le Blue  pour
  31     caler un paramêtre. La question importante est de choisir un dxparam pas
  32     trop grand.
  33
  34     On veut vérifier ici que l'hypothèse de linéarite du modèle par rapport  au
  35     paramêtre est valide sur l'intervalle du paramêtre [xm, xp]. Pour cela on
  36     s'assure que l'on peut retrouver la valeur Hxb par les développemenents de
  37     Taylor en xp et xm. Ainsi on calcule 2 estimations de Hxb, l'une à partir de
  38     Hxp (notee Hx1) et l'autre à partir de Hxm (notee Hx2), que l'on compare à
  39     la valeur calculée de Hxb. On s'intèresse ensuite a la distance entre Hxb et
  40     ses estimés Hx1 et Hx2. Si la distance est inférieure a un seuil de
  41     tolerance, l hypothese est valide.
  42 """
  43 __author__ = "Sophie RICCI - Septembre 2008"
  44
  45 import numpy
  46 from daCore import AssimilationStudy, BasicObjects, Persistence
  47
  48 # ==============================================================================
  49 class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
  50     def __init__(self, name="", unit="", basetype = None, parameters = {} ):
  51         BasicObjects.Diagnostic.__init__(self, name, parameters)
  52         Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
  53         if not self.parameters.has_key("tolerance"):
  54             raise ValueError("A parameter named \"tolerance\" is required.")
  55
  56     def formula(self, H, dxparam,  Hxp, Hxm, Hx):
  57         """
  58         Test sur la validite de l hypothese de linearite de H entre xp et xm
  59         """
  60         dimension = numpy.size(Hx)
  61         #
  62         # Reconstruit les valeurs Hx1 et Hx2 de Hx a partir de Hxm et Hxp
  63         # ---------------------------------------------------------------
  64         Hx1 = Hxm + H.T * dxparam
  65         Hx2 = Hxp - H.T * dxparam
  66         #
  67         # Calcul de l'ecart entre Hx1 et Hx et entre Hx2 et Hx
  68         # ----------------------------------------------------
  69         ADD = AssimilationStudy.AssimilationStudy()
  70         ADD.setDiagnostic("RMS",
  71             name = "Calcul de la RMS entre Hx1 et Hx et entre Hx2 et Hx")
  72         RMS = ADD.get("Calcul de la RMS entre Hx1 et Hx et entre Hx2 et Hx")
  73         RMS.calculate(Hx1,Hx)
  74         std1 = RMS.valueserie(0)
  75         RMS.calculate(Hx2,Hx)
  76         std2 = RMS.valueserie(1)
  77         #
  78         # Normalisation des écarts par Hx pour comparer a un pourcentage
  79         # --------------------------------------------------------------
  80         RMS.calculate(Hx,Hx-Hx)
  81         std = RMS.valueserie(2)
  82         err1=std1/std
  83         err2=std2/std
  84         #
  85         # Comparaison
  86         # -----------
  87         if ( (err1 < self.parameters["tolerance"]) and (err2 < self.parameters["tolerance"]) ):
  88             reponse = True
  89         else:
  90             reponse = False
  91         return reponse
  92
  93     def calculate(self, Hlin = None, deltaparam = None, Hxp = None,   Hxm = None, Hx = None, step = None):
  94         """
  95         Arguments :
  96             - Hlin : Operateur d obsevation lineaire
  97             - deltaparam : pas sur le parametre param
  98             - Hxp : calcul en xp = xb + deltaparam
  99             - Hxm : calcul en xm = xb - deltaparam
 100             - Hx : calcul en x (generalement xb)
 101         """
 102         value = self.formula(  Hlin,  deltaparam, Hxp, Hxm, Hx )
 103         #
 104         self.store( value = value,  step = step)
 105
 106 #===============================================================================
 107 if __name__ == "__main__":
 108     print '\n AUTODIAGNOSTIC \n'
 109
 110     print " Diagnotic de test sur la validité de l'hypothèse de linéarité de"
 111     print " l'opérateur H entre xp et xm."
 112     print
 113     #
 114     dimension = 3
 115     #
 116     # Définition des données
 117     # ----------------------
 118     Hx = numpy.array(([ 2., 4., 6.]))
 119     Hxp = numpy.array(([ 3., 5., 7.]))
 120     Hxm = numpy.array(([ 1., 3., 5.]))
 121     H =  (Hxp - Hxm)/(2.)
 122     dxparam = 1.
 123     #
 124     # Instanciation de l'objet diagnostic
 125     # -----------------------------------
 126     D = ElementaryDiagnostic("Mon TestHlin", parameters = {"tolerance": 0.1})
 127     #
 128     # Calcul
 129     # ------
 130     D.calculate( Hlin = H, deltaparam = dxparam, Hxp = Hxp, Hxm = Hxm, Hx = Hx)
 131
 132     # Validation du calcul
 133     # --------------------
 134     if not D.valueserie(0) :
 135         raise ValueError("La linearisation de H autour de x entre xm et xp est fausse pour ce cas test lineaire")
 136     else :
 137         print " La linéarisation de H autour de x entre xm et xp est valide pour ce cas-test linéaire."
 138         print