| # Gaudi Framework, version v20r3
Generated: 24 Nov 2008

# 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
    00029 // Declaration of the Tool Factory
    00030 DECLARE_TOOL_FACTORY(FuncMinimum)
    00031
    00032 // ============================================================================
    00033 // ============================================================================
    00034 FuncMinimum::FuncMinimumMisc::FuncMinimumMisc
    00035 ( const FuncMinimum::GenFunc& func  ,
    00036   FuncMinimum::Arg&           arg   )
    00037   : m_argum ( arg    )
    00038   , m_eq    ( &func  )
    00039   , m_grad  ()
    00040 {
    00041   const size_t N = func.dimensionality () ;
    00042
    00043   for( size_t i = 0 ; i < N ; ++i )
    00044     {
    00045       Genfun::GENFUNCTION fun = func.partial(i);
    00046       m_grad.push_back (fun.clone());
    00047     }
    00048 };
    00049 // ============================================================================
    00050
    00051 // ============================================================================
    00052 FuncMinimum::FuncMinimumMisc::~FuncMinimumMisc()
    00053 {
    00054   m_grad.clear();
    00055 };
    00056 // ============================================================================
    00057
    00058 //=============================================================================
    00059 // Standard constructor, initializes variables
    00060 //=============================================================================
    00061 FuncMinimum::FuncMinimum( const std::string& type,
    00062                           const std::string& name,
    00063                           const IInterface* parent )
    00064   : GaudiTool       ( type, name , parent )
    00065   , m_algType       ( "conjugate_fr"      )
    00066   , m_max_iter      ( 200                 )
    00067   , m_norm_gradient ( 1.0e-10             )
    00068   , m_step_size     ( 0.01                )
    00069   , m_tol           ( 1e-10               )
    00070   , m_type          ( 0                   )
    00071 {
    00072
    00074   declareInterface <IFuncMinimum>(this);
    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           break;
    00193         }
    00194
    00195       status = gsl_multimin_test_gradient (s->gradient,
    00196                                            m_norm_gradient);
    00197
    00198
    00199       if ( status != GSL_CONTINUE ) { break; }
    00200     }
    00201
    00202   for (unsigned int i = 0; i < vect.vector.size; ++i)
    00203     {
    00204       gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
    00205     }
    00206
    00207   if (status == GSL_SUCCESS)
    00208     {
    00209       log << MSG::DEBUG
    00210           << "We stopped in the method on the " << iter
    00211           << " iteration (we have maximum "     << m_max_iter
    00212           << " iterations)"                     << endreq;
    00213
    00214       log << "The Euclidean norm of gradient = "
    00215           << gsl_blas_dnrm2 (s->gradient)
    00216           << " by the absolute tolerance = "
    00217           << m_norm_gradient << endreq;
    00218     }
    00219   else if (status == GSL_CONTINUE && iter <= m_max_iter )
    00220     {
    00221       return Error ( "Method finished with '"
    00222                      + std::string(gsl_strerror(status))
    00223                      +  "' error" );
    00224     }
    00225   else
    00226     {
    00227       return Error ( "Method finished with '" +
    00228                      std::string(gsl_strerror(status))
    00229                      +  "' error" );
    00230     }
    00231
    00232   gsl_multimin_fdfminimizer_free (s);
    00233
    00234   if (status)
    00235     {
    00236       return Error ( "Method finished with '"
    00237                      + std::string(gsl_strerror(status))
    00238                      +  "' error" );
    00239     }
    00240
    00241   return StatusCode::SUCCESS;
    00242 };
    00243
    00244 //=============================================================================
    00249 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc&   func  ,
    00250                                  IFuncMinimum::Arg&             arg   ,
    00251                                  IFuncMinimum::Covariance&      covar ) const
    00252 //=============================================================================
    00253 {
    00254   MsgStream log( msgSvc(), name() );
    00255
    00257   StatusCode sc = minimum (func, arg) ;
    00258   char buffer[100];
    00259   sprintf(buffer, "%lu", sc.getCode());
    00260
    00261   if (sc.isFailure())
    00262     {
    00263       return Error ("MINIMUM IS NOT FOUND. StatusCode =  '"
    00264                     + std::string(buffer) + "'", sc);
    00265     }
    00266   else
    00267     {
    00268       HepSymMatrix cov(arg.dimension(), 0);
    00269       for (unsigned int i = 0; i < arg.dimension(); ++i)
    00270         {
    00271           Genfun::GENFUNCTION f = func.partial(i);
    00272           for (unsigned int j = i; j < arg.dimension(); ++j)
    00273             {
    00274               Genfun::GENFUNCTION fij = f.partial(j);
    00275               cov(i+1, j+1) = 0.5 * fij(arg);
    00276             }
    00277         }
    00278
    00279       int inv;
    00280       covar = cov.inverse(inv);
    00281       if ( inv != 0)
    00282         {
    00283           return Error
    00284             ("Matrix of Error is not complete successful");
    00285         }
    00286
    00287       return StatusCode::SUCCESS;
    00288     }
    00289
    00290 };
    00291
    00292 //=============================================================================
    00293 StatusCode  FuncMinimum::initialize()
    00294
    00295 {
    00296   StatusCode sc = GaudiTool::initialize() ;
    00297
    00298   MsgStream log( msgSvc() , name() ) ;
    00299
    00300   if ( sc.isFailure() )
    00301     {
    00302       return Error ("Could not initiliaze base class GaudiTool", sc);
    00303     }
    00304
    00306   if(       "conjugate_fr" == m_algType )
    00307     {
    00308       m_type = gsl_multimin_fdfminimizer_conjugate_fr     ;
    00309       log << MSG::DEBUG
    00310           << "Minimization algorithm to be used: "
    00311           << "'gsl_multimin_fdfminimizer_conjugate_fr'"
    00312           << endreq;
    00313     }
    00314   else if ( "conjugate_pr" == m_algType )
    00315     {
    00316       m_type = gsl_multimin_fdfminimizer_conjugate_pr     ;
    00317       log << MSG::DEBUG
    00318           << "Minimization algorithm to be used: "
    00319           << "'gsl_multimin_fdfminimizer_conjugate_pr'"
    00320           << endreq;
    00321     }
    00322   else if ( "vector_bfgs" == m_algType )
    00323     {
    00324       m_type = gsl_multimin_fdfminimizer_vector_bfgs      ;
    00325       log << MSG::DEBUG
    00326           << "Minimization algorithm to be used: " <<
    00327           "'gsl_multimin_fdfminimizer_vector_bfgs'" << endreq;
    00328     }
    00329   else if ( "steepest_descent" == m_algType )
    00330     {
    00331       m_type = gsl_multimin_fdfminimizer_steepest_descent ;
    00332       log << MSG::DEBUG
    00333           << "Minimization algorithm to be used: "
    00334           << "'gsl_multimin_fdfminimizer_steepest_descent'"
    00335           << endreq;
    00336     }
    00337   else
    00338     {
    00339       return Error(" Unknown algorith type '"
    00340                    + std::string(m_algType) + "'");
    00341     }
    00342
    00343   return StatusCode::SUCCESS;
    00344 }
    00345 //=============================================================================
    00346 StatusCode FuncMinimum::finalize   ()
    00347 {
    00348   StatusCode sc = GaudiTool::finalize() ;
    00349
    00350   MsgStream log( msgSvc() , name() ) ;
    00351
    00352   if ( sc.isFailure() )
    00353     {
    00354       return Error("Could not finaliaze base class GaudiTool", sc);
    00355     }
    00356   return StatusCode::SUCCESS;
    00357 };
    00358 //=============================================================================
    00359 FuncMinimum::~FuncMinimum( )
    00360
    00361 {}
    00362 //=============================================================================