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