| The Gaudi Framework v30r0 (c919700c)

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

FuncMinimum.cpp
Go to the documentation of this file.
1 // Include files
2
3 #include <cstdio>
4 #include <cstdlib>
5 // from Gaudi
6 #include "GaudiKernel/MsgStream.h"
7 // from GSL
8 #include "gsl/gsl_blas.h"
9 #include "gsl/gsl_errno.h"
10 #include "gsl/gsl_linalg.h"
11 #include "gsl/gsl_math.h"
12 #include "gsl/gsl_matrix.h"
13 #include "gsl/gsl_multimin.h"
14 #include "gsl/gsl_vector.h"
15
16 // local
17 #include "FuncMinimum.h"
18
19 // Handle CLHEP 2.0.x move to CLHEP namespace
20 namespace CLHEP
21 {
22 }
23 using namespace CLHEP;
24
32 // ============================================================================
33 // ============================================================================
34 FuncMinimum::FuncMinimumMisc::FuncMinimumMisc( const FuncMinimum::GenFunc& func, FuncMinimum::Arg& arg )
35 : m_argum( arg ), m_eq( &func ), m_grad()
36 {
37 const size_t N = func.dimensionality();
38
39 for ( size_t i = 0; i < N; ++i ) {
40 Genfun::GENFUNCTION fun = func.partial( i );
41 m_grad.push_back( fun.clone() );
42 }
43 }
44 //=============================================================================
45 namespace
46 {
47 using namespace Genfun;
49 double fun_gsl( const gsl_vector* v, void* params )
50 {
51 FuncMinimum::FuncMinimumMisc* local = static_cast<FuncMinimum::FuncMinimumMisc*>( params );
52 const FuncMinimum::GenFunc& eq = *( local->equation() );
53 Argument& arg = local->argument();
54
55 for ( unsigned int i = 0; i < v->size; ++i ) {
56 arg[i] = gsl_vector_get( v, i );
57 }
58
59 return eq( arg );
60 }
61
63 void dfun_gsl( const gsl_vector* v, void* params, gsl_vector* df )
64 {
65 FuncMinimum::FuncMinimumMisc* local = static_cast<FuncMinimum::FuncMinimumMisc*>( params );
66 const FuncMinimum::Gradient& grad = local->gradient();
67 Argument& arg = local->argument();
68
69 for ( unsigned int i = 0; i < v->size; ++i ) {
70 arg[i] = gsl_vector_get( v, i );
71 }
72
73 for ( unsigned int i = 0; i < df->size; ++i ) {
74 Genfun::GENFUNCTION f = *( grad[i] );
75 gsl_vector_set( df, i, f( arg ) );
76 }
77 }
78
81 void fdfun_gsl( const gsl_vector* v, void* params, double* f, gsl_vector* df )
82 {
83 *f = fun_gsl( v, params );
84 dfun_gsl( v, params, df );
85 }
86 }
87
88 //=============================================================================
93 StatusCode FuncMinimum::minimum( const IFuncMinimum::GenFunc& func, IFuncMinimum::Arg& arg ) const
94 //=============================================================================
95 {
96 using namespace Genfun;
97
98 gsl_vector_view vect = gsl_vector_view_array( &arg[0], arg.dimension() );
99 FuncMinimumMisc local( func, arg );
100
101 gsl_multimin_function_fdf function;
102
103 function.f = &fun_gsl;
104 function.df = &dfun_gsl;
105 function.fdf = &fdfun_gsl;
106 function.n = vect.vector.size;
107 function.params = (void*)&local;
108
109 size_t iter = 0;
110 int status = 0;
111 const gsl_multimin_fdfminimizer_type* T = m_type;
112
113 gsl_multimin_fdfminimizer* s;
114
115 s = gsl_multimin_fdfminimizer_alloc( T, vect.vector.size );
116
117 gsl_multimin_fdfminimizer_set( s, &function, &vect.vector, m_step_size, m_tol );
118
119 for ( iter = 0; iter < m_max_iter; ++iter ) {
120 status = gsl_multimin_fdfminimizer_iterate( s );
121
122 if ( status ) {
123 return Error( "Error from gsl_multimin_fdfminimizer_iterate '" + std::string( gsl_strerror( status ) ) + "'" );
124 }
125
126 status = gsl_multimin_test_gradient( s->gradient, m_norm_gradient );
127
128 if ( status != GSL_CONTINUE ) {
129 break;
130 }
131 }
132
133 for ( unsigned int i = 0; i < vect.vector.size; ++i ) {
134 gsl_vector_set( &vect.vector, i, gsl_vector_get( s->x, i ) );
135 }
136
137 if ( status == GSL_SUCCESS ) {
138 debug() << "We stopped in the method on the " << iter << " iteration (we have maximum " << m_max_iter
139 << " iterations)" << endmsg;
140
141 msgStream() << "The Euclidean norm of gradient = " << gsl_blas_dnrm2( s->gradient )
142 << " by the absolute tolerance = " << m_norm_gradient << endmsg;
143 } else if ( status == GSL_CONTINUE && iter <= m_max_iter ) {
144 return Error( "Method finished with '" + std::string( gsl_strerror( status ) ) + "' error" );
145 } else {
146 return Error( "Method finished with '" + std::string( gsl_strerror( status ) ) + "' error" );
147 }
148
149 gsl_multimin_fdfminimizer_free( s );
150
151 if ( status ) {
152 return Error( "Method finished with '" + std::string( gsl_strerror( status ) ) + "' error" );
153 }
154
155 return StatusCode::SUCCESS;
156 }
157
158 //=============================================================================
163 StatusCode FuncMinimum::minimum( const IFuncMinimum::GenFunc& func, IFuncMinimum::Arg& arg,
164 IFuncMinimum::Covariance& covar ) const
165 //=============================================================================
166 {
168 StatusCode sc = minimum( func, arg );
169
170 if ( sc.isFailure() ) {
171 return Error( "MINIMUM IS NOT FOUND. StatusCode = '" + std::to_string( sc.getCode() ) + '\'', sc );
172 }
173 HepSymMatrix cov( arg.dimension(), 0 );
174 for ( unsigned int i = 0; i < arg.dimension(); ++i ) {
175 auto f = func.partial( i );
176 for ( unsigned int j = i; j < arg.dimension(); ++j ) {
177 auto fij = f.partial( j );
178 cov( i + 1, j + 1 ) = 0.5 * fij( arg );
179 }
180 }
181
182 int inv;
183 covar = cov.inverse( inv );
184 return inv == 0 ? StatusCode::SUCCESS : Error( "Matrix of Error is not complete successful" );
185 }
186
187 //=============================================================================
188 StatusCode FuncMinimum::initialize()
189
190 {
191 StatusCode sc = GaudiTool::initialize();
192
193 if ( sc.isFailure() ) {
194 return Error( "Could not initialize base class GaudiTool", sc );
195 }
196
198 if ( "conjugate_fr" == m_algType ) {
199 m_type = gsl_multimin_fdfminimizer_conjugate_fr;
200 debug() << "Minimization algorithm to be used: "
201 << "'gsl_multimin_fdfminimizer_conjugate_fr'" << endmsg;
202 } else if ( "conjugate_pr" == m_algType ) {
203 m_type = gsl_multimin_fdfminimizer_conjugate_pr;
204 debug() << "Minimization algorithm to be used: "
205 << "'gsl_multimin_fdfminimizer_conjugate_pr'" << endmsg;
206 } else if ( "vector_bfgs" == m_algType ) {
207 m_type = gsl_multimin_fdfminimizer_vector_bfgs;
208 debug() << "Minimization algorithm to be used: "
209 << "'gsl_multimin_fdfminimizer_vector_bfgs'" << endmsg;
210 } else if ( "steepest_descent" == m_algType ) {
211 m_type = gsl_multimin_fdfminimizer_steepest_descent;
212 debug() << "Minimization algorithm to be used: "
213 << "'gsl_multimin_fdfminimizer_steepest_descent'" << endmsg;
214 } else {
215 return Error( " Unknown algorithm type '" + m_algType + "'" );
216 }
217
218 return StatusCode::SUCCESS;
219 }
220 //=============================================================================
221
222 // Declaration of the Tool Factory
223 DECLARE_COMPONENT( FuncMinimum )
StatusCode::SUCCESS
Definition: StatusCode.h:29
IFuncMinimum::Arg
Genfun::Argument Arg
Argument of function "GenFunc" (.
Definition: IFuncMinimum.h:35
FuncMinimum::FuncMinimumMisc
Definition: FuncMinimum.h:45
StatusCode::getCode
unsigned long getCode() const
Get the status code by value.
Definition: StatusCode.h:69
MsgStream.h
std::to_string
T to_string(T...args)
FuncMinimum::FuncMinimumMisc::argument
const Arg & argument() const
Definition: FuncMinimum.h:54
FuncMinimum::FuncMinimumMisc::equation
const GenFunc * equation() const
Definition: FuncMinimum.h:56
StatusCode::isFailure
bool isFailure() const
Test for a status code of FAILURE.
Definition: StatusCode.h:61
DECLARE_COMPONENT
#define DECLARE_COMPONENT(type)
Definition: PluginService.h:33
std::string
STL class.
EqSolver::m_algType
Gaudi::Property< std::string > m_algType
Definition: EqSolver.h:81
IOTest.N
int N
Definition: IOTest.py:101
FuncMinimum.h
EqSolver::m_max_iter
Gaudi::Property< double > m_max_iter
Definition: EqSolver.h:83
GaudiPluginService.cpluginsvc.func
func
Definition: cpluginsvc.py:229
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
IFuncMinimum::GenFunc
Genfun::AbsFunction GenFunc
Function which we minimize (.
Definition: IFuncMinimum.h:33
IFuncMinimum::Covariance
CLHEP::HepSymMatrix Covariance
Covariance matrix (matrix of error) (.
Definition: IFuncMinimum.h:37
FuncMinimum
The simplest concrete implementation of IFuncMinimum interface.
Definition: FuncMinimum.h:24
std::vector< const GenFunc * >
FuncMinimum::FuncMinimumMisc::gradient
const Gradient & gradient() const
Definition: FuncMinimum.h:57
CommonMessagingBase::debug
MsgStream & debug() const
shortcut for the method msgStream(MSG::DEBUG)
Definition: CommonMessaging.h:115
GaudiTool::Error
StatusCode Error(const std::string &msg, const StatusCode st=StatusCode::FAILURE, const size_t mx=10) const
Print the error message and return with the given StatusCode.
Definition: GaudiTool.h:671
CLHEP
Definition: IFuncMinimum.h:16
FuncMinimum::minimum
StatusCode minimum(const GenFunc &func, Arg &arg) const override
Find minimum of the function "GenFunc".
fun
double fun(const std::vector< double > &x)
Definition: PFuncTest.cpp:26
FuncMinimum::FuncMinimumMisc::FuncMinimumMisc
FuncMinimumMisc()=delete
CommonMessagingBase::msgStream
MsgStream & msgStream() const
Return an uninitialized MsgStream.
Definition: CommonMessaging.h:78
gaudirun.s
string s
Definition: gaudirun.py:253
EqSolver::m_type
const gsl_multiroot_fdfsolver_type * m_type
Definition: EqSolver.h:86
Genfun
CLHEP.
Definition: IEqSolver.h:13
endmsg
MsgStream & endmsg(MsgStream &s)
MsgStream Modifier: endmsg. Calls the output method of the MsgStream.
Definition: MsgStream.h:209
GaudiTool::initialize
StatusCode initialize() override
standard initialization method
Definition: GaudiTool.cpp:151
FuncMinimum::initialize
StatusCode initialize() override
Overriding initialize.
Definition: FuncMinimum.cpp:188