| The Gaudi Framework v28r2-patches (bc4ad17a)

- 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 namespace
49 {
50 using namespace Genfun;
52 double fun_gsl ( const gsl_vector* v ,
53 void * params)
54 {
55 FuncMinimum::FuncMinimumMisc * local =
56 static_cast <FuncMinimum::FuncMinimumMisc*> (params);
57 const FuncMinimum::GenFunc& eq = *(local->equation());
58 Argument& arg = local->argument() ;
59
60 for (unsigned int i = 0; i < v->size; ++i) {
61 arg[i] = gsl_vector_get (v, i);
62 }
63
64 return eq(arg);
65 }
66
68 void dfun_gsl ( const gsl_vector* v , void * params,
69 gsl_vector* df )
70 {
71 FuncMinimum::FuncMinimumMisc * local =
72 static_cast <FuncMinimum::FuncMinimumMisc*> (params);
73 const FuncMinimum::Gradient& grad = local->gradient();
74 Argument& arg = local->argument() ;
75
76 for ( unsigned int i = 0; i < v->size; ++i) {
77 arg[i] = gsl_vector_get (v, i);
78 }
79
80
81 for( unsigned int i = 0 ; i < df->size ; ++i )
82 {
83 Genfun::GENFUNCTION f = *(grad[i]);
84 gsl_vector_set ( df, i, f(arg) );
85 }
86
87
88 }
89
92 void fdfun_gsl ( const gsl_vector* v ,
93 void* params ,
94 double* f ,
95 gsl_vector* df )
96 {
97 *f = fun_gsl( v , params );
98 dfun_gsl ( v , params, df);
99 }
100
101 }
102
103 //=============================================================================
108 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc& func ,
109 IFuncMinimum::Arg& arg ) const
110 //=============================================================================
111 {
112 using namespace Genfun;
113
114 gsl_vector_view vect = gsl_vector_view_array ( &arg[0] ,
115 arg.dimension() );
116 FuncMinimumMisc local (func, arg);
117
118 gsl_multimin_function_fdf function;
119
120 function.f = &fun_gsl;
121 function.df = &dfun_gsl;
122 function.fdf = &fdfun_gsl;
123 function.n = vect.vector.size;
124 function.params = (void*) &local;
125
126 size_t iter = 0 ;
127 int status = 0 ;
128 const gsl_multimin_fdfminimizer_type *T = m_type ;
129
130 gsl_multimin_fdfminimizer *s;
131
132 s = gsl_multimin_fdfminimizer_alloc ( T, vect.vector.size);
133
134 gsl_multimin_fdfminimizer_set ( s, &function,
135 &vect.vector, m_step_size, m_tol);
136
137 for( iter = 0 ; iter < m_max_iter ; ++iter )
138 {
139 status = gsl_multimin_fdfminimizer_iterate (s);
140
141 if ( status )
142 {
143 return Error
144 ("Error from gsl_multimin_fdfminimizer_iterate '"
145 + std::string(gsl_strerror(status)) + "'") ;
146 }
147
148 status = gsl_multimin_test_gradient (s->gradient,
149 m_norm_gradient);
150
151
152 if ( status != GSL_CONTINUE ) { break; }
153 }
154
155 for (unsigned int i = 0; i < vect.vector.size; ++i)
156 {
157 gsl_vector_set (&vect.vector, i, gsl_vector_get (s->x, i));
158 }
159
160 if (status == GSL_SUCCESS)
161 {
162 debug()
163 << "We stopped in the method on the " << iter
164 << " iteration (we have maximum " << m_max_iter
165 << " iterations)" << endmsg;
166
167 msgStream() << "The Euclidean norm of gradient = "
168 << gsl_blas_dnrm2 (s->gradient)
169 << " by the absolute tolerance = "
170 << m_norm_gradient << endmsg;
171 }
172 else if (status == GSL_CONTINUE && iter <= m_max_iter )
173 {
174 return Error ( "Method finished with '"
175 + std::string(gsl_strerror(status))
176 + "' error" );
177 }
178 else
179 {
180 return Error ( "Method finished with '" +
181 std::string(gsl_strerror(status))
182 + "' error" );
183 }
184
185 gsl_multimin_fdfminimizer_free (s);
186
187 if (status)
188 {
189 return Error ( "Method finished with '"
190 + std::string(gsl_strerror(status))
191 + "' error" );
192 }
193
194 return StatusCode::SUCCESS;
195 }
196
197 //=============================================================================
202 StatusCode FuncMinimum::minimum (const IFuncMinimum::GenFunc& func ,
203 IFuncMinimum::Arg& arg ,
204 IFuncMinimum::Covariance& covar ) const
205 //=============================================================================
206 {
208 StatusCode sc = minimum (func, arg) ;
209
210 if (sc.isFailure())
211 {
212 return Error ( "MINIMUM IS NOT FOUND. StatusCode = '"
213 + std::to_string( sc.getCode() ) + '\'',
214 sc );
215 }
216 else
217 {
218 HepSymMatrix cov(arg.dimension(), 0);
219 for (unsigned int i = 0; i < arg.dimension(); ++i)
220 {
221 Genfun::GENFUNCTION f = func.partial(i);
222 for (unsigned int j = i; j < arg.dimension(); ++j)
223 {
224 Genfun::GENFUNCTION fij = f.partial(j);
225 cov(i+1, j+1) = 0.5 * fij(arg);
226 }
227 }
228
229 int inv;
230 covar = cov.inverse(inv);
231 if ( inv != 0)
232 {
233 return Error
234 ("Matrix of Error is not complete successful");
235 }
236
237 return StatusCode::SUCCESS;
238 }
239
240 }
241
242 //=============================================================================
243 StatusCode FuncMinimum::initialize()
244
245 {
246 StatusCode sc = GaudiTool::initialize() ;
247
248 if ( sc.isFailure() )
249 {
250 return Error ("Could not initialize base class GaudiTool", sc);
251 }
252
254 if( "conjugate_fr" == m_algType )
255 {
256 m_type = gsl_multimin_fdfminimizer_conjugate_fr ;
257 debug()
258 << "Minimization algorithm to be used: "
259 << "'gsl_multimin_fdfminimizer_conjugate_fr'"
260 << endmsg;
261 }
262 else if ( "conjugate_pr" == m_algType )
263 {
264 m_type = gsl_multimin_fdfminimizer_conjugate_pr ;
265 debug()
266 << "Minimization algorithm to be used: "
267 << "'gsl_multimin_fdfminimizer_conjugate_pr'"
268 << endmsg;
269 }
270 else if ( "vector_bfgs" == m_algType )
271 {
272 m_type = gsl_multimin_fdfminimizer_vector_bfgs ;
273 debug()
274 << "Minimization algorithm to be used: " <<
275 "'gsl_multimin_fdfminimizer_vector_bfgs'" << endmsg;
276 }
277 else if ( "steepest_descent" == m_algType )
278 {
279 m_type = gsl_multimin_fdfminimizer_steepest_descent ;
280 debug()
281 << "Minimization algorithm to be used: "
282 << "'gsl_multimin_fdfminimizer_steepest_descent'"
283 << endmsg;
284 }
285 else
286 {
287 return Error(" Unknown algorithm type '" + m_algType + "'");
288 }
289
290 return StatusCode::SUCCESS;
291 }
292 //=============================================================================
293 StatusCode FuncMinimum::finalize ()
294 {
295 StatusCode sc = GaudiTool::finalize() ;
296
297 if ( sc.isFailure() )
298 {
299 return Error("Could not finalize base class GaudiTool", sc);
300 }
301 return StatusCode::SUCCESS;
302 }
303 //=============================================================================
304
305 // Declaration of the Tool Factory
306 DECLARE_COMPONENT(FuncMinimum)
StatusCode::SUCCESS
Definition: StatusCode.h:30
FuncMinimum::FuncMinimumMisc
Definition: FuncMinimum.h:48
StatusCode::getCode
unsigned long getCode() const
Get the status code by value.
Definition: StatusCode.h:91
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:57
FuncMinimum::FuncMinimumMisc::equation
const GenFunc * equation() const
Definition: FuncMinimum.h:59
StatusCode::isFailure
bool isFailure() const
Test for a status code of FAILURE.
Definition: StatusCode.h:84
DECLARE_COMPONENT
#define DECLARE_COMPONENT(type)
Definition: PluginService.h:36
std::string
STL class.
EqSolver::m_algType
Gaudi::Property< std::string > m_algType
Definition: EqSolver.h:83
IOTest.N
int N
Definition: IOTest.py:90
FuncMinimum.h
EqSolver::m_max_iter
Gaudi::Property< double > m_max_iter
Definition: EqSolver.h:85
IFuncMinimum::GenFunc
Genfun::AbsFunction GenFunc
Function which we minimize (.
Definition: IFuncMinimum.h:33
GaudiPluginService.cpluginsvc.func
func
Definition: cpluginsvc.py:218
StatusCode
This class is used for returning status codes from appropriate routines.
Definition: StatusCode.h:26
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:60
CommonMessagingBase::debug
MsgStream & debug() const
shortcut for the method msgStream(MSG::DEBUG)
Definition: CommonMessaging.h:126
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: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
GaudiTool::finalize
StatusCode finalize() override
standard finalization method
Definition: GaudiTool.cpp:175
CommonMessagingBase::msgStream
MsgStream & msgStream() const
Return an uninitialized MsgStream.
Definition: CommonMessaging.h:88
gaudirun.s
string s
Definition: gaudirun.py:245
EqSolver::m_type
const gsl_multiroot_fdfsolver_type * m_type
Definition: EqSolver.h:88
FuncMinimum::finalize
StatusCode finalize() override
Definition: FuncMinimum.cpp:293
Genfun
CLHEP.
Definition: IEqSolver.h:13
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:157
FuncMinimum::initialize
StatusCode initialize() override
Overriding initialize.
Definition: FuncMinimum.cpp:243