| # Gaudi Framework, version v22r0 |

Home |  | Generated: 9 Feb 2011

# FuncMinimum.cpp

Go to the documentation of this file.

    00001 // $Id: FuncMinimum.cpp,v 1.4 2006/01/10 20:00:05 hmd Exp $
    00002 // ============================================================================
    00003 // Include files
    00004
    00005 #include <stdlib.h>
    00006 #include <stdio.h>
    00007 // from Gaudi
    00008 #include "GaudiKernel/ToolFactory.h"
    00009 #include "GaudiKernel/MsgStream.h"
    00010 //from GSL
    00011 #include "gsl/gsl_vector.h"
    00012 #include "gsl/gsl_multimin.h"
    00013 #include "gsl/gsl_math.h"
    00014 #include "gsl/gsl_matrix.h"
    00015 #include "gsl/gsl_linalg.h"
    00016 #include "gsl/gsl_blas.h"
    00017 #include "gsl/gsl_errno.h"
    00018
    00019 // local
    00020 #include "FuncMinimum.h"
    00021
    00022 // Handle CLHEP 2.0.x move to CLHEP namespace
    00023 namespace CLHEP { }
    00024 using namespace CLHEP;
    00025
    00033 // Declaration of the Tool Factory
    00034 DECLARE_TOOL_FACTORY(FuncMinimum)
    00035
    00036 // ============================================================================
    00037 // ============================================================================
    00038 FuncMinimum::FuncMinimumMisc::FuncMinimumMisc
    00039 ( const FuncMinimum::GenFunc& func  ,
    00040   FuncMinimum::Arg&           arg   )
    00041   : m_argum ( arg    )
    00042   , m_eq    ( &func  )
    00043   , m_grad  ()
    00044 {
    00045   const size_t N = func.dimensionality () ;
    00046
    00047   for( size_t i = 0 ; i < N ; ++i )
    00048     {
    00049       Genfun::GENFUNCTION fun = func.partial(i);
    00050       m_grad.push_back (fun.clone());
    00051     }
    00052 }
    00053 // ============================================================================
    00054
    00055 // ============================================================================
    00056 FuncMinimum::FuncMinimumMisc::~FuncMinimumMisc()
    00057 {
    00058   m_grad.clear();
    00059 }
    00060 // ============================================================================
    00061
    00062 //=============================================================================
    00063 // Standard constructor, initializes variables
    00064 //=============================================================================
    00065 FuncMinimum::FuncMinimum( const std::string& type,
    00066                           const std::string& name,
    00067                           const IInterface* parent )
    00068   : base_class      ( type, name , parent )
    00069   , m_algType       ( "conjugate_fr"      )
    00070   , m_max_iter      ( 200                 )
    00071   , m_norm_gradient ( 1.0e-10             )
    00072   , m_step_size     ( 0.01                )
    00073   , m_tol           ( 1e-10               )
    00074   , m_type          ( 0                   )
    00075 {
    00077   declareProperty ( "Algorithm", m_algType       );
    00079   declareProperty ( "Iteration", m_max_iter      );
    00081   declareProperty ( "Gradient" , m_norm_gradient );
    00083   declareProperty ( "Step_size", m_step_size     );
    00085   declareProperty ( "Tol"      , m_tol           );
    00086 }
    00087 //=============================================================================
    00088 namespace
    00089 {
    00090   using namespace Genfun;
    00092   double fun_gsl ( const gsl_vector* v ,
    00093                    void * params)
    00094   {
    00095     FuncMinimum::FuncMinimumMisc * local =
    00096       static_cast <FuncMinimum::FuncMinimumMisc*> (params);
    00097     const FuncMinimum::GenFunc& eq  = *(local->equation());
    00098     Argument&                   arg = local->argument()   ;
    00099
    00100     {for (unsigned int i = 0; i < v->size; ++i)
    00101       {
    00102         arg[i] = gsl_vector_get (v, i);
    00103       }
    00104     }
    00105
    00106     return eq(arg);
    00107   }
    00108
    00110   void dfun_gsl ( const gsl_vector* v , void * params,
    00111                   gsl_vector* df )
    00112   {
    00113     FuncMinimum::FuncMinimumMisc * local =
    00114       static_cast <FuncMinimum::FuncMinimumMisc*> (params);
    00115     const FuncMinimum::Gradient& grad  = local->gradient();
    00116     Argument&                    arg   = local->argument()  ;
    00117
    00118     {for ( unsigned int i = 0; i < v->size; ++i)
    00119       {
    00120         arg[i] = gsl_vector_get (v, i);
    00121       }
    00122     }
    00123
    00124     {for( unsigned int i = 0 ; i < df->size ; ++i )
    00125
    00126       {
    00127         Genfun::GENFUNCTION f = *(grad[i]);
    00128         gsl_vector_set ( df, i, f(arg) );
    00129       }
    00130     }
    00131
    00132   }
    00133
    00136   void fdfun_gsl ( const gsl_vector* v ,
    00137                    void*       params  ,
    00138                    double*     f       ,
    00139                    gsl_vector* df )
    00140   {
    00141     *f = fun_gsl( v , params );
    00142     dfun_gsl ( v , params, df);
    00143   }
    00144
    00145 }
    00146
    00147 //=============================================================================
    00152 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc&   func  ,
    00153                                  IFuncMinimum::Arg&             arg   ) const
    00154 //=============================================================================
    00155 {
    00156   using namespace Genfun;
    00157
    00158   gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
    00159                                                  arg.dimension() );
    00160   MsgStream log( msgSvc(), name() );
    00161
    00162   FuncMinimumMisc local (func, arg);
    00163
    00164   gsl_multimin_function_fdf function;
    00165
    00166   function.f = &fun_gsl;
    00167   function.df = &dfun_gsl;
    00168   function.fdf = &fdfun_gsl;
    00169   function.n = vect.vector.size;
    00170   function.params = (void*) &local;
    00171
    00172   size_t iter   = 0 ;
    00173   int    status = 0 ;
    00174   const gsl_multimin_fdfminimizer_type *T = m_type ;
    00175
    00176   gsl_multimin_fdfminimizer *s;
    00177
    00178   s = gsl_multimin_fdfminimizer_alloc ( T, vect.vector.size);
    00179
    00180   gsl_multimin_fdfminimizer_set ( s, &function,
    00181                                   &vect.vector, m_step_size, m_tol);
    00182
    00183   for( iter = 0 ; iter < m_max_iter ; ++iter )
    00184     {
    00185       status = gsl_multimin_fdfminimizer_iterate (s);
    00186
    00187       if ( status )
    00188         {
    00189           return Error
    00190             ("Error from gsl_multimin_fdfminimizer_iterate '"
    00191              + std::string(gsl_strerror(status)) + "'") ;
    00192         }
    00193
    00194       status = gsl_multimin_test_gradient (s->gradient,
    00195                                            m_norm_gradient);
    00196
    00197
    00198       if ( status != GSL_CONTINUE ) { break; }
    00199     }
    00200
    00201   for (unsigned int i = 0; i < vect.vector.size; ++i)
    00202     {
    00203       gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
    00204     }
    00205
    00206   if (status == GSL_SUCCESS)
    00207     {
    00208       log << MSG::DEBUG
    00209           << "We stopped in the method on the " << iter
    00210           << " iteration (we have maximum "     << m_max_iter
    00211           << " iterations)"                     << endmsg;
    00212
    00213       log << "The Euclidean norm of gradient = "
    00214           << gsl_blas_dnrm2 (s->gradient)
    00215           << " by the absolute tolerance = "
    00216           << m_norm_gradient << endmsg;
    00217     }
    00218   else if (status == GSL_CONTINUE && iter <= m_max_iter )
    00219     {
    00220       return Error ( "Method finished with '"
    00221                      + std::string(gsl_strerror(status))
    00222                      +  "' error" );
    00223     }
    00224   else
    00225     {
    00226       return Error ( "Method finished with '" +
    00227                      std::string(gsl_strerror(status))
    00228                      +  "' error" );
    00229     }
    00230
    00231   gsl_multimin_fdfminimizer_free (s);
    00232
    00233   if (status)
    00234     {
    00235       return Error ( "Method finished with '"
    00236                      + std::string(gsl_strerror(status))
    00237                      +  "' error" );
    00238     }
    00239
    00240   return StatusCode::SUCCESS;
    00241 }
    00242
    00243 //=============================================================================
    00248 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc&   func  ,
    00249                                  IFuncMinimum::Arg&             arg   ,
    00250                                  IFuncMinimum::Covariance&      covar ) const
    00251 //=============================================================================
    00252 {
    00253   MsgStream log( msgSvc(), name() );
    00254
    00256   StatusCode sc = minimum (func, arg) ;
    00257   char buffer[100];
    00258   sprintf(buffer, "%lu", sc.getCode());
    00259
    00260   if (sc.isFailure())
    00261     {
    00262       return Error ("MINIMUM IS NOT FOUND. StatusCode =  '"
    00263                     + std::string(buffer) + "'", sc);
    00264     }
    00265   else
    00266     {
    00267       HepSymMatrix cov(arg.dimension(), 0);
    00268       for (unsigned int i = 0; i < arg.dimension(); ++i)
    00269         {
    00270           Genfun::GENFUNCTION f = func.partial(i);
    00271           for (unsigned int j = i; j < arg.dimension(); ++j)
    00272             {
    00273               Genfun::GENFUNCTION fij = f.partial(j);
    00274               cov(i+1, j+1) = 0.5 * fij(arg);
    00275             }
    00276         }
    00277
    00278       int inv;
    00279       covar = cov.inverse(inv);
    00280       if ( inv != 0)
    00281         {
    00282           return Error
    00283             ("Matrix of Error is not complete successful");
    00284         }
    00285
    00286       return StatusCode::SUCCESS;
    00287     }
    00288
    00289 }
    00290
    00291 //=============================================================================
    00292 StatusCode  FuncMinimum::initialize()
    00293
    00294 {
    00295   StatusCode sc = GaudiTool::initialize() ;
    00296
    00297   MsgStream log( msgSvc() , name() ) ;
    00298
    00299   if ( sc.isFailure() )
    00300     {
    00301       return Error ("Could not initiliaze base class GaudiTool", sc);
    00302     }
    00303
    00305   if(       "conjugate_fr" == m_algType )
    00306     {
    00307       m_type = gsl_multimin_fdfminimizer_conjugate_fr     ;
    00308       log << MSG::DEBUG
    00309           << "Minimization algorithm to be used: "
    00310           << "'gsl_multimin_fdfminimizer_conjugate_fr'"
    00311           << endmsg;
    00312     }
    00313   else if ( "conjugate_pr" == m_algType )
    00314     {
    00315       m_type = gsl_multimin_fdfminimizer_conjugate_pr     ;
    00316       log << MSG::DEBUG
    00317           << "Minimization algorithm to be used: "
    00318           << "'gsl_multimin_fdfminimizer_conjugate_pr'"
    00319           << endmsg;
    00320     }
    00321   else if ( "vector_bfgs" == m_algType )
    00322     {
    00323       m_type = gsl_multimin_fdfminimizer_vector_bfgs      ;
    00324       log << MSG::DEBUG
    00325           << "Minimization algorithm to be used: " <<
    00326           "'gsl_multimin_fdfminimizer_vector_bfgs'" << endmsg;
    00327     }
    00328   else if ( "steepest_descent" == m_algType )
    00329     {
    00330       m_type = gsl_multimin_fdfminimizer_steepest_descent ;
    00331       log << MSG::DEBUG
    00332           << "Minimization algorithm to be used: "
    00333           << "'gsl_multimin_fdfminimizer_steepest_descent'"
    00334           << endmsg;
    00335     }
    00336   else
    00337     {
    00338       return Error(" Unknown algorith type '"
    00339                    + std::string(m_algType) + "'");
    00340     }
    00341
    00342   return StatusCode::SUCCESS;
    00343 }
    00344 //=============================================================================
    00345 StatusCode FuncMinimum::finalize   ()
    00346 {
    00347   StatusCode sc = GaudiTool::finalize() ;
    00348
    00349   MsgStream log( msgSvc() , name() ) ;
    00350
    00351   if ( sc.isFailure() )
    00352     {
    00353       return Error("Could not finaliaze base class GaudiTool", sc);
    00354     }
    00355   return StatusCode::SUCCESS;
    00356 }
    00357 //=============================================================================
    00358 FuncMinimum::~FuncMinimum( )
    00359
    00360 {}
    00361 //=============================================================================

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines