29 #ifndef CONSTRAINED_SOLVER_HPP
30 #define CONSTRAINED_SOLVER_HPP
43 template<
class>
class Problem = problem
78 typedef typename problem_type::f_type f_type;
79 typedef typename problem_type::g_type g_type;
80 typedef typename problem_type::h_type h_type;
104 rho_decrease ( 0.1 ),
106 rho_increase_big ( 5.0 ),
108 rho_increase_small ( 2.0 ),
110 Primal_Dual (
true ),
135 double rho_increase_big;
137 double rho_increase_small;
177 void push_x(
const double *__x );
180 void push_all(
double mu,
double E,
double E0,
182 bool truss_exit,
bool SOC,
183 double Delta,
double ared,
double pred,
double gamma );
189 void show(
bool verbose =
false )
const;
203 std::vector<double> x_hstr, mu_hstr, E_hstr, E0_hstr;
204 std::vector<int> n_CGiter_hstr;
205 std::vector<bool> truss_exit_hstr, SOC_hstr;
206 std::vector<double> Delta_hstr, ared_hstr, pred_hstr, gamma_hstr;
226 M_prob.define_problem( x_definitions );
244 return M_solver_stats;
252 void initialize_solver_with_x(
double *_x,
double *_s,
double *_lambda_h,
double *_lambda_g,
double MU );
254 void update_opt_objs(
double *_x,
double *_s,
double *_lambda_h,
255 double *_lambda_g,
double MU,
bool update_Hess );
260 void dogleg_solve(
double *_vCP_til,
double *_vN_til,
double *_s,
double Delta,
263 double model_m(
double *_s,
double *_v_til );
265 double get_vpred(
double *_s,
double *_v_til );
268 double get_normal_steps(
double *_s,
double Delta,
269 double *_dx_normal,
double *_ds_normal,
double *_v_til );
271 void get_Pr(
double *_r,
double *_Pr );
273 bool slack_feasible(
double *_w_til,
double *_v_til );
275 double get_q(
double *_v_til,
double *_w_til,
double MU );
278 double get_tangen_steps(
double *_s,
double Delta,
double *_v_til,
double MU,
279 double *_dx_tangen,
double *_ds_tangen,
280 bool &truss_exit,
int &n_CGiter );
282 double get_steps(
double *_s,
double _Delta,
double MU,
283 double *_dx,
double *_ds,
double &_vpred,
double &_q,
284 bool &truss_exit,
int &n_CGiter );
286 double get_phi(
double *_x,
double *_s,
double MU,
double nu,
double *rhs );
287 double get_phi(
double *_x,
double *_s,
double MU,
double nu );
288 double get_phi(
double *_s,
double MU,
double nu );
290 double get_SOC(
double *_s,
double *_rhs_xs_p_d,
double *_xy,
double *_sy );
297 void evaluate_Ag(
double *x );
299 void evaluate_Ah(
double *x );
301 void evaluate_A_hat(
double *_x,
double *_s );
303 void evaluate_Hessian_L(
double *_x,
double *_lambda_h,
double *_lambda_g );
305 void evaluate_G(
double *_x,
double *_s,
double *_lambda_h,
double *_lambda_g,
306 double MU,
bool evaluate__Hessians_flag );
307 double evaluate_E_current(
double *_s,
double *_lambda_h,
double *_lambda_g,
double MU );
309 void evaluate_gradf(
double *_x );
311 void evaluate_fxgxhx(
double *_x );
322 void print_lambda_LS(
double *_lambda_g,
double *_lambda_h );
328 void print_Hessian_L();
343 double fx, gradf[
_E_n ];
353 static const double Max_Allowed_Delta = 1e8;
364 template<
typename Data,
template<
class>
class Problem>
367 M_prob( x_definitions ),
368 M_solver_stats( this, true ),
379 std::cerr <<
"SolverConstrained::SolverConstrained()\n";
383 template<
typename Data,
template<
class>
class Problem>
388 template<
typename Data,
template<
class>
class Problem>
393 std::ifstream fin(
"N_options.inp" );
398 options.Primal_Dual =
true;
404 if ( !strcmp( str,
"Delta_init" ) ) options.Delta_init = num;
406 if ( !strcmp( str,
"eta" ) ) options.eta = num;
408 if ( !strcmp( str,
"tau" ) ) options.tau = num;
410 if ( !strcmp( str,
"theta_N" ) ) options.theta_N = num;
412 if ( !strcmp( str,
"zeta_N" ) ) options.zeta_N = num;
414 if ( !strcmp( str,
"eTOL" ) ) options.eTOL = num;
416 if ( !strcmp( str,
"rho_N" ) ) options.rho_N = num;
418 if ( !strcmp( str,
"s_theta_init" ) ) options.s_theta_init = num;
420 if ( !strcmp( str,
"eMU_init" ) ) options.eMU_init = num;
422 if ( !strcmp( str,
"MU_init" ) ) options.MU_init = num;
424 if ( !strcmp( str,
"nu0_init" ) ) options.nu0_init = num;
426 if ( !strcmp( str,
"Primal_Dual" ) ) options.Primal_Dual = (
bool )num;
428 if ( !strcmp( str,
"max_TR_iter" ) ) options.max_TR_iter = (
int )num;
430 if ( !strcmp( str,
"max_hom_iter" ) ) options.max_hom_iter = (
int )num;
432 if ( !strcmp( str,
"rho_decrease" ) ) options.rho_decrease = num;
434 if ( !strcmp( str,
"rho_big" ) ) options.rho_big = num;
436 if ( !strcmp( str,
"rho_increase_big" ) ) options.rho_increase_big = num;
438 if ( !strcmp( str,
"rho_small" ) ) options.rho_small = num;
440 if ( !strcmp( str,
"rho_increase_small" ) ) options.rho_increase_small = num;
447 template<
typename Data,
template<
class>
class Problem>
448 void SolverConstrained<Data,Problem>::print_Hessian_L()
452 for (
int i = 0; i < n; i++ )
456 for (
int j = 0; j < n; j++ )
457 fprintf( stderr,
"%10.3f", Hessian_L[i][j] );
463 template<
typename Data,
template<
class>
class Problem>
464 void SolverConstrained<Data,Problem>::evaluate_Hessian_L(
double *_x,
double *_lambda_h,
double *_lambda_g )
469 M_prob.evaluate ( _x, __fx, diff_order<2>() );
470 M_prob.evaluate ( _x, __gx, diff_order<2>() );
471 M_prob.evaluate ( _x, __hx, diff_order<2>() );
473 for (
int i = 0; i < _E_n; i++ )
476 for (
int j = i; j < _E_n; j++ )
480 for (
int __k = 0; __k < _E_f; ++__k )
482 val += __fx.hessian( __k, i, j );
485 for (
int m = 0; m < _E_h; m++ )
487 val += _lambda_h[ m ] * __hx.hessian( m,i,j );
490 for (
int m = 0; m < _E_g; m++ )
492 val += _lambda_g[ m ] * __gx.hessian( m,i,j );
495 Hessian_L[ i ][ j ] =
val;
496 Hessian_L[ j ][ i ] =
val;
501 template<
typename Data,
template<
class>
class Problem>
502 void SolverConstrained<Data,Problem>::print_G()
504 int n = M_prob.n(), nv = n+_E_g+2*n;
506 for (
int i = 0; i < nv; i++ )
510 for (
int j = 0; j < nv; j++ )
511 fprintf( stderr,
"%12.5f", G[i][j] );
517 template<
typename Data,
template<
class>
class Problem>
518 void SolverConstrained<Data,Problem>::evaluate_G(
double *_x,
double *_s,
double *_lambda_h,
519 double *_lambda_g,
double MU,
bool evaluate_Hessians_flag )
524 if ( evaluate_Hessians_flag )
525 this->evaluate_Hessian_L( _x, _lambda_h, _lambda_g );
527 for (
int i = 0; i < nv; i++ )
528 for (
int j = 0; j < nv; j++ )
531 for (
int i = 0; i < n; i++ )
532 for (
int j = 0; j < n; j++ )
533 G[ i ][ j ] = Hessian_L[ i ][ j ];
535 for (
int i = 0; i < _E_g+2*n; i++ )
536 if ( options.Primal_Dual && _lambda_g[ i ] > 0 )
537 G[ n+i ][ n+i ] = _lambda_g[ i ] * _s[ i ];
540 G[ n+i ][ n+i ] = MU;
558 template<
typename Data,
template<
class>
class Problem>
559 void SolverConstrained<Data,Problem>::evaluate_Ag(
double *_x )
562 M_prob.evaluate ( _x, __gx, diff_order<1>() );
566 assert ( n == _E_n );
569 for (
int i = 0; i < _E_n; i++ )
570 for (
int m = 0; m < _E_g+2*_E_n; m++ )
573 for (
int i = 0; i < _E_n; i++ )
575 for (
int m = 0; m < _E_g; m++ )
577 Ag[ i ][ m ] = __gx.gradient( m,i );
580 int __index = _E_g+i;
581 Ag[ i ][ __index ] = -1.0;
582 Ag[ i ][ __index + n] = +1.0;
586 template<
typename Data,
template<
class>
class Problem>
587 void SolverConstrained<Data,Problem>::evaluate_Ah(
double *_x )
590 M_prob.evaluate ( _x, __hx, diff_order<1>() );
595 for (
int i = 0; i < n; i++ )
596 for (
int m = 0; m < _E_h; m++ )
599 for (
int i = 0; i < n; i++ )
600 for (
int m = 0; m < _E_h; m++ )
601 Ah[ i ][ m ] = __hx.gradient( m,i );
605 template<
typename Data,
template<
class>
class Problem>
606 void SolverConstrained<Data,Problem>::print_Ag()
610 for (
int i = 0; i < n; i++ )
614 for (
int j = 0; j < _E_g+2*n; j++ )
615 fprintf( stderr,
"%10.1f", Ag[i][j] );
621 template<
typename Data,
template<
class>
class Problem>
622 void SolverConstrained<Data,Problem>::print_Ah()
626 for (
int i = 0; i < n; i++ )
630 for (
int j = 0; j < _E_h; j++ )
631 fprintf( stderr,
"%10.1f", Ah[i][j] );
637 template<
typename Data,
template<
class>
class Problem>
638 void SolverConstrained<Data,Problem>::evaluate_fxgxhx(
double *_x )
645 M_prob.evaluate ( _x, __fx, diff_order<0>() );
646 M_prob.evaluate ( _x, __gx, diff_order<0>() );
647 M_prob.evaluate ( _x, __hx, diff_order<0>() );
649 fx = __fx.value( 0 );
651 for ( uint __i = 0; __i < __gx.value().size(); ++__i )
652 gx[__i] = __gx.value( __i );
654 for ( uint __i = 0; __i < _E_h; ++__i )
655 hx[__i] = __hx.value( __i );
659 template<
typename Data,
template<
class>
class Problem>
661 SolverConstrained<Data,Problem>::initialize_solver_with_x(
double *_x,
663 double *_lambda_h,
double *_lambda_g,
667 M_prob.evaluate ( _x, __gx, diff_order<0>() );
669 for ( uint __i = 0; __i < __gx.value().size(); ++__i )
671 gx[__i] = __gx.value( __i );
672 _s[ __i ] = std::max( -gx[ __i ], options.s_theta_init );
680 this->update_opt_objs( _x, _s, _lambda_h, _lambda_g, MU,
true );
684 template<
typename Data,
template<
class>
class Problem>
686 SolverConstrained<Data,Problem>::evaluate_A_hat(
double *_x,
double *_s )
690 this->evaluate_Ag( _x );
691 this->evaluate_Ah( _x );
694 for (
int i = 0; i < n+_E_g+2*n; i++ )
695 for (
int m = 0; m < _E_h+_E_g+2*n; m++ )
698 for (
int i = 0; i < n; i++ )
700 for (
int m = 0; m < _E_h; m++ )
701 A_hat[ i ][ m ] = Ah[ i ][ m ];
703 for (
int m = 0; m < _E_g+2*n; m++ )
704 A_hat[ i ][ _E_h+m ] = Ag[ i ][ m ];
707 for (
int m = 0; m < _E_g+2*n; m++ )
708 A_hat[ n+m ][ _E_h+m ] = _s[ m ];
712 template<
typename Data,
template<
class>
class Problem>
713 void SolverConstrained<Data,Problem>::print_A_hat()
717 for (
int i = 0; i < n+_E_g+2*n; i++ )
721 for (
int j = 0; j < _E_h+_E_g+2*n; j++ )
722 fprintf( stderr,
"%14.9f", A_hat[i][j] );
728 template<
typename Data,
template<
class>
class Problem>
729 void SolverConstrained<Data,Problem>::solve_AtA_system(
double _r[ _E_n + _E_g+2*_E_n ],
double _b[ _E_h + _E_g+2*_E_n ],
730 double _g[ _E_n + _E_g+2*_E_n],
double _l[ _E_h + _E_g+2*_E_n ] )
734 int n_A_hat = n + _E_g+2*n;
735 int m_A_hat = _E_h + _E_g+2*n;
736 int n_big = n_A_hat + m_A_hat;
737 int big_N = _E_n + _E_g+2*_E_n + _E_h + _E_g+2*_E_n;
739 double bigA[ big_N*big_N ], rhs[ big_N ];
741 for (
int i = 0; i < n_big*n_big; i++ )
744 for (
int i = 0; i < n_A_hat; i++ )
746 bigA[ i + i*n_big ] = 1;
752 for (
int i = 0; i < n_A_hat; i++ )
754 for (
int j = 0; j < m_A_hat; j++ )
756 val = A_hat[ i ][ j ];
757 bigA[ ( i ) + ( n_A_hat+j )*n_big ] =
val;
758 bigA[ ( n_A_hat+j ) + ( i )*n_big ] =
val;
762 for (
int j = 0; j < m_A_hat; j++ )
763 rhs[ n_A_hat + j ] = -_b[ j ];
781 #warning disabled Lusolve
784 for (
int i = 0; i < n_A_hat; i++ )
787 for (
int i = 0; i < m_A_hat; i++ )
788 _l[ i ] = sol[ n_A_hat+i ];
792 template<
typename Data,
template<
class>
class Problem>
793 void SolverConstrained<Data,Problem>::evaluate_gradf(
double *_x )
796 M_prob.evaluate ( _x, __fx, diff_order<1>() );
798 for ( uint __i = 0; __i < _E_n; ++__i )
800 gradf[__i] = __fx.gradient( 0, __i );
804 template<
typename Data,
template<
class>
class Problem>
805 void SolverConstrained<Data,Problem>::update_opt_objs(
double *_x,
double *_s,
806 double *_lambda_h,
double *_lambda_g,
807 double MU,
bool update_Hess )
810 double _r[ _E_n + _E_g+2*_E_n ], _b[ _E_h + _E_g+2*_E_n ],
811 _dummy[ _E_n + _E_g+2*_E_n ];
812 double _lambda_LS[ _E_h + _E_g+2*_E_n ];
814 this->evaluate_fxgxhx( _x );
815 this->evaluate_gradf( _x );
817 for (
int i = 0; i < n; i++ )
818 _r[ i ] = -gradf[ i ];
820 for (
int i = n; i < n+_E_g+2*n; i++ )
823 for (
int i = 0; i < _E_h + _E_g+2*n; i++ )
826 this->evaluate_A_hat( _x, _s );
828 this->solve_AtA_system( _r, _b, _dummy, _lambda_LS );
830 for (
int i = 0; i < _E_h; i++ )
831 _lambda_h[ i ] = _lambda_LS[ i ];
833 for (
int i = 0; i < _E_g+2*n; i++ )
834 _lambda_g[ i ] = _lambda_LS[ _E_h+i ];
836 this->evaluate_G( _x, _s, _lambda_h, _lambda_g, MU, update_Hess );
847 template<
typename Data,
template<
class>
class Problem>
848 void SolverConstrained<Data,Problem>::print_lambda_LS(
double *_lambda_h,
double *_lambda_g )
851 std::cerr <<
"\nlambda_h = ";
852 __print_vec( _E_h, _lambda_h );
853 std::cerr <<
"\nlambda_g = ";
854 __print_vec( _E_g+2*n, _lambda_g );
858 template<
typename Data,
template<
class>
class Problem>
859 double SolverConstrained<Data,Problem>::evaluate_E_current(
double *_s,
double *_lambda_h,
double *_lambda_g,
double MU )
862 double Ahlh[ _E_n ], Aglg[ _E_n ], dfAhlhAglg[ _E_n ], Slg_mue[ _E_g+2*_E_n], gps[ _E_g+2*_E_n];
865 for (
int i = 0; i < n; i++ )
869 for (
int m = 0; m < _E_h; m++ )
870 Ahlh[ i ] += Ah[ i ][ m ] * _lambda_h[ m ];
874 for (
int m = 0; m < _E_g+2*n; m++ )
875 Aglg[ i ] += Ag[ i ][ m ] * _lambda_g[ m ];
877 dfAhlhAglg[ i ] = gradf[ i ] + Ahlh[ i ] + Aglg[ i ];
880 E = std::max( E, __Norm_infty( n, dfAhlhAglg ) );
882 for (
int m = 0; m < _E_g+2*n; m++ )
883 Slg_mue[ m ] = _s[ m ] * _lambda_g[ m ] - MU;
885 E = std::max( E, __Norm_infty( _E_g+2*n, Slg_mue ) );
887 E = std::max( E, __Norm_infty( _E_h, hx ) );
889 for (
int m = 0; m < _E_g+2*n; m++ )
890 gps[ m ] = gx[ m ] + _s[ m ];
892 E = std::max( E, __Norm_infty( _E_g+2*n, gps ) );
913 template<
typename Data,
template<
class>
class Problem>
914 double SolverConstrained<Data,Problem>::model_m(
double *_s,
double *_v_til )
916 int n = M_prob.n(), nv = n+_E_g+2*n;
917 double A_hatTxv_til[ _E_h + _E_g+2*_E_n ], rhs[ _E_h + _E_g+2*_E_n ];
919 for (
int i = 0; i < _E_h+_E_g+2*n; i++ )
921 A_hatTxv_til[ i ] = 0;
923 for (
int j = 0; j < nv; j++ )
924 A_hatTxv_til[ i ] += A_hat[ j ][ i ] * _v_til[ j ];
927 for (
int m = 0; m < _E_h; m++ )
930 for (
int m = 0; m < _E_g+2*n; m++ )
931 rhs[ _E_h + m ] = gx[ m ] + _s[ m ];
933 double model = 2. * __Dot( _E_h+_E_g+2*n, rhs, A_hatTxv_til )
934 + __Dot( _E_h+_E_g+2*n, A_hatTxv_til, A_hatTxv_til );
948 template<
typename Data,
template<
class>
class Problem>
949 void SolverConstrained<Data,Problem>::dogleg_solve(
double *_vCP_til,
double *_vN_til,
double *_s,
double _Delta,
952 int n = M_prob.n(), nv = n+_E_g+2*n;
953 double v_diff[ _E_n+_E_g+2*_E_n ];
954 double zeta_N = options.zeta_N, tau = options.tau;
955 double theta1 = 1, theta2 = 1, theta3 = 1;
957 theta1 = std::min( theta1, zeta_N * _Delta / __Norm_2( nv, _vN_til ) );
959 for (
int i = n; i < nv; i++ )
960 if ( _vN_til[ i ] < 0 )
961 theta1 = std::min( theta1, -tau/2./_vN_til[ i ] );
963 assert( theta1 > 0 );
966 for (
int i = 0; i < nv; i++ )
967 _v_til[ i ] = theta1*_vN_til[ i ];
971 for (
int i = 0; i < nv; i++ )
972 v_diff[ i ] = _vN_til[ i ] - _vCP_til[ i ];
974 double a = __Dot( nv, v_diff, v_diff );
975 double b = 2. * __Dot( nv, v_diff, _vCP_til );
976 double c = __Dot( nv, _vCP_til, _vCP_til ) - ( zeta_N*_Delta )*( zeta_N*_Delta );
977 double root = b*b - 4.*a*c;
982 double r2 = ( ( -b + std::sqrt( root ) ) / ( 2.*a ) );
997 for (
int i = n; i < nv; i++ )
998 if ( v_diff[ i ] < 0 )
999 theta2 = std::min( theta2, ( -tau/2. - _vCP_til[ i ] ) / ( v_diff[ i ] ) );
1007 theta3 = std::min( theta3, zeta_N * _Delta / __Norm_2( nv, _vCP_til ) );
1009 for (
int i = n; i < nv; i++ )
1010 if ( _vCP_til[ i ] < 0 )
1011 theta3 = std::min( theta3, -tau/2./_vCP_til[ i ] );
1013 assert( theta3 > 0 );
1015 for (
int i = 0; i < nv; i++ )
1016 _v_til[ i ] = theta3 * _vCP_til[ i ];
1020 for (
int i = 0; i < nv; i++ )
1021 _v_til[ i ] = ( 1-theta2 ) * _vCP_til[ i ] + theta2 * _vN_til[ i ];
1023 double theta1xvN_til[ _E_n+_E_g+2*_E_n ];
1025 for (
int i = 0; i < nv; i++ )
1026 theta1xvN_til[ i ] = theta1 * _vN_til[ i ];
1028 if ( this->model_m( _s, _v_til ) >= this->model_m( _s, theta1xvN_til ) )
1029 for (
int i = 0; i < nv; i++ )
1030 _v_til[ i ] = theta1xvN_til[ i ];
1046 template<
typename Data,
template<
class>
class Problem>
1047 double SolverConstrained<Data,Problem>::get_vpred(
double *_rhs,
double *_v_til )
1049 int n = M_prob.n(), nv = n + _E_g+2*n;
1050 double A_hatTxv_til[ _E_h+_E_g+2*_E_n ], rhspAv[ _E_h + _E_g+2*_E_n ];
1052 for (
int i = 0; i < _E_h+_E_g+2*n; i++ )
1054 A_hatTxv_til[ i ] = 0;
1056 for (
int j = 0; j < nv; j++ )
1057 A_hatTxv_til[ i ] += A_hat[ j ][ i ] * _v_til[ j ];
1059 rhspAv[ i ] = _rhs[ i ] + A_hatTxv_til[ i ];
1062 return __Norm_2( _E_h+_E_g+2*n, _rhs ) - __Norm_2( _E_h+_E_g+2*n, rhspAv );
1065 template<
typename Data,
template<
class>
class Problem>
1066 double SolverConstrained<Data,Problem>::get_normal_steps(
double *_s,
double _Delta,
1067 double *_dx_normal,
double *_ds_normal,
1070 int n = M_prob.n(), nv = n + _E_g+2*n;
1071 double vCP_til[ _E_n + _E_g+2*_E_n ], vN_til[ _E_n + _E_g+2*_E_n ],
1072 rhs[ _E_h + _E_g+2*_E_n ], A_hatxrhs[ _E_n + _E_g+2*_E_n ],
1073 A_hatTxA_hat[ _E_h + _E_g+2*_E_n ][ _E_h + _E_g+2*_E_n ],
1074 A_hatTxA_hat2[ _E_h + _E_g+2*_E_n ][ _E_h + _E_g+2*_E_n ],
1075 A_hat2xrhs [ _E_h+_E_g+2*_E_n ];
1076 double zeros[ _E_n + _E_g+2*_E_n ], dummy[ _E_n + _E_g+2*_E_n ], x_til[ _E_h + _E_g+2*_E_n ];
1078 for (
int m = 0; m < _E_h; m++ )
1081 for (
int m = 0; m < _E_g+2*n; m++ )
1082 rhs[ _E_h + m ] = gx[ m ] + _s[ m ];
1085 for (
int mi = 0; mi < _E_h+_E_g+2*n; mi++ )
1087 for (
int mj = 0; mj < _E_h+_E_g+2*n; mj++ )
1089 A_hatTxA_hat[ mi ][ mj ] = 0;
1091 for (
int i = 0; i < nv; i++ )
1092 A_hatTxA_hat[ mi ][ mj ] += A_hat[ i ][ mi ] * A_hat[ i ][ mj ];
1097 for (
int mi = 0; mi < _E_h+_E_g+2*n; mi++ )
1099 for (
int mj = 0; mj < _E_h+_E_g+2*n; mj++ )
1101 A_hatTxA_hat2[ mi ][ mj ] = 0;
1103 for (
int i = 0; i < _E_h+_E_g+2*n; i++ )
1104 A_hatTxA_hat2[ mi ][ mj ] += A_hatTxA_hat[ i ][ mi ] * A_hatTxA_hat[ i ][ mj ];
1108 for (
int i = 0; i < nv; i++ )
1113 for (
int m = 0; m < _E_h+_E_g+2*n; m++ )
1114 A_hatxrhs[ i ] += A_hat[ i ][ m ] * rhs[ m ];
1117 for (
int i = 0; i < _E_h+_E_g+2*n; i++ )
1119 A_hat2xrhs[ i ] = 0;
1121 for (
int m = 0; m < _E_h+_E_g+2*n; m++ )
1122 A_hat2xrhs[ i ] += A_hatTxA_hat2[ i ][ m ] * rhs[ m ];
1125 double norm = __Norm_2( nv, A_hatxrhs );
1126 double alpha = norm*norm / __Dot( _E_h+_E_g+2*n, rhs, A_hat2xrhs );
1128 for (
int i = 0; i < nv; i++ )
1129 vCP_til[ i ] = -alpha * A_hatxrhs[ i ];
1131 this->solve_AtA_system( zeros, rhs, dummy, x_til );
1133 for (
int i = 0; i < n+_E_g+2*n; i++ )
1137 for (
int m = 0; m < _E_h+_E_g+2*n; m++ )
1138 vN_til[ i ] += -A_hat[ i ][ m ] * x_til[ m ];
1141 this->dogleg_solve( vCP_til, vN_til, _s, _Delta, _v_til );
1143 for (
int i = 0; i < n; i++ )
1144 _dx_normal[ i ] = _v_til[ i ];
1146 for (
int i = 0; i < _E_g+2*n; i++ )
1147 _ds_normal[ i ] = _s[ i ] * _v_til[ n+i ];
1183 return this->get_vpred( rhs, _v_til );
1187 template<
typename Data,
template<
class>
class Problem>
1188 double SolverConstrained<Data,Problem>::get_q(
double *_v_til,
double *_w_til,
double MU )
1190 int n = M_prob.n(), nd = n + _E_g+2*n;
1191 double d_til[ _E_n+_E_g+2*_E_n ], gr[ _E_n+_E_g+2*_E_n ], Gd_til[ _E_n+_E_g+2*_E_n ];
1193 __SumVecs( nd, _v_til, _w_til, d_til );
1195 for (
int i = 0; i < n; i++ )
1196 gr[ i ] = gradf[ i ];
1198 for (
int i = n; i < nd; i++ )
1201 for (
int i = 0; i < nd; i++ )
1205 for (
int j = 0; j < nd; j++ )
1206 Gd_til[ i ] += G[ i ][ j ] * d_til[ j ];
1209 return __Dot( nd, gr, d_til ) + 0.5 * __Dot( nd, d_til, Gd_til );
1212 template<
typename Data,
template<
class>
class Problem>
1213 void SolverConstrained<Data,Problem>::get_Pr(
double *_r,
double *_Pr )
1215 int n = M_prob.n(), nw = n + _E_g+2*n;
1216 double zeros[ _E_n + _E_g+2*_E_n ], dummy[ _E_n + _E_g+2*_E_n ],
1217 Atr[ _E_h + _E_g+2*_E_n ], x_Atr[ _E_h + _E_g+2*_E_n ];
1219 for (
int i = 0; i < nw; i++ )
1224 for (
int i = 0; i < _E_h+_E_g+2*n; i++ )
1228 for (
int j = 0; j < nw; j++ )
1229 Atr[ i ] += A_hat[ j ][ i ] * _r[ j ];
1232 this->solve_AtA_system( zeros, Atr, dummy, x_Atr );
1234 for (
int i = 0; i < nw; i++ )
1238 for (
int m = 0; m < _E_h+_E_g+2*n; m++ )
1239 _Pr[ i ] += -A_hat[ i ][ m ] * x_Atr[ m ];
1243 template<
typename Data,
template<
class>
class Problem>
1244 bool SolverConstrained<Data,Problem>::slack_feasible(
double *_w_til,
double *_v_til )
1247 double tau = options.tau;
1248 bool feasible =
true;
1250 for (
int i = n; i < n+_E_g+2*n; i++ )
1252 if ( _w_til[ i ] + _v_til[ i ] + tau < 1e-13 )
1261 template<
typename Data,
template<
class>
class Problem>
1262 double SolverConstrained<Data,Problem>::get_tangen_steps(
double *_s,
double _Delta,
double *_v_til,
double MU,
1263 double *_dx_tangen,
double *_ds_tangen,
1264 bool &truss_exit,
int &n_CGiter )
1267 int nw = n + _E_g+2*n;
1268 double last_slack_feasible_rg;
1269 double w_til[ _E_n+_E_g+2*_E_n ], wP_til[ _E_n+_E_g+2*_E_n ],
1270 last_slack_feasible_w_til[ _E_n+_E_g+2*_E_n ],
1271 r[ _E_n+_E_g+2*_E_n ], rP[ _E_n+_E_g+2*_E_n ], Gp[ _E_n+_E_g+2*_E_n ],
1272 g[ _E_n+_E_g+2*_E_n ], gP[ _E_n+_E_g+2*_E_n ],
1273 p[ _E_n+_E_g+2*_E_n ], pP[ _E_n+_E_g+2*_E_n ],
1274 last_slack_feasible_p[ _E_n+_E_g+2*_E_n ];
1279 for (
int i = 0; i < nw; i++ )
1283 for (
int j = 0; j < nw; j++ )
1284 Gp[ i ] += G[ i ][ j ] * _v_til[ j ];
1287 for (
int i = 0; i < n; i++ )
1288 r[ i ] = gradf[ i ] + Gp[ i ];
1290 for (
int i = n; i < nw; i++ )
1291 r[ i ] = -MU + Gp[ i ];
1302 for (
int i = 0; i < nw; i++ )
1306 last_slack_feasible_w_til[ i ] = 0;
1309 last_slack_feasible_p[ i ] = -g[ i ];
1312 last_slack_feasible_rg = __Dot( nw, g, r );
1314 double tol = 0.01 * std::sqrt( __Dot( nw, g, r ) );
1317 double gamma, alpha, beta;
1319 while ( iter < 2*( n-_E_h ) && !done )
1322 std::cerr <<
"\nCGiter = " << iter <<
",\t r'g = " << __Dot( nw, g, r );
1325 for (
int i = 0; i < nw; i++ )
1329 for (
int j = 0; j < nw; j++ )
1330 Gp[ i ] += G[ i ][ j ] * p[ j ];
1333 gamma = __Dot( nw, p, Gp );
1337 double a = __Dot( nw, p, p );
1338 double b = 2. * __Dot( nw, w_til, p );
1339 double c = __Dot( nw, w_til, w_til ) + __Dot( nw, _v_til, _v_til ) - _Delta*_Delta;
1340 double thetaa = ( -b + std::sqrt( b*b - 4.*a*c ) ) / ( 2.*a );
1341 assert( thetaa > 0 );
1343 for (
int i = 0; i < nw; i++ )
1344 wP_til[ i ] = w_til[ i ] + thetaa * p[ i ];
1349 std::cerr <<
"\nCG indefinitness EXIT: ||wP_til||^2 + ||v_til||^2 = "
1350 << __Dot( nw, wP_til, wP_til ) + __Dot( nw, _v_til, _v_til )
1351 <<
" ,\t Delta^2 = "
1352 << _Delta*_Delta <<
"\n";
1358 alpha = __Dot( nw, r, g ) / gamma;
1360 for (
int i = 0; i < nw; i++ )
1361 wP_til[ i ] = w_til[ i ] + alpha * p[ i ];
1363 if ( __Dot( nw, wP_til, wP_til ) + __Dot( nw, _v_til, _v_til ) > _Delta*_Delta )
1365 double a = __Dot( nw, p, p );
1366 double b = 2. * __Dot( nw, w_til, p );
1367 double c = __Dot( nw, w_til, w_til ) + __Dot( nw, _v_til, _v_til ) - _Delta*_Delta;
1368 double thetab = ( -b + std::sqrt( b*b - 4.*a*c ) ) / ( 2.*a );
1369 assert( thetab > 0 );
1371 for (
int i = 0; i < nw; i++ )
1372 wP_til[ i ] = w_til[ i ] + thetab * p[ i ];
1377 std::cerr <<
"\nCG truss EXIT: ||wP_til||^2 + ||v_til||^2 = "
1378 << __Dot( nw, wP_til, wP_til ) + __Dot( nw, _v_til, _v_til )
1379 <<
" ,\t Delta^2 = "
1380 << _Delta*_Delta <<
"\n";
1387 for (
int i = 0; i < nw; i++ )
1388 rP[ i ] = r[ i ] + alpha * Gp[ i ];
1397 if ( __Dot( nw, gP, rP ) < tol )
1401 std::cerr <<
"\nCG tol EXIT: \trP'gP = "
1402 << __Dot( nw, gP, rP )
1411 beta = __Dot( nw, gP, rP ) / __Dot( nw, g, r );
1413 for (
int i = 0; i < nw; i++ )
1415 pP[ i ] = -gP[ i ] + beta * p[ i ];
1416 w_til[ i ] = wP_til[ i ];
1422 if ( this->slack_feasible( wP_til, _v_til ) )
1424 for (
int i = 0; i < nw; i++ )
1426 last_slack_feasible_w_til[ i ] = wP_til[ i ];
1427 last_slack_feasible_p[ i ] = p[ i ];
1430 last_slack_feasible_rg = __Dot( nw, g, r );
1437 if ( !( this->slack_feasible( wP_til, _v_til ) ) )
1439 double a = __Dot( nw, last_slack_feasible_p, last_slack_feasible_p );
1440 double b = 2. * __Dot( nw, last_slack_feasible_w_til, last_slack_feasible_p );
1441 double c = __Dot( nw, last_slack_feasible_w_til, last_slack_feasible_w_til )
1442 + __Dot( nw, _v_til, _v_til ) - _Delta*_Delta;
1443 double thetac = ( -b + std::sqrt( b*b - 4.*a*c ) ) / ( 2.*a );
1444 assert( thetac > 0 );
1445 double tau = options.tau;
1448 for (
int i = n; i < nw; i++ )
1449 if ( last_slack_feasible_p[ i ] < 0 )
1451 thetac = std::min( thetac,
1452 ( -tau - _v_til[ i ] - last_slack_feasible_w_til[ i ] )
1453 / last_slack_feasible_p[ i ] );
1456 assert( thetac > 0 );
1458 for (
int i = 0; i < nw; i++ )
1459 wP_til[ i ] = last_slack_feasible_w_til[ i ] + thetac * last_slack_feasible_p[ i ];
1462 std::cerr <<
"\nCorrected slack feasibility violation of wP_til: tau = " << tau <<
", wP_til_s + v_til_s = " ;
1464 for (
int i = n; i < nw; i++ )
1465 std::cerr <<
"\n " << wP_til[ i ] + _v_til[ i ];
1467 std::cerr<<
"\nr'g = " << last_slack_feasible_rg <<
"\n";
1473 for (
int i = 0; i < n; i++ )
1474 _dx_tangen[ i ] = wP_til[ i ];
1476 for (
int i = 0; i < _E_g+2*n; i++ )
1477 _ds_tangen[ i ] = _s[ i ] * wP_til[ n+i ];
1507 return this->get_q( _v_til, wP_til, MU );
1511 template<
typename Data,
template<
class>
class Problem>
1512 double SolverConstrained<Data,Problem>::get_steps(
double *_s,
double _Delta,
double MU,
1513 double *_dx,
double *_ds,
double &_vpred,
double &_q,
1514 bool &truss_exit,
int &n_CGiter )
1517 double dx_normal[ _E_n ], ds_normal[ _E_g+2*_E_n ], v_til[ _E_n+_E_g+2*_E_n ];
1518 double dx_tangen[ _E_n ], ds_tangen[ _E_g+2*_E_n ], d[ _E_n+_E_g+2*_E_n ];
1520 _vpred = this->get_normal_steps( _s, _Delta, dx_normal, ds_normal, v_til );
1521 _q = this->get_tangen_steps( _s, _Delta, v_til, MU, dx_tangen, ds_tangen, truss_exit, n_CGiter );
1523 __SumVecs( n, dx_normal, dx_tangen, _dx );
1524 __SumVecs( _E_g+2*n, ds_normal, ds_tangen, _ds );
1538 for (
int i = 0; i < n; i++ )
1541 for (
int i = 0; i < _E_g+2*n; i++ )
1542 d[ n+i ] = _ds[ i ];
1544 return __Norm_2( n+_E_g+2*n, d );
1547 template<
typename Data,
template<
class>
class Problem>
1548 double SolverConstrained<Data,Problem>::get_phi(
double *_x,
double *_s,
double MU,
1549 double nu,
double *rhs )
1556 double f_at_x, g_at_x[ _E_g+2*_E_n ], h_at_x[ _E_h ];
1559 M_prob.fn_double( _x, f_at_x );
1560 M_prob.gn_double( _x, g_at_x );
1562 for (
int i = 0; i < n; i++ )
1564 g_at_x[ _E_g+i ] = M_prob.x_l( i ) - _x[i];
1565 g_at_x[ _E_g+n+i ] = _x[i] - M_prob.x_u( i );
1568 M_prob.hn_double( _x, h_at_x );
1570 for (
int m = 0; m < _E_h; m++ )
1571 rhs[ m ] = h_at_x[ m ];
1573 for (
int m = 0; m < _E_g+2*n; m++ )
1574 rhs[ _E_h + m ] = g_at_x[ m ] + _s[ m ];
1578 for (
int i = 0; i < _E_g+2*n; i++ )
1579 phi_val += -MU * std::log( _s[ i ] );
1581 phi_val += nu * __Norm_2( _E_g+2*n, rhs );
1588 M_prob.evaluate ( _x, __fx, diff_order<0>() );
1589 M_prob.evaluate ( _x, __gx, diff_order<0>() );
1590 M_prob.evaluate ( _x, __hx, diff_order<0>() );
1592 double __f_at_x = __fx.value( 0 );
1593 double __g_at_x[ _E_g+2*_E_n ];
1595 for ( uint __i = 0; __i < __gx.value().size(); ++__i )
1596 __g_at_x[__i] = __gx.value( __i );
1598 double __h_at_x[ _E_h ];
1600 for ( uint __i = 0; __i < _E_h; ++__i )
1601 __h_at_x[__i] = __hx.value( __i );
1603 for (
int m = 0; m < _E_h; m++ )
1604 rhs[ m ] = __h_at_x[ m ];
1606 for (
int m = 0; m < _E_g+2*_E_n; m++ )
1607 rhs[ _E_h + m ] = __g_at_x[ m ] + _s[ m ];
1611 for (
int i = 0; i < _E_g+2*_E_n; i++ )
1612 phi_val += -MU * std::log( _s[ i ] );
1614 phi_val += nu * __Norm_2( _E_g+2*_E_n, rhs );
1619 template<
typename Data,
template<
class>
class Problem>
1620 double SolverConstrained<Data,Problem>::get_phi(
double *_x,
double *_s,
double MU,
double nu )
1622 double rhs[ _E_h + _E_g+2*_E_n ];
1623 double phi_val = this->get_phi( _x, _s, MU, nu, rhs );
1628 template<
typename Data,
template<
class>
class Problem>
1629 double SolverConstrained<Data,Problem>::get_SOC(
double *_s,
double *_rhs_xs_pd,
1630 double *_xy,
double *_sy )
1635 double zeros[ _E_n + _E_g+2*_E_n ], dummy[ _E_n + _E_g+2*_E_n ],
1636 x_til[ _E_h + _E_g+2*_E_n ], y[ _E_n + _E_g+2*_E_n ];
1638 for (
int i = 0; i < n+_E_g+2*_E_n; i++ )
1641 this->solve_AtA_system( zeros, _rhs_xs_pd, dummy, x_til );
1643 for (
int i = 0; i < n+_E_g+2*n; i++ )
1647 for (
int m = 0; m < _E_h+_E_g+2*n; m++ )
1648 y[ i ] += A_hat[ i ][ m ] * x_til[ m ];
1651 for (
int i = 0; i < n; i++ )
1654 for (
int i = 0; i < _E_g+2*n; i++ )
1655 _sy[ i ] = y[ n+i ];
1657 return __Norm_2( n+_E_g+2*n, y );
1660 template<
typename Data,
template<
class>
class Problem>
1661 double SolverConstrained<Data,Problem>::get_phi(
double *_s,
double MU,
double nu )
1664 double rhs[ _E_h + _E_g+2*_E_n ];
1667 for (
int m = 0; m < _E_h; m++ )
1670 for (
int m = 0; m < _E_g+2*n; m++ )
1671 rhs[ _E_h + m ] = gx[ m ] + _s[ m ];
1675 for (
int i = 0; i < _E_g+2*n; i++ )
1676 phi_val += -MU * std::log( _s[ i ] );
1678 phi_val += nu * __Norm_2( _E_g+2*n, rhs );
1696 template<
typename Data,
template<
class>
class Problem>
1702 M_solver_stats.clear();
1705 std::cout <<
"\n[SolverConstrained<Data,Problem>::optimize()]: Optimizing...\n";
1708 double Delta = std::min( options.Delta_init, (
double )Max_Allowed_Delta );
1709 double eMU = options.eMU_init;
1710 double MU = options.MU_init;
1711 double MU_last = MU;
1712 double nu = options.nu0_init;
1713 double rho_N = options.rho_N;
1715 double rhs_xs_p_d[ _E_h + _E_g+2*_E_n ];
1718 double x_p_dx[ _E_n ];
1719 double s[ _E_g+2*_E_n ];
1720 double ds[ _E_g+2*_E_n ];
1721 double s_p_ds[ _E_g+2*_E_n ];
1722 double lambda_h[ _E_h ], lambda_g[ _E_g+2*_E_n ];
1723 double vpred, q, pred, ared, gamma, norm_d;
1725 M_prob.copy_x0_to_x( __x );
1726 this->initialize_solver_with_x( __x, s, lambda_h, lambda_g, MU );
1728 double E = this->evaluate_E_current( s, lambda_h, lambda_g, MU );
1730 int err_iter = 0, hom_iter = 0;
1732 while ( err_iter < 500 && hom_iter < options.max_hom_iter && E > options.eTOL )
1735 double max_TR_iter = options.max_TR_iter;
1737 if ( hom_iter+1 >= options.max_hom_iter )
1738 max_TR_iter = 2 * max_TR_iter;
1743 std::cerr <<
"\n================================================================================"
1744 <<
"\n homotopy iter = " << hom_iter
1746 <<
", E(x,y,MU) = " << this->evaluate_E_current( s, lambda_h, lambda_g, MU )
1747 <<
"\n================================================================================\n";
1752 while ( iter < max_TR_iter && E > eMU )
1755 E = this->evaluate_E_current( s, lambda_h, lambda_g, MU );
1757 double Delta_used = Delta;
1758 bool truss_exit =
false, SOC =
false;
1761 bool evaluate_Hess =
false;
1764 evaluate_Hess =
true;
1766 norm_d = this->get_steps( s, Delta, MU, dx, ds, vpred, q, truss_exit, n_CGiter );
1769 nu = std::max( 1.5*nu, q / ( ( 1 - rho_N ) * vpred ) );
1771 pred = -q + nu * vpred;
1772 __SumVecs( n, __x, dx, x_p_dx );
1773 __SumVecs( _E_g+2*n, s, ds, s_p_ds );
1774 ared = get_phi( s, MU, nu ) - get_phi( x_p_dx, s_p_ds, MU, nu, rhs_xs_p_d );
1775 gamma = ared / pred;
1778 std::cerr <<
"\n-------------------------\n"
1779 <<
"\nTRiter = " << iter
1780 <<
"\n E(x,s,MU) = " << E
1781 <<
"\n gamma = " << gamma
1783 <<
"\n vpred = " << vpred
1784 <<
"\n pred = " << pred
1786 <<
"\n ared = " << ared
1787 <<
"\n norm_d = " << norm_d
1788 <<
"\n Delta = " << Delta
1791 __print_vec( n, dx );
1792 __print_vec( _E_g+2*n, ds );
1795 if ( gamma >= options.eta )
1797 __CopyVecs( n, x_p_dx, __x );
1798 __CopyVecs( _E_g+2*n, s_p_ds, s );
1800 if ( gamma >= options.rho_big )
1801 Delta = std::max( options.rho_increase_big * norm_d , Delta );
1803 else if ( gamma >= options.rho_small )
1804 Delta = std::max( options.rho_increase_small * norm_d , Delta );
1806 Delta = std::min( Delta, (
double )Max_Allowed_Delta );
1812 double xy[ _E_n ], sy[ _E_g+2*_E_n ];
1813 double norm_xypsy = this->get_SOC( s_p_ds, rhs_xs_p_d, xy, sy );
1815 std::cerr <<
"\n SECOND ORDER CORRECTION:"
1816 <<
"\n ||xy||+||sy|| = " << norm_xypsy;
1819 if ( norm_xypsy > 0 )
1821 bool sy_feasible =
true;
1822 double tau = options.tau;
1824 for (
int i = 0; i < _E_g+2*n; i++ )
1825 if ( ds[ i ] + sy[ i ] + tau*s[ i ] < 0 )
1826 sy_feasible =
false;
1830 __SumVecs( n, __x, dx, x_p_dx );
1831 __SumVecs( _E_g+2*n, s, ds, s_p_ds );
1833 __SumVecs( n, x_p_dx, xy, x_p_dx );
1834 __SumVecs( _E_g+2*n, s_p_ds, sy, s_p_ds );
1836 ared = get_phi( s, MU, nu ) - get_phi( x_p_dx, s_p_ds, MU, nu );
1837 gamma = ared / pred;
1839 std::cerr <<
"\n gamma = " << gamma;
1842 if ( gamma >= options.eta )
1845 __CopyVecs( n, x_p_dx, __x );
1846 __CopyVecs( _E_g+2*n, s_p_ds, s );
1853 Delta = options.rho_decrease * Delta;
1856 if ( gamma >= options.eta )
1858 if ( !options.Primal_Dual )
1861 this->update_opt_objs( __x, s, lambda_h, lambda_g, MU_last, evaluate_Hess );
1863 M_solver_stats.push_all( MU, E, this->evaluate_E_current( s, lambda_h, lambda_g, 0 ),
1864 n_CGiter, truss_exit, SOC,
1865 Delta_used, ared, pred, gamma );
1871 std::cerr <<
"\n E(x+Dx,s+ds,MU) = " << E <<
" ( MU = " << MU <<
")\n";
1875 MU = MU * options.theta_N;
1883 eMU = eMU * options.theta_N * eMU;
1886 Delta = options.Delta_init;
1889 nu = options.nu0_init;
1892 E = this->evaluate_E_current( s, lambda_h, lambda_g, 0 );
1895 std::cerr <<
"\n E(x+Dx,s+ds,0) = " << E <<
"\n";
1901 M_solver_stats.show();
1922 template<
typename Data,
template<
class>
class Problem>
1931 n_CGiter_hstr.clear();
1932 truss_exit_hstr.clear();
1941 template<
typename Data,
template<
class>
class Problem>
1946 x_hstr.push_back( __x[i] );
1949 template<
typename Data,
template<
class>
class Problem>
1953 bool truss_exit,
bool SOC,
1954 double Delta,
double ared,
double pred,
double gamma )
1957 mu_hstr.push_back( mu );
1958 E_hstr.push_back( E );
1959 E0_hstr.push_back( E0 );
1960 n_CGiter_hstr.push_back( n_CGiter );
1961 truss_exit_hstr.push_back( truss_exit );
1962 SOC_hstr.push_back( SOC );
1963 Delta_hstr.push_back( Delta );
1964 ared_hstr.push_back( ared );
1965 pred_hstr.push_back( pred );
1966 gamma_hstr.push_back( gamma );
1970 template<
typename Data,
template<
class>
class Problem>
1974 iter = mu_hstr.size();
1975 std::cerr <<
"\nScaled Trust-Region Method Statistics:\n\n";
1977 std::cerr <<
" k mu E(x,s;0) E(x,s;mu) CGiters trustEx SOC Delta^k ared^k pred^k gamma^k \n";
1978 std::cerr <<
"------------------------------------------------------------------------------------------------\n";
1982 for (
int k = 0; k < iter; k++ )
1984 if ( mu_hstr[ k ] != mu_last )
1986 fprintf( stderr,
"%3d ", k+1 );
1987 fprintf ( stderr,
" %3.1e", mu_hstr[ k ] );
1988 fprintf ( stderr,
" %10.3e", E0_hstr[ k ] );
1989 fprintf ( stderr,
" %10.3e", E_hstr[ k ] );
1990 fprintf ( stderr,
" %5d", n_CGiter_hstr[ k ] );
1992 if ( truss_exit_hstr[ k ] ==
true )
1993 fprintf ( stderr,
" -Y- " );
1996 fprintf ( stderr,
" N " );
1998 if ( SOC_hstr[ k ] ==
true )
1999 fprintf ( stderr,
"-Y- " );
2002 fprintf ( stderr,
" N " );
2004 fprintf ( stderr,
" %3.1e", Delta_hstr[ k ] );
2005 fprintf ( stderr,
" %9.2e", ared_hstr[ k ] );
2006 fprintf ( stderr,
" %9.2e", pred_hstr[ k ] );
2007 fprintf ( stderr,
" %7.1f", gamma_hstr[ k ] );
2008 std::cerr << std::endl;
2009 mu_last = mu_hstr[ k ];
2012 else if ( !verbose )
2014 fprintf( stderr,
"%3d ", k+1 );
2015 fprintf ( stderr,
" " );
2016 fprintf ( stderr,
" %10.3e", E0_hstr[ k ] );
2017 fprintf ( stderr,
" %10.3e", E_hstr[ k ] );
2018 fprintf ( stderr,
" %5d", n_CGiter_hstr[ k ] );
2020 if ( truss_exit_hstr[ k ] ==
true )
2021 fprintf ( stderr,
" -Y- " );
2024 fprintf ( stderr,
" N " );
2026 if ( SOC_hstr[ k ] ==
true )
2027 fprintf ( stderr,
"-Y- " );
2030 fprintf ( stderr,
" N " );
2032 fprintf ( stderr,
" %3.1e", Delta_hstr[ k ] );
2033 fprintf ( stderr,
" %9.2e", ared_hstr[ k ] );
2034 fprintf ( stderr,
" %9.2e", pred_hstr[ k ] );
2035 fprintf ( stderr,
" %7.1f", gamma_hstr[ k ] );
2036 std::cerr << std::endl;
2040 std::cerr <<
"\nNumber of accepted steps: " << iter <<
"\n";
holds the statistics for SolverConstrained
Definition: solverconstrained.hpp:153
Stats(solver_type *__s, bool __c=false)
default constructor
Definition: solverconstrained.hpp:159
void clear()
clear all statistics
Definition: solverconstrained.hpp:1924
Polynomial< Poly, Type > dx(Polynomial< Poly, Type > const &p)
compute of a polynomial p
Definition: operations.hpp:63
Problem< Data > problem_type
problem data type
Definition: solverconstrained.hpp:53
bool optimize(double __x[_E_n])
optimizes the problem defined by problem_type
Definition: solverconstrained.hpp:1698
Expr< Val< typename mpl::if_< boost::is_arithmetic< ExprT1 >, mpl::identity< Cst< ExprT1 > >, mpl::identity< ExprT1 > >::type::type > > val(ExprT1 const &__e1)
precompute expression tensor
Definition: val.hpp:304
statistics_type const & stats() const
Definition: solverconstrained.hpp:242
number of multipliers
Definition: solverconstrained.hpp:62
size of the matrix
Definition: solverconstrained.hpp:61
problem_type::ad_0_type ad_0_type
automatic differentiation type of order one
Definition: solverconstrained.hpp:67
Definition: solverconstrained.hpp:45
void push_all(double mu, double E, double E0, int n_CGiter, bool truss_exit, bool SOC, double Delta, double ared, double pred, double gamma)
insert all statistics
Definition: solverconstrained.hpp:1951
problem_type & problem()
Definition: solverconstrained.hpp:236
ad_0_type::value_type value_type
numerical type
Definition: solverconstrained.hpp:76
bool collectStats() const
tells if we collect statistics regarding the solver convergence
Definition: solverconstrained.hpp:168
void show(bool verbose=false) const
show statistics
Definition: solverconstrained.hpp:1972
void redefine_problem(double x_definitions[_E_n][3])
redefined the control variables bounds and if they are avtive
Definition: solverconstrained.hpp:224
SolverConstrained< Data, Problem > solver_type
this solver type
Definition: solverconstrained.hpp:50
~SolverConstrained()
destructor
Definition: solverconstrained.hpp:384
number of objective functions
Definition: solverconstrained.hpp:58
number of equality constraints
Definition: solverconstrained.hpp:60
problem_type::ad_1_type ad_1_type
automatic differentiation type of order one
Definition: solverconstrained.hpp:70
COptions options_type
option data type
Definition: solverconstrained.hpp:147
holds the options for the SolverConstrained
Definition: solverconstrained.hpp:87
problem_type::ad_2_type ad_2_type
automatic differentiation type of order two
Definition: solverconstrained.hpp:73
size of the multipliers matrix
Definition: solverconstrained.hpp:63
number of inequality constraints
Definition: solverconstrained.hpp:59
Stats statistics_type
statistics data type
Definition: solverconstrained.hpp:210
number of control variables
Definition: solverconstrained.hpp:57
void push_x(const double *__x)
insert x in statistics
Definition: solverconstrained.hpp:1943