aimsalgo 6.0.0
Neuroimaging image processing
powell.h
Go to the documentation of this file.
1/* This software and supporting documentation are distributed by
2 * Institut Federatif de Recherche 49
3 * CEA/NeuroSpin, Batiment 145,
4 * 91191 Gif-sur-Yvette cedex
5 * France
6 *
7 * This software is governed by the CeCILL-B license under
8 * French law and abiding by the rules of distribution of free software.
9 * You can use, modify and/or redistribute the software under the
10 * terms of the CeCILL-B license as circulated by CEA, CNRS
11 * and INRIA at the following URL "http://www.cecill.info".
12 *
13 * As a counterpart to the access to the source code and rights to copy,
14 * modify and redistribute granted by the license, users are provided only
15 * with a limited warranty and the software's author, the holder of the
16 * economic rights, and the successive licensors have only limited
17 * liability.
18 *
19 * In this respect, the user's attention is drawn to the risks associated
20 * with loading, using, modifying and/or developing or reproducing the
21 * software by the user in light of its specific status of free software,
22 * that may mean that it is complicated to manipulate, and that also
23 * therefore means that it is reserved for developers and experienced
24 * professionals having in-depth computer knowledge. Users are therefore
25 * encouraged to load and test the software's suitability as regards their
26 * requirements in conditions enabling the security of their systems and/or
27 * data to be ensured and, more generally, to use and operate it in the
28 * same conditions as regards security.
29 *
30 * The fact that you are presently reading this means that you have had
31 * knowledge of the CeCILL-B license and that you accept its terms.
32 */
33
34
35#ifndef AIMS_OPTIMIZATION_POWELL_H
36#define AIMS_OPTIMIZATION_POWELL_H
37
38#include <aims/optimization/optimizer.h>
39#include <aims/optimization/objective.h>
40#include <aims/vector/vector.h>
41#include <algorithm>
42#include <list>
43#include <math.h>
44
45
46#define GOLD 1.618034
47#define GLIMIT 100.0
48#define TINY 1.0e-20
49#define SIGN( a, b ) ( (b) >= 0.0 ? fabs( a ) : -fabs( a ) )
50#define ITMAX 100
51#define TOL 2.0e-4 // tolerance passed to brent
52#define ZEPS 1.0e-10
53#define CGOLD 0.3819660
54#define SQR( a ) ( ( a ) * ( a ) )
55
56
57//
58// class PowellOptimizer
59//
60template < class T, int D >
61class PowellOptimizer : public Optimizer< T, D >
62{
63
64 public:
65 PowellOptimizer( const ObjectiveFunc< T, D >& func, T error );
66 virtual ~PowellOptimizer();
67
68 AimsVector< T, D > doit( const AimsVector< T, D > & pinit,
69 const AimsVector< T, D > & deltaP );
70
71 const AimsVector< T, D >& getDirection( int which ) const;
72 const T& getOptCost() const;
73 const std::list< T >& getCostEvolution() const;
74
75 private:
76
77 void linmin( AimsVector< T, D >& p, AimsVector< T, D >& xi );
78 void mnbrak( T& ax, T& bx, T& cx, T& fa, T& fb, T& fc );
79 T brent( T ax, T bx, T cx, T& fb, T& xmin );
80
81 AimsVector< T, D >* _xi;
82 T _optCost;
83 std::list< T > _costs;
84 T _fret;
85 AimsVector< T, D > _pcom;
86 AimsVector< T, D > _xicom;
87
88};
89
90
91template < class T, int D >
92inline
93PowellOptimizer< T, D >::PowellOptimizer( const ObjectiveFunc< T, D >& func,
94 T error )
95 : Optimizer< T, D >( func, error )
96{
97
98 _xi = new AimsVector< T, D >[ D ];
99
100}
101
102
103template < class T, int D >
104inline
105PowellOptimizer< T, D >::~PowellOptimizer()
106{
107
108 delete [] _xi;
109
110}
111
112
113template < class T, int D >
114inline
115AimsVector< T, D >
116PowellOptimizer< T, D >::doit( const AimsVector< T, D > & pinit,
117 const AimsVector< T, D > & deltaP )
118{
119
120 AimsVector< T, D > p( pinit ), pt, ptt, xit;
121 T fp, del, fptt, t;
122 int ibig;
123
124 _optCost = ( T )0;
125 _costs.clear();
126
127 for ( int i = 0; i < D; i++ )
128 for ( int j = 0; j < D; j++ )
129 _xi[ i ][ j ] = ( i == j ? deltaP[ i ] : ( T )0 );
130
131 _fret = this->_func.eval( p );
132
133 carto::AttributedObject monitoredParams("monitored_params") ;
134 this->_func.fillMonitoredParameters( monitoredParams ) ;
135 monitoredParams.setProperty("info", std::string( "first evaluation") ) ;
136 this->_probe->test( p, monitoredParams, &_fret) ;
137 pt = p;
138
139 for ( int iter = 0; ; iter++ )
140 {
141
142 fp = _fret;
143 ibig = 0;
144 del = ( T )0;
145 for ( int i = 0; i < D; i++ )
146 {
147
148 xit = _xi[ i ];
149 fptt = _fret;
150 linmin( p, xit );
151 if ( ( T )fabs( fptt - _fret ) > del )
152 {
153
154 del = ( T )fabs( fptt - _fret );
155 ibig = i;
156
157 }
158
159 }
160
161 this->_func.fillMonitoredParameters( monitoredParams ) ;
162 monitoredParams.setProperty("info", std::string( "iteration step") ) ;
163 this->_probe->iteration( p, monitoredParams, &_fret) ;
164
165
166 if ( 2.0 * fabs( fp - _fret ) <= this->_error * ( fabs( fp ) + fabs( _fret ) ) )
167 {
168
169 _optCost = _fret;
170 this->_probe->end() ;
171 return p;
172
173 }
174 if( iter == 200 ){
175 std::cerr << "Reached maximal iteration" << std::endl ;
176 _optCost = _fret;
177 this->_probe->end() ;
178 return p;
179 }
180 ptt = p * ( T )2 - pt;
181 xit = p - pt;
182 pt = p;
183
184 fptt = this->_func.eval( ptt );
185 this->_func.fillMonitoredParameters( monitoredParams ) ;
186 monitoredParams.setProperty( "info",
187 std::string( "conjugate directions investigation") );
188 this->_probe->test( ptt, monitoredParams, &_fret) ;
189 if ( fptt < fp )
190 {
191 t = ( T )( 2.0 * ( fp - 2.0 * _fret + fptt ) * SQR( fp - _fret - del )
192 - del * SQR( fp - fptt ) );
193 if ( t < ( T )0 )
194 {
195
196 _fret = this->_func.eval( p );
197 this->_func.fillMonitoredParameters( monitoredParams ) ;
198 monitoredParams.setProperty( "info",
199 std::string( "conjugate directions computing") ) ;
200 this->_probe->test( p, monitoredParams, &_fret) ;
201 linmin( p, xit );
202 _xi[ ibig ] = _xi[ D - 1 ];
203 _xi[ D - 1 ] = xit;
204
205 }
206
207 }
208 _costs.push_back( _fret );
209 }
210
211 this->_probe->end() ;
212
213 return p;
214
215}
216
217
218template < class T, int D >
219inline
220const AimsVector< T, D >&
221PowellOptimizer< T, D >::getDirection( int which ) const
222{
223
224 return _xi[ which ];
225
226}
227
228
229template < class T, int D >
230inline
231const T& PowellOptimizer< T, D >::getOptCost() const
232{
233
234 return _optCost;
235
236}
237
238
239template < class T, int D >
240inline
241const std::list< T >& PowellOptimizer< T, D >::getCostEvolution() const
242{
243
244 return _costs;
245
246}
247
248
249template < class T, int D >
250inline
251void PowellOptimizer< T, D >::linmin( AimsVector<T,D>& p, AimsVector<T,D>& xi )
252{
253
254 T xx, xmin, fx, fb, fa, bx, ax;
255
256 _pcom = p;
257 _xicom = xi;
258 ax = ( T )0;
259 xx = ( T )1;
260 fa = _fret;
261 mnbrak( ax, xx, bx, fa, fx, fb );
262 _fret = brent( ax, xx, bx, fx, xmin );
263 xi *= xmin;
264 p += xi;
265
266}
267
268
269template < class T, int D >
270inline
271void PowellOptimizer< T, D >::mnbrak( T& ax, T& bx, T& cx, T& fa, T& fb, T& fc )
272{
273
274 T ulim, u, r, q, fu, tmp;
275
276 fb = this->_func.eval( _pcom + _xicom * bx );
277
278 carto::AttributedObject monitoredParams("monitored_params") ;
279 this->_func.fillMonitoredParameters( monitoredParams ) ;
280 monitoredParams.setProperty("info", std::string( "mnbrak_bx") ) ;
281 this->_probe->test( _pcom + _xicom * bx, monitoredParams, &fb) ;
282
283 if ( fb > fa )
284 {
285
286 tmp = ax;
287 ax = bx;
288 bx = tmp;
289 tmp = fb;
290 fb = fa;
291 fa = tmp;
292
293 }
294 cx = bx + ( T )GOLD * ( bx - ax );
295 fc = this->_func.eval( _pcom + _xicom * cx );
296
297 this->_func.fillMonitoredParameters( monitoredParams ) ;
298 monitoredParams.setProperty("info", std::string( "mnbrak_cx") ) ;
299 this->_probe->test( _pcom + _xicom * cx, monitoredParams, &fc) ;
300
301
302 while ( fb > fc )
303 {
304
305 r = ( bx - ax ) * ( fb - fc );
306 q = ( bx - cx ) * ( fb - fa );
307 u = bx - ( ( bx - cx ) * q - ( bx - ax ) * r ) /
308 ( 2.0 * SIGN( std::max( ( T )fabs( q - r ), ( T )TINY ) ,
309 q - r ) );
310 ulim = bx + ( T )GLIMIT * ( cx - bx );
311 if ( ( bx - u ) * ( u - cx ) > ( T )0 )
312 {
313
314 fu = this->_func.eval( _pcom + _xicom * u );
315
316 this->_func.fillMonitoredParameters( monitoredParams ) ;
317 monitoredParams.setProperty("info", std::string( "mnbrak_u") ) ;
318 this->_probe->test( _pcom + _xicom * u, monitoredParams, &fu) ;
319
320 if ( fu < fc )
321 {
322
323 ax = bx;
324 bx = u;
325 fa = fb;
326 fb = fu;
327 return;
328
329 }
330 else if ( fu > fb )
331 {
332
333 cx = u;
334 fc = fu;
335 return;
336
337 }
338 u = cx + ( T )GOLD * ( cx - bx );
339 fu = this->_func.eval( _pcom + _xicom * u );
340
341 this->_func.fillMonitoredParameters( monitoredParams ) ;
342 monitoredParams.setProperty("info", std::string( "mnbrak_u") ) ;
343 this->_probe->test( _pcom + _xicom * u, monitoredParams, &fu) ;
344
345 }
346 else if ( ( cx - u ) * ( u - ulim ) > ( T )0 )
347 {
348
349 fu = this->_func.eval( _pcom + _xicom * u );
350
351 this->_func.fillMonitoredParameters( monitoredParams ) ;
352 monitoredParams.setProperty("info", std::string( "mnbrak_u") ) ;
353 this->_probe->test( _pcom + _xicom * u, monitoredParams, &fu) ;
354
355 if ( fu < fc )
356 {
357
358 bx = cx;
359 cx = u;
360 u = cx + ( T )GOLD * ( cx - bx );
361 fb = fc;
362 fc = fu;
363 fu = this->_func.eval( _pcom + _xicom * u );
364
365 }
366
367 }
368 else if ( ( u - ulim ) * ( ulim - cx ) >= ( T )0 )
369 {
370
371 u = ulim;
372 fu = this->_func.eval( _pcom + _xicom * u );
373
374 this->_func.fillMonitoredParameters( monitoredParams ) ;
375 monitoredParams.setProperty("info", std::string( "mnbrak_u") ) ;
376 this->_probe->test( _pcom + _xicom * u, monitoredParams, &fu) ;
377
378 }
379 else
380 {
381
382 u = cx + ( T )GOLD * ( cx - bx );
383 fu = this->_func.eval( _pcom + _xicom * u );
384
385 this->_func.fillMonitoredParameters( monitoredParams ) ;
386 monitoredParams.setProperty("info", std::string( "mnbrak_u") ) ;
387 this->_probe->test( _pcom + _xicom * u, monitoredParams, &fu) ;
388
389
390 }
391 ax = bx;
392 bx = cx;
393 cx = u;
394 fa = fb;
395 fb = fc;
396 fc = fu;
397
398 }
399
400}
401
402
403template < class T, int D >
404inline
405T PowellOptimizer< T, D >::brent( T ax, T bx, T cx, T& fb, T& xmin )
406{
407
408 T a, b, d, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, x, xm, e;
409
410 d = ( T )0;
411 e = ( T )0;
412 a = ( ax < cx ? ax : cx );
413 b = ( ax > cx ? ax : cx );
414 x = w = v = bx;
415 fw = fv = fx = fb;
416 for ( int iter = 0; iter < ITMAX; iter++ )
417 {
418
419 xm = ( T )( 0.5 * ( a + b ) );
420 tol1 = ( T )( TOL * fabs( x ) + ZEPS );
421 tol2 = ( T )2 * tol1;
422 if ( fabs( x - xm ) <= tol2 - 0.5 * ( b - a ) )
423 {
424
425 xmin = x;
426 return fx;
427
428 }
429 if ( fabs( e ) > tol1 )
430 {
431
432 r = ( x - w ) * ( fx - fv );
433 q = ( x - v ) * ( fx - fw );
434 p = ( x - v ) * q - ( x - w ) * r;
435 q = ( T )( 2.0 * ( q - r ) );
436 if ( q > ( T )0 )
437 p = -p;
438 q = ( T )fabs( q );
439 etemp = e;
440 e = d;
441 if ( fabs( p ) >= fabs( 0.5 * q * etemp ) ||
442 p <= q * ( a - x ) || p >= q * ( b - x ) )
443 {
444
445 e = ( x >= xm ? a - x : b - x );
446 d =( T )CGOLD * e;
447
448 }
449 else
450 {
451
452 d = p / q;
453 u = x + d;
454 if ( u - a < tol2 || b - u < tol2 )
455 d = SIGN( tol1, xm - x );
456
457 }
458
459 }
460 else
461 {
462
463 e = ( x >= xm ? a - x : b - x );
464 d = ( T )CGOLD * e;
465
466 }
467 u = ( T )( fabs( d ) >= tol1 ? x + d : x + SIGN( tol1, d ) );
468 fu = this->_func.eval( _pcom + _xicom * u );
469
470 carto::AttributedObject monitoredParams("monitored_params") ;
471 this->_func.fillMonitoredParameters( monitoredParams ) ;
472 monitoredParams.setProperty("info",
473 std::string( "mnbrak_bx") ) ;
474 this->_probe->test( _pcom + _xicom * u, monitoredParams, &fu) ;
475
476 if ( fu <= fx )
477 {
478
479 if ( u >= x )
480 a = x;
481 else
482 b = x;
483 v = w;
484 w = x;
485 x = u;
486 fv = fw;
487 fw = fx;
488 fx = fu;
489
490 }
491 else
492 {
493
494 if ( u < x )
495 a = u;
496 else
497 b = u;
498 if ( fu <= fw || w == x )
499 {
500
501 v = w;
502 w = u;
503 fv = fw;
504 fw = fu;
505
506 }
507 else if ( fu <= fv || v == x || v == w )
508 {
509
510 v = u;
511 fv = fu;
512
513 }
514
515 }
516
517 }
518 std::cerr << "brent: too many iterations" << std::endl;
519 xmin = x;
520
521 return ( T )-1;
522}
523
524
525#endif
Optimizer::_probe
OptimizerProbe< T, D > * _probe
Definition optimizer.h:91
Optimizer::_func
const ObjectiveFunc< T, D > & _func
Definition optimizer.h:88
Optimizer::Optimizer
Optimizer(const ObjectiveFunc< T, D > &func, T error, OptimizerProbe< T, D > *probe=0)
Definition optimizer.h:68
PowellOptimizer::~PowellOptimizer
virtual ~PowellOptimizer()
Definition powell.h:105
PowellOptimizer::doit
AimsVector< T, D > doit(const AimsVector< T, D > &pinit, const AimsVector< T, D > &deltaP)
Definition powell.h:116
PowellOptimizer::getCostEvolution
const std::list< T > & getCostEvolution() const
Definition powell.h:241
PowellOptimizer::PowellOptimizer
PowellOptimizer(const ObjectiveFunc< T, D > &func, T error)
Definition powell.h:93
PowellOptimizer::getOptCost
const T & getOptCost() const
Definition powell.h:231
PowellOptimizer::getDirection
const AimsVector< T, D > & getDirection(int which) const
Definition powell.h:221
carto::AttributedObject
SyntaxedObject< PropertySet > AttributedObject
ZEPS
#define ZEPS
Definition powell.h:52
TOL
#define TOL
Definition powell.h:51
ITMAX
#define ITMAX
Definition powell.h:50
CGOLD
#define CGOLD
Definition powell.h:53
GLIMIT
#define GLIMIT
Definition powell.h:47
SIGN
#define SIGN(a, b)
Definition powell.h:49
TINY
#define TINY
Definition powell.h:48
SQR
#define SQR(a)
Definition powell.h:54
GOLD
#define GOLD
Definition powell.h:46