| # Gaudi Framework, version v25r0 |

Home |  | Generated: Mon Feb 17 2014

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

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
FuncMinimum.cpp
Go to the documentation of this file.
1 // $Id: FuncMinimum.cpp,v 1.4 2006/01/10 20:00:05 hmd Exp $
2 // ============================================================================
3 // Include files
4
5 #include <stdlib.h>
6 #include <stdio.h>
7 // from Gaudi
8 #include "GaudiKernel/ToolFactory.h"
9 #include "GaudiKernel/MsgStream.h"
10 //from GSL
11 #include "gsl/gsl_vector.h"
12 #include "gsl/gsl_multimin.h"
13 #include "gsl/gsl_math.h"
14 #include "gsl/gsl_matrix.h"
15 #include "gsl/gsl_linalg.h"
16 #include "gsl/gsl_blas.h"
17 #include "gsl/gsl_errno.h"
18
19 #include <sstream>
20
21 // local
22 #include "FuncMinimum.h"
23
24 // Handle CLHEP 2.0.x move to CLHEP namespace
25 namespace CLHEP { }
26 using namespace CLHEP;
27
35 // ============================================================================
36 // ============================================================================
37 FuncMinimum::FuncMinimumMisc::FuncMinimumMisc
38 ( const FuncMinimum::GenFunc& func ,
39 FuncMinimum::Arg& arg )
40 : m_argum ( arg )
41 , m_eq ( &func )
42 , m_grad ()
43 {
44 const size_t N = func.dimensionality () ;
45
46 for( size_t i = 0 ; i < N ; ++i )
47 {
48 Genfun::GENFUNCTION fun = func.partial(i);
49 m_grad.push_back (fun.clone());
50 }
51 }
52 // ============================================================================
53
54 // ============================================================================
55 FuncMinimum::FuncMinimumMisc::~FuncMinimumMisc()
56 {
57 m_grad.clear();
58 }
59 // ============================================================================
60
61 //=============================================================================
62 // Standard constructor, initializes variables
63 //=============================================================================
64 FuncMinimum::FuncMinimum( const std::string& type,
65 const std::string& name,
66 const IInterface* parent )
67 : base_class ( type, name , parent )
68 , m_algType ( "conjugate_fr" )
69 , m_max_iter ( 200 )
70 , m_norm_gradient ( 1.0e-10 )
71 , m_step_size ( 0.01 )
72 , m_tol ( 1e-10 )
73 , m_type ( 0 )
74 {
76 declareProperty ( "Algorithm", m_algType );
78 declareProperty ( "Iteration", m_max_iter );
80 declareProperty ( "Gradient" , m_norm_gradient );
82 declareProperty ( "Step_size", m_step_size );
84 declareProperty ( "Tol" , m_tol );
85 }
86 //=============================================================================
87 namespace
88 {
89 using namespace Genfun;
91 double fun_gsl ( const gsl_vector* v ,
92 void * params)
93 {
94 FuncMinimum::FuncMinimumMisc * local =
95 static_cast <FuncMinimum::FuncMinimumMisc*> (params);
96 const FuncMinimum::GenFunc& eq = *(local->equation());
97 Argument& arg = local->argument() ;
98
99 {for (unsigned int i = 0; i < v->size; ++i)
100 {
101 arg[i] = gsl_vector_get (v, i);
102 }
103 }
104
105 return eq(arg);
106 }
107
109 void dfun_gsl ( const gsl_vector* v , void * params,
110 gsl_vector* df )
111 {
112 FuncMinimum::FuncMinimumMisc * local =
113 static_cast <FuncMinimum::FuncMinimumMisc*> (params);
114 const FuncMinimum::Gradient& grad = local->gradient();
115 Argument& arg = local->argument() ;
116
117 {for ( unsigned int i = 0; i < v->size; ++i)
118 {
119 arg[i] = gsl_vector_get (v, i);
120 }
121 }
122
123 {for( unsigned int i = 0 ; i < df->size ; ++i )
124
125 {
126 Genfun::GENFUNCTION f = *(grad[i]);
127 gsl_vector_set ( df, i, f(arg) );
128 }
129 }
130
131 }
132
135 void fdfun_gsl ( const gsl_vector* v ,
136 void* params ,
137 double* f ,
138 gsl_vector* df )
139 {
140 *f = fun_gsl( v , params );
141 dfun_gsl ( v , params, df);
142 }
143
144 }
145
146 //=============================================================================
151 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc& func ,
152 IFuncMinimum::Arg& arg ) const
153 //=============================================================================
154 {
155 using namespace Genfun;
156
157 gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
158 arg.dimension() );
159 MsgStream log( msgSvc(), name() );
160
161 FuncMinimumMisc local (func, arg);
162
163 gsl_multimin_function_fdf function;
164
165 function.f = &fun_gsl;
166 function.df = &dfun_gsl;
167 function.fdf = &fdfun_gsl;
168 function.n = vect.vector.size;
169 function.params = (void*) &local;
170
171 size_t iter = 0 ;
172 int status = 0 ;
173 const gsl_multimin_fdfminimizer_type *T = m_type ;
174
175 gsl_multimin_fdfminimizer *s;
176
177 s = gsl_multimin_fdfminimizer_alloc ( T, vect.vector.size);
178
179 gsl_multimin_fdfminimizer_set ( s, &function,
180 &vect.vector, m_step_size, m_tol);
181
182 for( iter = 0 ; iter < m_max_iter ; ++iter )
183 {
184 status = gsl_multimin_fdfminimizer_iterate (s);
185
186 if ( status )
187 {
188 return Error
189 ("Error from gsl_multimin_fdfminimizer_iterate '"
190 + std::string(gsl_strerror(status)) + "'") ;
191 }
192
193 status = gsl_multimin_test_gradient (s->gradient,
194 m_norm_gradient);
195
196
197 if ( status != GSL_CONTINUE ) { break; }
198 }
199
200 for (unsigned int i = 0; i < vect.vector.size; ++i)
201 {
202 gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
203 }
204
205 if (status == GSL_SUCCESS)
206 {
207 log << MSG::DEBUG
208 << "We stopped in the method on the " << iter
209 << " iteration (we have maximum " << m_max_iter
210 << " iterations)" << endmsg;
211
212 log << "The Euclidean norm of gradient = "
213 << gsl_blas_dnrm2 (s->gradient)
214 << " by the absolute tolerance = "
215 << m_norm_gradient << endmsg;
216 }
217 else if (status == GSL_CONTINUE && iter <= m_max_iter )
218 {
219 return Error ( "Method finished with '"
220 + std::string(gsl_strerror(status))
221 + "' error" );
222 }
223 else
224 {
225 return Error ( "Method finished with '" +
226 std::string(gsl_strerror(status))
227 + "' error" );
228 }
229
230 gsl_multimin_fdfminimizer_free (s);
231
232 if (status)
233 {
234 return Error ( "Method finished with '"
235 + std::string(gsl_strerror(status))
236 + "' error" );
237 }
238
239 return StatusCode::SUCCESS;
240 }
241
242 //=============================================================================
247 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc& func ,
248 IFuncMinimum::Arg& arg ,
249 IFuncMinimum::Covariance& covar ) const
250 //=============================================================================
251 {
252 MsgStream log( msgSvc(), name() );
253
255 StatusCode sc = minimum (func, arg) ;
256
257 if (sc.isFailure())
258 {
259 std::ostringstream buffer;
260 buffer << "MINIMUM IS NOT FOUND. StatusCode = '" << sc.getCode() << '\'';
261 return Error (buffer.str(), sc);
262 }
263 else
264 {
265 HepSymMatrix cov(arg.dimension(), 0);
266 for (unsigned int i = 0; i < arg.dimension(); ++i)
267 {
268 Genfun::GENFUNCTION f = func.partial(i);
269 for (unsigned int j = i; j < arg.dimension(); ++j)
270 {
271 Genfun::GENFUNCTION fij = f.partial(j);
272 cov(i+1, j+1) = 0.5 * fij(arg);
273 }
274 }
275
276 int inv;
277 covar = cov.inverse(inv);
278 if ( inv != 0)
279 {
280 return Error
281 ("Matrix of Error is not complete successful");
282 }
283
284 return StatusCode::SUCCESS;
285 }
286
287 }
288
289 //=============================================================================
290 StatusCode FuncMinimum::initialize()
291
292 {
293 StatusCode sc = GaudiTool::initialize() ;
294
295 MsgStream log( msgSvc() , name() ) ;
296
297 if ( sc.isFailure() )
298 {
299 return Error ("Could not initialize base class GaudiTool", sc);
300 }
301
303 if( "conjugate_fr" == m_algType )
304 {
305 m_type = gsl_multimin_fdfminimizer_conjugate_fr ;
306 log << MSG::DEBUG
307 << "Minimization algorithm to be used: "
308 << "'gsl_multimin_fdfminimizer_conjugate_fr'"
309 << endmsg;
310 }
311 else if ( "conjugate_pr" == m_algType )
312 {
313 m_type = gsl_multimin_fdfminimizer_conjugate_pr ;
314 log << MSG::DEBUG
315 << "Minimization algorithm to be used: "
316 << "'gsl_multimin_fdfminimizer_conjugate_pr'"
317 << endmsg;
318 }
319 else if ( "vector_bfgs" == m_algType )
320 {
321 m_type = gsl_multimin_fdfminimizer_vector_bfgs ;
322 log << MSG::DEBUG
323 << "Minimization algorithm to be used: " <<
324 "'gsl_multimin_fdfminimizer_vector_bfgs'" << endmsg;
325 }
326 else if ( "steepest_descent" == m_algType )
327 {
328 m_type = gsl_multimin_fdfminimizer_steepest_descent ;
329 log << MSG::DEBUG
330 << "Minimization algorithm to be used: "
331 << "'gsl_multimin_fdfminimizer_steepest_descent'"
332 << endmsg;
333 }
334 else
335 {
336 return Error(" Unknown algorithm type '"
337 + std::string(m_algType) + "'");
338 }
339
340 return StatusCode::SUCCESS;
341 }
342 //=============================================================================
343 StatusCode FuncMinimum::finalize ()
344 {
345 StatusCode sc = GaudiTool::finalize() ;
346
347 MsgStream log( msgSvc() , name() ) ;
348
349 if ( sc.isFailure() )
350 {
351 return Error("Could not finalize base class GaudiTool", sc);
352 }
353 return StatusCode::SUCCESS;
354 }
355 //=============================================================================
356 FuncMinimum::~FuncMinimum( )
357
358 {}
359 //=============================================================================
360
361 // Declaration of the Tool Factory
362 DECLARE_COMPONENT(FuncMinimum)