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