13 #include "msdevstudio/MSconfig.h"
23 #define _GLIBCPP_USE_C99 1
56 using namespace hippodraw;
58 using namespace Numeric;
63 m_grad_cutoff( 1e-6 ),
64 m_step_cutoff( 1e-6 ),
65 m_proj_cutoff( 1e-6 ),
94 vector < double > init_params;
99 vector< double > xold =
m_xinit;
100 vector< double > xnew =
m_xinit;
101 vector< double > gk =
gradient( xold );
102 vector< double > gkPlusUn =
gradient( xnew );
104 vector< double > pk(
m_xinit.size() );
105 vector< double > s(
m_xinit.size() );
106 vector< double > y(
m_xinit.size() );
110 m_M = ( 1.0 /
norm( gk ) ) * m_M;
112 vector< vector< double > > t1 , t2;
116 double fx =
function( xnew );
127 xnew = xold + Alpha_star * pk;
128 fx =
function( xnew );
131 while( isnan( fx ) );
139 double yy =
norm( y );
140 double ss =
norm( s );
144 ( abs( Alpha_star ) < DBL_EPSILON ) ||
152 double temp = ( 1 +
innerProduct( y, m_M * y ) / ys ) / ys;
163 m_fcn -> setFreeParameters ( xold );
172 const std::vector< double >& p )
const
174 double step_fac = sqrt(2.0);
176 double phi0 =
function( x0 );
177 double dphi0 =
gradp( x0, p );
188 double Alphaim = Alpha0;
198 phii =
function( x0 + Alphai * p );
200 if ( (phii > (phi0 +
m_c1 * Alphai * dphi0)) ||
201 ( (phii >= phiim) && ( i > 1)) )
202 return zoom( x0, p, phi0, dphi0, Alphaim, Alphai );
204 dphii =
gradp( x0 + Alphai * p , p );
206 if ( abs( dphii ) <= -
m_c2 * dphi0 )
210 return zoom( x0, p, phi0, dphi0, Alphai, Alphaim );
235 const std::vector< double >& p,
236 double phi0,
double dphi0,
237 double Alphalo,
double Alphahi )
const
246 double Alpha_star = 0.0;
248 while ( !done && iter < MaxIter )
252 philo =
function( x0 + Alphalo * p );
258 phij =
function( x0 + Alphaj * p );
260 if( (phij > phi0 +
m_c1 * Alphaj * dphi0) || (phij >= philo) )
264 dphij =
gradp( x0 + Alphaj * p , p );
266 if ( abs( dphij ) <= -
m_c2 * dphi0 )
270 if ( dphij * ( Alphahi - Alphalo ) >= 0 )
280 Alpha_star = 0.5 * ( Alphahi + Alphalo );
288 const std::vector< double >& p,
290 double Alphai )
const
293 if ( Alphaim > Alphai )
294 swap( Alphaim, Alphai);
296 double phiim =
function( x0 + Alphaim * p );
297 double phii =
function( x0 + Alphai * p );
299 double dphiim =
gradp( x0 + Alphaim * p, p );
300 double dphii =
gradp( x0 + Alphai * p, p );
302 double d1 = dphiim + dphii - 3 * ( (phiim - phii)/(Alphaim - Alphai) );
303 double d2 = sqrt( d1 * d1 - dphiim * dphii);
305 double Alphaip = Alphai - (Alphai - Alphaim) *
306 ( (dphiim + d2 - d1) / (dphii - dphiim + 2 * d2) );
308 double lth = abs(Alphai - Alphaim);
310 if( abs(Alphaip - Alphai) < 0.05 * lth ||
311 abs(Alphaip - Alphaim) < 0.05 * lth ||
314 Alphaip = 0.5 * (Alphai + Alphaim);
323 vector< double > x( u.size() );
325 for(
unsigned int i = 0; i < u.size(); i++ )
328 m_fcn -> setFreeParameters ( x );
354 vector< double > x( u.size(), 0.0 );
355 vector< double > xph( u.size(), 0.0 );
357 vector< double > g( u.size() );
359 for(
unsigned int i = 0; i < u.size(); i++ )
363 m_fcn -> setFreeParameters ( x );
366 for(
unsigned int i = 0; i < u.size(); i++ )
368 for(
unsigned int j = 0; j < u.size(); j++ ) {
369 xph[j] = ( i == j ) ? ( x[j] + h ) : x[j];
371 m_fcn -> setFreeParameters ( xph );
373 g[i] = ( fxph - fx ) / h;
405 const std::vector< double > & p )
const
408 vector< double > x( u.size() );
411 for (
unsigned int i = 0; i < u.size(); i++ ) {
415 m_fcn -> setFreeParameters ( x );
418 for (
unsigned int i = 0; i < u.size(); i++ ) {
421 m_fcn -> setFreeParameters ( x );
424 return ( fxph - fx ) / h;
460 m_xinit.resize( xinit.size() );
475 for(
unsigned int i = 0; i <
m_xinit.size(); i++ )
476 cov[i].resize(
m_xinit.size(), 0.0 );
478 for(
unsigned int i = 0; i <
m_xinit.size(); i++ )
479 for(
unsigned int j = 0; j <
m_xinit.size(); j++ )
480 cov[i][j] =
m_M[i][j];
486 for(
unsigned int i = 0; i <
m_xinit.size(); i++ )
487 for(
unsigned int j = 0; j <
m_xinit.size(); j++ )
488 cov[i][j] =
m_M[i][j];
498 if( name ==
"max_iterations" )
503 map< string, double * >::const_iterator it
507 cout << name <<
" is not a valid iteration parameter name" << endl;
520 if( name ==
"max_iterations" )
529 map< string, double * >::const_iterator it
534 cout << name <<
" is not a valid iteration parameter name" << endl;