src/daComposant/daDiagnostics/ComputeCostFunction.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2009  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     Calcul de la fonction coût
  23 """
  24 __author__ = "Sophie RICCI - Octobre 2008"
  25
  26 import sys ; sys.path.insert(0, "../daCore")
  27
  28 import numpy
  29 import Persistence
  30 from BasicObjects import Diagnostic
  31 from AssimilationStudy import AssimilationStudy
  32 import logging
  33
  34 # ==============================================================================
  35 class ElementaryDiagnostic(Diagnostic,Persistence.OneScalar):
  36     def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
  37         Diagnostic.__init__(self, name)
  38         Persistence.OneScalar.__init__( self, name, unit, basetype = float)
  39
  40     def _formula(self, X, HX, Xb, Y, R, B):
  41         """
  42         Calcul de la fonction cout
  43         """
  44         Jb = 1./2. * (X - Xb).T * B.I * (X - Xb)
  45         logging.info( "Partial cost function : Jb = %s"%Jb )
  46         #
  47         Jo = 1./2. * (Y - HX).T  * R.I * (Y - HX)
  48         logging.info( "Partial cost function : Jo = %s"%Jo )
  49         #
  50         J = Jb + Jo
  51         logging.info( "Total cost function : J = Jo + Jb = %s"%J )
  52         return J
  53
  54     def calculate(self, x = None, Hx = None, xb = None, yo = None, R = None, B = None , step = None):
  55         """
  56         Teste les arguments, active la formule de calcul et stocke le résultat
  57         """
  58         if (x is None) or (xb is None) or (yo is None) :
  59             raise ValueError("Vectors x, xb and yo must be given to compute J")
  60 #        if (type(x) is not float) and (type(x) is not numpy.float64) :
  61 #            if (x.size < 1) or (xb.size < 1) or (yo.size < 1):
  62 #                raise ValueError("Vectors x, xb and yo must not be empty")
  63         if hasattr(numpy.matrix(x),'A1') :
  64             X = numpy.matrix(x).A1
  65         if hasattr(numpy.matrix(xb),'A1') :
  66             Xb = numpy.matrix(xb).A1
  67         if hasattr(numpy.matrix(yo),'A1') :
  68             Y = numpy.matrix(yo).A1
  69         B = numpy.matrix(B)
  70         R = numpy.matrix(R)
  71         if (Hx is None ) :
  72             raise ValueError("The given vector must be given")
  73 #        if (Hx.size < 1) :
  74 #            raise ValueError("The given vector must not be empty")
  75         HX = Hx.A1
  76         if (B is None ) or (R is None ):
  77             raise ValueError("The matrices B and R must be given")
  78 #        if (B.size < 1) or (R.size < 1) :
  79 #            raise ValueError("The matrices B and R must not be empty")
  80         #
  81         value = self._formula(X, HX, Xb, Y, R, B)
  82         #
  83         self.store( value = value,  step = step )
  84
  85 #===============================================================================
  86 if __name__ == "__main__":
  87     print "\nAUTOTEST\n"
  88     #
  89     D = ElementaryDiagnostic("Ma fonction cout")
  90     #
  91     # Vecteur de type array
  92     # ---------------------
  93     x = numpy.array([1., 2.])
  94     xb = numpy.array([2., 2.])
  95     yo = numpy.array([5., 6.])
  96     H = numpy.matrix(numpy.identity(2))
  97     Hx = H*x
  98     Hx = Hx.T
  99     B =  numpy.matrix(numpy.identity(2))
 100     R =  numpy.matrix(numpy.identity(2))
 101     #
 102     D.calculate( x = x, Hx = Hx, xb = xb, yo = yo, R = R, B = B)
 103     print "Le vecteur x choisi est...:", x
 104     print "L ebauche xb choisie est...:", xb
 105     print "Le vecteur d observation est...:", yo
 106     print "B = ", B
 107     print "R = ", R
 108     print "La fonction cout J vaut ...: %.2e"%D.valueserie(0)
 109     print "La fonction cout J vaut ...: ",D.valueserie(0)
 110
 111     if (abs(D.valueserie(0) - 16.5) > 1.e-6)  :
 112         raise ValueError("The computation of the cost function is NOT correct")
 113     else :
 114         print "The computation of the cost function is OK"
 115     print
 116     #
 117     # float simple
 118     # ------------
 119     x = 1.
 120     print type(x)
 121     xb = 2.
 122     yo = 5.
 123     H = numpy.matrix(numpy.identity(1))
 124     Hx = numpy.dot(H,x)
 125     Hx = Hx.T
 126     B =  1.
 127     R =  1.
 128     #
 129     D.calculate( x = x, Hx = Hx, xb = xb, yo = yo, R = R, B = B)
 130     print "Le vecteur x choisi est...:", x
 131     print "L ebauche xb choisie est...:", xb
 132     print "Le vecteur d observation est...:", yo
 133     print "B = ", B
 134     print "R = ", R
 135     print "La fonction cout J vaut ...: %.2e"%D.valueserie(1)
 136     if (abs(D.valueserie(1) - 8.5) > 1.e-6)  :
 137         raise ValueError("The computation of the cost function is NOT correct")
 138     else :
 139         print "The computation of the cost function is OK"
 140     print
 141