%
%  trust-region solution for non-convex quadratic (easy case)
%
%  Nick Gould

lambda_min = -0.6;
lambda_max = 1;
psi_min = -0.3;
psi_max = 0.7;

tick = 0.03 ;
lambda_mint = -0.63 ;
psi_mint = -0.32;

inc = .005;

H = [-0.5 0.5];
g = [0.5 1];
Delta = 4 ;
lambda_crit = 0.5 ;
lambda_fade = lambda_min:inc:lambda_crit;
lambda = lambda_crit:inc:lambda_max;

hold off

s_fade = g(1) ./ ( H(1) + lambda_fade ) ;
psi_fade = s_fade.^2 ;
s_fade = g(2) ./ ( H(2) + lambda_fade ) ;
psi_fade = psi_fade + s_fade.^2 ;

%plot( lambda_fade, psi_fade, '--b' );
phi_fade = 1 ./ sqrt( psi_fade ) ;
phi_fade = phi_fade - ( 1 ./ Delta ) ;
plot( lambda_fade, phi_fade, '--b' );

axis( [ lambda_mint lambda_max psi_mint psi_max ] );

hold on

s = g(1) ./ ( H(1) + lambda ) ;
psi = s.^2 ;
s = g(2) ./ ( H(2) + lambda ) ;
psi = psi + s.^2 ;

%plot( lambda, psi, '-r', 'LineWidth', 2 );

phi = 1 ./ sqrt( psi ) ;
phi = phi - ( 1 / Delta ) ;
plot( lambda, phi, '-r', 'LineWidth', 2 );


crit = [ lambda_crit lambda_crit ] ;
vertical = [ psi_min psi_max ] ;
%plot( crit, vertical, ':k' );

axis off

lambda_star = 0.62818687 ;

psi_star = ( 1 / sqrt(  ( g(1) / ( H(1) + lambda_star ) )^2 + ...
                        ( g(2) / ( H(2) + lambda_star ) )^2 ) ) ...
             - ( 1 / Delta ) ;

plot( [ lambda_min lambda_max ], [ 0 0 ], ':k' );
plot( [ lambda_min lambda_max ], [ psi_min psi_min ], '-k' );
plot( [ lambda_min lambda_min ], [ psi_min psi_max ], '-k' );
plot( [ 0 0 ], [ psi_mint psi_min ], '-k' );
plot( [ lambda_crit lambda_crit ], [ psi_mint psi_min ], '-k' );
plot( [ lambda_star lambda_star ], [ psi_mint psi_min ], '-k' );
plot( [ lambda_mint lambda_min ], [ 0 0 ], '-k' ) ;

hold off

print -depsc eps/tr4

