| # Gaudi Framework, version v21r7 |

Home |  | Generated: 22 Jan 2010

# 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   : base_class      ( 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 {
    00073   declareProperty ( "Algorithm", m_algType       );
    00075   declareProperty ( "Iteration", m_max_iter      );
    00077   declareProperty ( "Gradient" , m_norm_gradient );
    00079   declareProperty ( "Step_size", m_step_size     );
    00081   declareProperty ( "Tol"      , m_tol           );
    00082 }
    00083 //=============================================================================
    00084 namespace
    00085 {
    00086   using namespace Genfun;
    00088   double fun_gsl ( const gsl_vector* v ,
    00089                    void * params)
    00090   {
    00091     FuncMinimum::FuncMinimumMisc * local =
    00092       static_cast <FuncMinimum::FuncMinimumMisc*> (params);
    00093     const FuncMinimum::GenFunc& eq  = *(local->equation());
    00094     Argument&                   arg = local->argument()   ;
    00095
    00096     {for (unsigned int i = 0; i < v->size; ++i)
    00097       {
    00098         arg[i] = gsl_vector_get (v, i);
    00099       }
    00100     }
    00101
    00102     return eq(arg);
    00103   }
    00104
    00106   void dfun_gsl ( const gsl_vector* v , void * params,
    00107                   gsl_vector* df )
    00108   {
    00109     FuncMinimum::FuncMinimumMisc * local =
    00110       static_cast <FuncMinimum::FuncMinimumMisc*> (params);
    00111     const FuncMinimum::Gradient& grad  = local->gradient();
    00112     Argument&                    arg   = local->argument()  ;
    00113
    00114     {for ( unsigned int i = 0; i < v->size; ++i)
    00115       {
    00116         arg[i] = gsl_vector_get (v, i);
    00117       }
    00118     }
    00119
    00120     {for( unsigned int i = 0 ; i < df->size ; ++i )
    00121
    00122       {
    00123         Genfun::GENFUNCTION f = *(grad[i]);
    00124         gsl_vector_set ( df, i, f(arg) );
    00125       }
    00126     }
    00127
    00128   }
    00129
    00132   void fdfun_gsl ( const gsl_vector* v ,
    00133                    void*       params  ,
    00134                    double*     f       ,
    00135                    gsl_vector* df )
    00136   {
    00137     *f = fun_gsl( v , params );
    00138     dfun_gsl ( v , params, df);
    00139   }
    00140
    00141 };
    00142
    00143 //=============================================================================
    00148 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc&   func  ,
    00149                                  IFuncMinimum::Arg&             arg   ) const
    00150 //=============================================================================
    00151 {
    00152   using namespace Genfun;
    00153
    00154   gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
    00155                                                  arg.dimension() );
    00156   MsgStream log( msgSvc(), name() );
    00157
    00158   FuncMinimumMisc local (func, arg);
    00159
    00160   gsl_multimin_function_fdf function;
    00161
    00162   function.f = &fun_gsl;
    00163   function.df = &dfun_gsl;
    00164   function.fdf = &fdfun_gsl;
    00165   function.n = vect.vector.size;
    00166   function.params = (void*) &local;
    00167
    00168   size_t iter   = 0 ;
    00169   int    status = 0 ;
    00170   const gsl_multimin_fdfminimizer_type *T = m_type ;
    00171
    00172   gsl_multimin_fdfminimizer *s;
    00173
    00174   s = gsl_multimin_fdfminimizer_alloc ( T, vect.vector.size);
    00175
    00176   gsl_multimin_fdfminimizer_set ( s, &function,
    00177                                   &vect.vector, m_step_size, m_tol);
    00178
    00179   for( iter = 0 ; iter < m_max_iter ; ++iter )
    00180     {
    00181       status = gsl_multimin_fdfminimizer_iterate (s);
    00182
    00183       if ( status )
    00184         {
    00185           return Error
    00186             ("Error from gsl_multimin_fdfminimizer_iterate '"
    00187              + std::string(gsl_strerror(status)) + "'") ;
    00188           break;
    00189         }
    00190
    00191       status = gsl_multimin_test_gradient (s->gradient,
    00192                                            m_norm_gradient);
    00193
    00194
    00195       if ( status != GSL_CONTINUE ) { break; }
    00196     }
    00197
    00198   for (unsigned int i = 0; i < vect.vector.size; ++i)
    00199     {
    00200       gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
    00201     }
    00202
    00203   if (status == GSL_SUCCESS)
    00204     {
    00205       log << MSG::DEBUG
    00206           << "We stopped in the method on the " << iter
    00207           << " iteration (we have maximum "     << m_max_iter
    00208           << " iterations)"                     << endmsg;
    00209
    00210       log << "The Euclidean norm of gradient = "
    00211           << gsl_blas_dnrm2 (s->gradient)
    00212           << " by the absolute tolerance = "
    00213           << m_norm_gradient << endmsg;
    00214     }
    00215   else if (status == GSL_CONTINUE && iter <= m_max_iter )
    00216     {
    00217       return Error ( "Method finished with '"
    00218                      + std::string(gsl_strerror(status))
    00219                      +  "' error" );
    00220     }
    00221   else
    00222     {
    00223       return Error ( "Method finished with '" +
    00224                      std::string(gsl_strerror(status))
    00225                      +  "' error" );
    00226     }
    00227
    00228   gsl_multimin_fdfminimizer_free (s);
    00229
    00230   if (status)
    00231     {
    00232       return Error ( "Method finished with '"
    00233                      + std::string(gsl_strerror(status))
    00234                      +  "' error" );
    00235     }
    00236
    00237   return StatusCode::SUCCESS;
    00238 };
    00239
    00240 //=============================================================================
    00245 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc&   func  ,
    00246                                  IFuncMinimum::Arg&             arg   ,
    00247                                  IFuncMinimum::Covariance&      covar ) const
    00248 //=============================================================================
    00249 {
    00250   MsgStream log( msgSvc(), name() );
    00251
    00253   StatusCode sc = minimum (func, arg) ;
    00254   char buffer[100];
    00255   sprintf(buffer, "%lu", sc.getCode());
    00256
    00257   if (sc.isFailure())
    00258     {
    00259       return Error ("MINIMUM IS NOT FOUND. StatusCode =  '"
    00260                     + std::string(buffer) + "'", sc);
    00261     }
    00262   else
    00263     {
    00264       HepSymMatrix cov(arg.dimension(), 0);
    00265       for (unsigned int i = 0; i < arg.dimension(); ++i)
    00266         {
    00267           Genfun::GENFUNCTION f = func.partial(i);
    00268           for (unsigned int j = i; j < arg.dimension(); ++j)
    00269             {
    00270               Genfun::GENFUNCTION fij = f.partial(j);
    00271               cov(i+1, j+1) = 0.5 * fij(arg);
    00272             }
    00273         }
    00274
    00275       int inv;
    00276       covar = cov.inverse(inv);
    00277       if ( inv != 0)
    00278         {
    00279           return Error
    00280             ("Matrix of Error is not complete successful");
    00281         }
    00282
    00283       return StatusCode::SUCCESS;
    00284     }
    00285
    00286 };
    00287
    00288 //=============================================================================
    00289 StatusCode  FuncMinimum::initialize()
    00290
    00291 {
    00292   StatusCode sc = GaudiTool::initialize() ;
    00293
    00294   MsgStream log( msgSvc() , name() ) ;
    00295
    00296   if ( sc.isFailure() )
    00297     {
    00298       return Error ("Could not initiliaze base class GaudiTool", sc);
    00299     }
    00300
    00302   if(       "conjugate_fr" == m_algType )
    00303     {
    00304       m_type = gsl_multimin_fdfminimizer_conjugate_fr     ;
    00305       log << MSG::DEBUG
    00306           << "Minimization algorithm to be used: "
    00307           << "'gsl_multimin_fdfminimizer_conjugate_fr'"
    00308           << endmsg;
    00309     }
    00310   else if ( "conjugate_pr" == m_algType )
    00311     {
    00312       m_type = gsl_multimin_fdfminimizer_conjugate_pr     ;
    00313       log << MSG::DEBUG
    00314           << "Minimization algorithm to be used: "
    00315           << "'gsl_multimin_fdfminimizer_conjugate_pr'"
    00316           << endmsg;
    00317     }
    00318   else if ( "vector_bfgs" == m_algType )
    00319     {
    00320       m_type = gsl_multimin_fdfminimizer_vector_bfgs      ;
    00321       log << MSG::DEBUG
    00322           << "Minimization algorithm to be used: " <<
    00323           "'gsl_multimin_fdfminimizer_vector_bfgs'" << endmsg;
    00324     }
    00325   else if ( "steepest_descent" == m_algType )
    00326     {
    00327       m_type = gsl_multimin_fdfminimizer_steepest_descent ;
    00328       log << MSG::DEBUG
    00329           << "Minimization algorithm to be used: "
    00330           << "'gsl_multimin_fdfminimizer_steepest_descent'"
    00331           << endmsg;
    00332     }
    00333   else
    00334     {
    00335       return Error(" Unknown algorith type '"
    00336                    + std::string(m_algType) + "'");
    00337     }
    00338
    00339   return StatusCode::SUCCESS;
    00340 }
    00341 //=============================================================================
    00342 StatusCode FuncMinimum::finalize   ()
    00343 {
    00344   StatusCode sc = GaudiTool::finalize() ;
    00345
    00346   MsgStream log( msgSvc() , name() ) ;
    00347
    00348   if ( sc.isFailure() )
    00349     {
    00350       return Error("Could not finaliaze base class GaudiTool", sc);
    00351     }
    00352   return StatusCode::SUCCESS;
    00353 };
    00354 //=============================================================================
    00355 FuncMinimum::~FuncMinimum( )
    00356
    00357 {}
    00358 //=============================================================================