| # Gaudi Framework, version v21r9 |

Home |  | Generated: 3 May 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         }
    00189
    00190       status = gsl_multimin_test_gradient (s->gradient,
    00191                                            m_norm_gradient);
    00192
    00193
    00194       if ( status != GSL_CONTINUE ) { break; }
    00195     }
    00196
    00197   for (unsigned int i = 0; i < vect.vector.size; ++i)
    00198     {
    00199       gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
    00200     }
    00201
    00202   if (status == GSL_SUCCESS)
    00203     {
    00204       log << MSG::DEBUG
    00205           << "We stopped in the method on the " << iter
    00206           << " iteration (we have maximum "     << m_max_iter
    00207           << " iterations)"                     << endmsg;
    00208
    00209       log << "The Euclidean norm of gradient = "
    00210           << gsl_blas_dnrm2 (s->gradient)
    00211           << " by the absolute tolerance = "
    00212           << m_norm_gradient << endmsg;
    00213     }
    00214   else if (status == GSL_CONTINUE && iter <= m_max_iter )
    00215     {
    00216       return Error ( "Method finished with '"
    00217                      + std::string(gsl_strerror(status))
    00218                      +  "' error" );
    00219     }
    00220   else
    00221     {
    00222       return Error ( "Method finished with '" +
    00223                      std::string(gsl_strerror(status))
    00224                      +  "' error" );
    00225     }
    00226
    00227   gsl_multimin_fdfminimizer_free (s);
    00228
    00229   if (status)
    00230     {
    00231       return Error ( "Method finished with '"
    00232                      + std::string(gsl_strerror(status))
    00233                      +  "' error" );
    00234     }
    00235
    00236   return StatusCode::SUCCESS;
    00237 }
    00238
    00239 //=============================================================================
    00244 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc&   func  ,
    00245                                  IFuncMinimum::Arg&             arg   ,
    00246                                  IFuncMinimum::Covariance&      covar ) const
    00247 //=============================================================================
    00248 {
    00249   MsgStream log( msgSvc(), name() );
    00250
    00252   StatusCode sc = minimum (func, arg) ;
    00253   char buffer[100];
    00254   sprintf(buffer, "%lu", sc.getCode());
    00255
    00256   if (sc.isFailure())
    00257     {
    00258       return Error ("MINIMUM IS NOT FOUND. StatusCode =  '"
    00259                     + std::string(buffer) + "'", sc);
    00260     }
    00261   else
    00262     {
    00263       HepSymMatrix cov(arg.dimension(), 0);
    00264       for (unsigned int i = 0; i < arg.dimension(); ++i)
    00265         {
    00266           Genfun::GENFUNCTION f = func.partial(i);
    00267           for (unsigned int j = i; j < arg.dimension(); ++j)
    00268             {
    00269               Genfun::GENFUNCTION fij = f.partial(j);
    00270               cov(i+1, j+1) = 0.5 * fij(arg);
    00271             }
    00272         }
    00273
    00274       int inv;
    00275       covar = cov.inverse(inv);
    00276       if ( inv != 0)
    00277         {
    00278           return Error
    00279             ("Matrix of Error is not complete successful");
    00280         }
    00281
    00282       return StatusCode::SUCCESS;
    00283     }
    00284
    00285 }
    00286
    00287 //=============================================================================
    00288 StatusCode  FuncMinimum::initialize()
    00289
    00290 {
    00291   StatusCode sc = GaudiTool::initialize() ;
    00292
    00293   MsgStream log( msgSvc() , name() ) ;
    00294
    00295   if ( sc.isFailure() )
    00296     {
    00297       return Error ("Could not initiliaze base class GaudiTool", sc);
    00298     }
    00299
    00301   if(       "conjugate_fr" == m_algType )
    00302     {
    00303       m_type = gsl_multimin_fdfminimizer_conjugate_fr     ;
    00304       log << MSG::DEBUG
    00305           << "Minimization algorithm to be used: "
    00306           << "'gsl_multimin_fdfminimizer_conjugate_fr'"
    00307           << endmsg;
    00308     }
    00309   else if ( "conjugate_pr" == m_algType )
    00310     {
    00311       m_type = gsl_multimin_fdfminimizer_conjugate_pr     ;
    00312       log << MSG::DEBUG
    00313           << "Minimization algorithm to be used: "
    00314           << "'gsl_multimin_fdfminimizer_conjugate_pr'"
    00315           << endmsg;
    00316     }
    00317   else if ( "vector_bfgs" == m_algType )
    00318     {
    00319       m_type = gsl_multimin_fdfminimizer_vector_bfgs      ;
    00320       log << MSG::DEBUG
    00321           << "Minimization algorithm to be used: " <<
    00322           "'gsl_multimin_fdfminimizer_vector_bfgs'" << endmsg;
    00323     }
    00324   else if ( "steepest_descent" == m_algType )
    00325     {
    00326       m_type = gsl_multimin_fdfminimizer_steepest_descent ;
    00327       log << MSG::DEBUG
    00328           << "Minimization algorithm to be used: "
    00329           << "'gsl_multimin_fdfminimizer_steepest_descent'"
    00330           << endmsg;
    00331     }
    00332   else
    00333     {
    00334       return Error(" Unknown algorith type '"
    00335                    + std::string(m_algType) + "'");
    00336     }
    00337
    00338   return StatusCode::SUCCESS;
    00339 }
    00340 //=============================================================================
    00341 StatusCode FuncMinimum::finalize   ()
    00342 {
    00343   StatusCode sc = GaudiTool::finalize() ;
    00344
    00345   MsgStream log( msgSvc() , name() ) ;
    00346
    00347   if ( sc.isFailure() )
    00348     {
    00349       return Error("Could not finaliaze base class GaudiTool", sc);
    00350     }
    00351   return StatusCode::SUCCESS;
    00352 }
    00353 //=============================================================================
    00354 FuncMinimum::~FuncMinimum( )
    00355
    00356 {}
    00357 //=============================================================================