| # Gaudi Framework, version v22r4 |

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

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