%
% Script to describe some Trust-region ideas by
% examining the structure of the quadratic sub-problem.
% We generate level curves for a quadratic functional,
% plot the level curves and then add
% some features of interest in Trust - Region formulations.
% Specifically, we show the Newton and Cauchy steps and
% where appropriate, the Powell dog-leg.
%
 global G   c
%
  G = [ 2 1 ; 1 1 ]; c = [-2; -2];
%	  
 figure % open a new figure
 grid
 axis('square')
 hold
%
% Find the Newton and the Cauchy steps
% These are used in scaling the contours
% for the quadratic
%
 Newt  = -G\c;
 Cs    =  -((c'*c)/(c'*G*c))*c; 
%
% the extent of these arrays should be tuned for the problem
%
  half = 4; inc = 0.1;  
  xv = [Newt(1) - half:inc:Newt(1)+half]; xv = xv(:); 
  yv = [Newt(2) - half:inc:Newt(2)+half]; yv = yv(:); 
  nxv = length(xv); nyv = length(yv);
%
  zv = zeros(nxv,1); Z = zeros(nxv,nyv);
%
% generate data for the contour plot of the quadratic functional
%
  for j=1:nyv
     vec(2) = yv(j);
	 for i=1:nxv
	     vec(1) = xv(i);
	     zv(i) = quadratic(G,c,vec);
     end
	 Z(:,j) = zv(:);
  end
%
% display the contour lines
%
% Find the min-value for the quadratic
%
 v_min = quadratic(G,c,Newt); v_max = max(max(Z));
 v_cs  = quadratic(G,c,Cs);
%
% We step from the Newton-point to the Cauchy-point and then along the
% extended line until the quadratic exceeds v_max
%
  vec = Newt; step = Cs - Newt; vals = []; val = v_min;
  while val < v_max
     vec = vec + step ; val = quadratic(G,c,vec);
	 vals = [vals val];
  end
 contour(xv,yv,Z',vals)
%
% Add the Cauchy and the Newton steps
%
  plot([0 ; Cs(1)],[0 ;Cs(2)],'g')
 plot(Cs(1),Cs(2),'ro')
 text(0.7*Cs(1)+.1,0.7*Cs(2),'Cauchy')  
%
 plot([0; Newt(1)],[0; Newt(2)],'b')
 plot(Newt(1),Newt(2),'ro')
 text(.7*Newt(1) - .7,.7*Newt(2),'Newton')
%
% Add the dog-leg from Cauchy to Newton
%
  tmp = [Cs(:)'; Newt(:)']';
  plot(tmp(1,:),tmp(2,:),'g--')
% Show a trust-region 
% a circle of specified radius 
%
 rad = 1.5;
 thet = [0: pi/50: 2*pi]; xt = rad*cos(thet); yt = rad*sin(thet);	
 plot(xt,yt,'w')
%
% generate the locus of points p-hat(lambda)
%
 np = 20;  Ph=[0;0]; lambda_0 = 0.1; scl = 4;
 for i=1:np
    lambda = lambda_0 * exp(i/scl);
    Ph = -(G + lambda*eye(2))\c;
	plot(Ph(1),Ph(2),'+r')
%	if (i ==np/4)
%	   text(Ph(1),Ph(2),'P-hat locus')
%	end
 end
% 
  
