| The Gaudi Framework v27r1

- 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 <stdlib.h>
4 #include <stdio.h>
5 // from Gaudi
6 #include "GaudiKernel/MsgStream.h"
7 //from GSL
8 #include "gsl/gsl_vector.h"
9 #include "gsl/gsl_multimin.h"
10 #include "gsl/gsl_math.h"
11 #include "gsl/gsl_matrix.h"
12 #include "gsl/gsl_linalg.h"
13 #include "gsl/gsl_blas.h"
14 #include "gsl/gsl_errno.h"
15
16 // local
17 #include "FuncMinimum.h"
18
19 // Handle CLHEP 2.0.x move to CLHEP namespace
20 namespace CLHEP { }
21 using namespace CLHEP;
22
30 // ============================================================================
31 // ============================================================================
32 FuncMinimum::FuncMinimumMisc::FuncMinimumMisc
33 ( const FuncMinimum::GenFunc& func ,
34 FuncMinimum::Arg& arg )
35 : m_argum ( arg )
36 , m_eq ( &func )
37 , m_grad ()
38 {
39 const size_t N = func.dimensionality () ;
40
41 for( size_t i = 0 ; i < N ; ++i )
42 {
43 Genfun::GENFUNCTION fun = func.partial(i);
44 m_grad.push_back (fun.clone());
45 }
46 }
47 // ============================================================================
48
49 //=============================================================================
50 // Standard constructor, initializes variables
51 //=============================================================================
52 FuncMinimum::FuncMinimum( const std::string& type,
53 const std::string& name,
54 const IInterface* parent )
55 : base_class ( type, name , parent )
56 {
58 declareProperty ( "Algorithm", m_algType );
60 declareProperty ( "Iteration", m_max_iter );
62 declareProperty ( "Gradient" , m_norm_gradient );
64 declareProperty ( "Step_size", m_step_size );
66 declareProperty ( "Tol" , m_tol );
67 }
68 //=============================================================================
69 namespace
70 {
71 using namespace Genfun;
73 double fun_gsl ( const gsl_vector* v ,
74 void * params)
75 {
76 FuncMinimum::FuncMinimumMisc * local =
77 static_cast <FuncMinimum::FuncMinimumMisc*> (params);
78 const FuncMinimum::GenFunc& eq = *(local->equation());
79 Argument& arg = local->argument() ;
80
81 for (unsigned int i = 0; i < v->size; ++i) {
82 arg[i] = gsl_vector_get (v, i);
83 }
84
85 return eq(arg);
86 }
87
89 void dfun_gsl ( const gsl_vector* v , void * params,
90 gsl_vector* df )
91 {
92 FuncMinimum::FuncMinimumMisc * local =
93 static_cast <FuncMinimum::FuncMinimumMisc*> (params);
94 const FuncMinimum::Gradient& grad = local->gradient();
95 Argument& arg = local->argument() ;
96
97 for ( unsigned int i = 0; i < v->size; ++i) {
98 arg[i] = gsl_vector_get (v, i);
99 }
100
101
102 for( unsigned int i = 0 ; i < df->size ; ++i )
103 {
104 Genfun::GENFUNCTION f = *(grad[i]);
105 gsl_vector_set ( df, i, f(arg) );
106 }
107
108
109 }
110
113 void fdfun_gsl ( const gsl_vector* v ,
114 void* params ,
115 double* f ,
116 gsl_vector* df )
117 {
118 *f = fun_gsl( v , params );
119 dfun_gsl ( v , params, df);
120 }
121
122 }
123
124 //=============================================================================
129 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc& func ,
130 IFuncMinimum::Arg& arg ) const
131 //=============================================================================
132 {
133 using namespace Genfun;
134
135 gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
136 arg.dimension() );
137 FuncMinimumMisc local (func, arg);
138
139 gsl_multimin_function_fdf function;
140
141 function.f = &fun_gsl;
142 function.df = &dfun_gsl;
143 function.fdf = &fdfun_gsl;
144 function.n = vect.vector.size;
145 function.params = (void*) &local;
146
147 size_t iter = 0 ;
148 int status = 0 ;
149 const gsl_multimin_fdfminimizer_type *T = m_type ;
150
151 gsl_multimin_fdfminimizer *s;
152
153 s = gsl_multimin_fdfminimizer_alloc ( T, vect.vector.size);
154
155 gsl_multimin_fdfminimizer_set ( s, &function,
156 &vect.vector, m_step_size, m_tol);
157
158 for( iter = 0 ; iter < m_max_iter ; ++iter )
159 {
160 status = gsl_multimin_fdfminimizer_iterate (s);
161
162 if ( status )
163 {
164 return Error
165 ("Error from gsl_multimin_fdfminimizer_iterate '"
166 + std::string(gsl_strerror(status)) + "'") ;
167 }
168
169 status = gsl_multimin_test_gradient (s->gradient,
170 m_norm_gradient);
171
172
173 if ( status != GSL_CONTINUE ) { break; }
174 }
175
176 for (unsigned int i = 0; i < vect.vector.size; ++i)
177 {
178 gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
179 }
180
181 if (status == GSL_SUCCESS)
182 {
183 debug()
184 << "We stopped in the method on the " << iter
185 << " iteration (we have maximum " << m_max_iter
186 << " iterations)" << endmsg;
187
188 msgStream() << "The Euclidean norm of gradient = "
189 << gsl_blas_dnrm2 (s->gradient)
190 << " by the absolute tolerance = "
191 << m_norm_gradient << endmsg;
192 }
193 else if (status == GSL_CONTINUE && iter <= m_max_iter )
194 {
195 return Error ( "Method finished with '"
196 + std::string(gsl_strerror(status))
197 + "' error" );
198 }
199 else
200 {
201 return Error ( "Method finished with '" +
202 std::string(gsl_strerror(status))
203 + "' error" );
204 }
205
206 gsl_multimin_fdfminimizer_free (s);
207
208 if (status)
209 {
210 return Error ( "Method finished with '"
211 + std::string(gsl_strerror(status))
212 + "' error" );
213 }
214
215 return StatusCode::SUCCESS;
216 }
217
218 //=============================================================================
223 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc& func ,
224 IFuncMinimum::Arg& arg ,
225 IFuncMinimum::Covariance& covar ) const
226 //=============================================================================
227 {
229 StatusCode sc = minimum (func, arg) ;
230
231 if (sc.isFailure())
232 {
233 return Error ( "MINIMUM IS NOT FOUND. StatusCode = '"
234 + std::to_string( sc.getCode() ) + '\'',
235 sc );
236 }
237 else
238 {
239 HepSymMatrix cov(arg.dimension(), 0);
240 for (unsigned int i = 0; i < arg.dimension(); ++i)
241 {
242 Genfun::GENFUNCTION f = func.partial(i);
243 for (unsigned int j = i; j < arg.dimension(); ++j)
244 {
245 Genfun::GENFUNCTION fij = f.partial(j);
246 cov(i+1, j+1) = 0.5 * fij(arg);
247 }
248 }
249
250 int inv;
251 covar = cov.inverse(inv);
252 if ( inv != 0)
253 {
254 return Error
255 ("Matrix of Error is not complete successful");
256 }
257
258 return StatusCode::SUCCESS;
259 }
260
261 }
262
263 //=============================================================================
264 StatusCode FuncMinimum::initialize()
265
266 {
267 StatusCode sc = GaudiTool::initialize() ;
268
269 if ( sc.isFailure() )
270 {
271 return Error ("Could not initialize base class GaudiTool", sc);
272 }
273
275 if( "conjugate_fr" == m_algType )
276 {
277 m_type = gsl_multimin_fdfminimizer_conjugate_fr ;
278 debug()
279 << "Minimization algorithm to be used: "
280 << "'gsl_multimin_fdfminimizer_conjugate_fr'"
281 << endmsg;
282 }
283 else if ( "conjugate_pr" == m_algType )
284 {
285 m_type = gsl_multimin_fdfminimizer_conjugate_pr ;
286 debug()
287 << "Minimization algorithm to be used: "
288 << "'gsl_multimin_fdfminimizer_conjugate_pr'"
289 << endmsg;
290 }
291 else if ( "vector_bfgs" == m_algType )
292 {
293 m_type = gsl_multimin_fdfminimizer_vector_bfgs ;
294 debug()
295 << "Minimization algorithm to be used: " <<
296 "'gsl_multimin_fdfminimizer_vector_bfgs'" << endmsg;
297 }
298 else if ( "steepest_descent" == m_algType )
299 {
300 m_type = gsl_multimin_fdfminimizer_steepest_descent ;
301 debug()
302 << "Minimization algorithm to be used: "
303 << "'gsl_multimin_fdfminimizer_steepest_descent'"
304 << endmsg;
305 }
306 else
307 {
308 return Error(" Unknown algorithm type '" + m_algType + "'");
309 }
310
311 return StatusCode::SUCCESS;
312 }
313 //=============================================================================
314 StatusCode FuncMinimum::finalize ()
315 {
316 StatusCode sc = GaudiTool::finalize() ;
317
318 if ( sc.isFailure() )
319 {
320 return Error("Could not finalize base class GaudiTool", sc);
321 }
322 return StatusCode::SUCCESS;
323 }
324 //=============================================================================
325
326 // Declaration of the Tool Factory
327 DECLARE_COMPONENT(FuncMinimum)
FuncMinimum::m_norm_gradient
double m_norm_gradient
Definition: FuncMinimum.h:102
FuncMinimum::FuncMinimumMisc
Definition: FuncMinimum.h:51
gaudiComponentHelp.name
name
Definition: gaudiComponentHelp.py:28
StatusCode::getCode
unsigned long getCode() const
Get the status code by value.
Definition: StatusCode.h:93
MsgStream.h
IFuncMinimum::Arg
Genfun::Argument Arg
Argument of function "GenFunc" (.
Definition: IFuncMinimum.h:35
std::to_string
T to_string(T...args)
FuncMinimum::FuncMinimumMisc::argument
const Arg & argument() const
Definition: FuncMinimum.h:61
FuncMinimum::m_tol
double m_tol
Definition: FuncMinimum.h:104
FuncMinimum::FuncMinimumMisc::equation
const GenFunc * equation() const
Definition: FuncMinimum.h:63
merge_files.sc
int sc
Definition: merge_files.py:167
FuncMinimum::m_type
const gsl_multimin_fdfminimizer_type * m_type
Definition: FuncMinimum.h:105
StatusCode::isFailure
bool isFailure() const
Test for a status code of FAILURE.
Definition: StatusCode.h:86
DECLARE_COMPONENT
#define DECLARE_COMPONENT(type)
Definition: PluginService.h:36
FuncMinimum::minimum
StatusCode minimum(const GenFunc &func, Arg &arg) const override
Find minimum of the function "GenFunc".
std::string
STL class.
FuncMinimum::m_step_size
double m_step_size
Definition: FuncMinimum.h:103
IOTest.N
int N
Definition: IOTest.py:90
FuncMinimum.h
gaudirun.type
string type
Definition: gaudirun.py:151
IFuncMinimum::GenFunc
Genfun::AbsFunction GenFunc
Function which we minimize (.
Definition: IFuncMinimum.h:33
FuncMinimum::m_algType
std::string m_algType
Definition: FuncMinimum.h:100
GaudiPluginService.cpluginsvc.func
tuple func
Definition: cpluginsvc.py:218
StatusCode::SUCCESS
Definition: StatusCode.h:30
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
IInterface
Definition of the basic interface.
Definition: IInterface.h:234
GaudiTool::declareProperty
Property * declareProperty(const std::string &name, DataObjectHandle< T > &hndl, const std::string &doc="none") const
Definition: GaudiTool.h:834
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:64
CommonMessagingBase::debug
MsgStream & debug() const
shortcut for the method msgStream(MSG::DEBUG)
Definition: CommonMessaging.h:127
IFuncMinimum::Covariance
CLHEP::HepSymMatrix Covariance
Covariance matrix (matrix of error) (.
Definition: IFuncMinimum.h:37
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:689
CLHEP
Definition: IFuncMinimum.h:16
fun
double fun(const std::vector< double > &x)
Definition: PFuncTest.cpp:26
AlgTool
Base class from which all the concrete tool classes should be derived.
Definition: AlgTool.h:45
FuncMinimum::FuncMinimumMisc::FuncMinimumMisc
FuncMinimumMisc()=delete
GaudiTool::finalize
StatusCode finalize() override
standard finalization method
Definition: GaudiTool.cpp:179
gaudirun.s
string s
Definition: gaudirun.py:245
FuncMinimum::finalize
StatusCode finalize() override
Definition: FuncMinimum.cpp:314
Genfun
CLHEP.
Definition: IEqSolver.h:13
FuncMinimum::FuncMinimum
FuncMinimum()=delete
default constructor is private
FuncMinimum::m_max_iter
double m_max_iter
Definition: FuncMinimum.h:101
ana.i
list i
Definition: ana.py:128
endmsg
MsgStream & endmsg(MsgStream &s)
MsgStream Modifier: endmsg. Calls the output method of the MsgStream.
Definition: MsgStream.h:244
GaudiTool::initialize
StatusCode initialize() override
standard initialization method
Definition: GaudiTool.cpp:161
FuncMinimum::initialize
StatusCode initialize() override
Overriding initialize.
Definition: FuncMinimum.cpp:264