A.I.M.S algorithms


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 
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 //
60 template < class T, int D >
61 class PowellOptimizer : public Optimizer< T, D >
62 {
63 
64  public:
65  PowellOptimizer( const ObjectiveFunc< T, D >& func, T error );
66  virtual ~PowellOptimizer();
67 
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 
91 template < class T, int D >
92 inline
94  T error )
95  : Optimizer< T, D >( func, error )
96 {
97 
98  _xi = new AimsVector< T, D >[ D ];
99 
100 }
101 
102 
103 template < class T, int D >
104 inline
106 {
107 
108  delete [] _xi;
109 
110 }
111 
112 
113 template < class T, int D >
114 inline
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 
218 template < class T, int D >
219 inline
220 const AimsVector< T, D >&
222 {
223 
224  return _xi[ which ];
225 
226 }
227 
228 
229 template < class T, int D >
230 inline
232 {
233 
234  return _optCost;
235 
236 }
237 
238 
239 template < class T, int D >
240 inline
241 const std::list< T >& PowellOptimizer< T, D >::getCostEvolution() const
242 {
243 
244  return _costs;
245 
246 }
247 
248 
249 template < class T, int D >
250 inline
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 
269 template < class T, int D >
270 inline
271 void 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 
403 template < class T, int D >
404 inline
405 T 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
float max(float x, float y)
Definition: thickness.h:97
const AimsVector< T, D > & getDirection(int which) const
Definition: powell.h:221
#define ZEPS
Definition: powell.h:52
const std::list< T > & getCostEvolution() const
Definition: powell.h:241
AimsVector< T, D > doit(const AimsVector< T, D > &pinit, const AimsVector< T, D > &deltaP)
Definition: powell.h:116
PowellOptimizer(const ObjectiveFunc< T, D > &func, T error)
Definition: powell.h:93
#define TOL
Definition: powell.h:51
virtual ~PowellOptimizer()
Definition: powell.h:105
#define SIGN(a, b)
Definition: powell.h:49
#define GOLD
Definition: powell.h:46
#define SQR(a)
Definition: powell.h:54
const T & getOptCost() const
Definition: powell.h:231
#define ITMAX
Definition: powell.h:50
#define TINY
Definition: powell.h:48
#define GLIMIT
Definition: powell.h:47
#define CGOLD
Definition: powell.h:53