% trcurve.m
% Script to draw a trust region curve
% Adjust contour position in lines marked by ***
% Note automatic plotting in final line may need file name change

% (c) Philip D. Loewen, 5 Mar 2001


% Choices for part 1
cp = [0;-1];   % Centre Point
gs = 3;        % Grid Size

% Choices for part 2
cp = [0;0.5];  % Centre Point
gs = 3;        % Grid Size: extent above, below, left, right of cp.

Dmax = 2;     % Maximum TR radius

plot(cp(1),cp(2),'.','MarkerSize',12,'color','g');
hold on;

[X,Y] = meshgrid(linspace(cp(1)-gs,cp(1)+gs,103),...
                 linspace(cp(2)-gs,cp(2)+gs,103));
f = 10*(Y-X.^2).^2 + (1-X).^2;   % given function values

x = cp(1); y = cp(2);
fcp = 10*(y-x.^2).^2 + (1-x).^2;   % given function values
gcp = [-40*x*(y-x^2)-2*(1-x), 20*(y-x^2)];
Hcp = [120*x^2-40*y+2 , -40*x;
        -40*x,           20];
M = fcp + gcp(1)*(X-x) + gcp(2)*(Y-y) + ...
      0.5*Hcp(1,1)*(X-x).^2 + Hcp(1,2)*(X-x).*(Y-y) + 0.5*Hcp(2,2)*(Y-y).^2;

[zzz,fcon] = contour(X,Y,f,f(86,1:16:103));   % *** Chg 86 if needed
set(fcon,'EdgeColor','blue');
[zzz,fcon] = contour(X,Y,f,[fcp fcp]);
set(fcon,'EdgeColor','Red','LineWidth',1);

Mmin = min(min(M)); Mmax = max(max(M)); Mvals=linspace(Mmin,Mmax,12);
[zzz,Mcon] = contour(X,Y,M,M(86,1:20:52));    % *** Chg first 86 if needed
set(Mcon,'EdgeColor','cyan');
[zzz,Mcon] = contour(X,Y,M,[fcp fcp]);
set(Mcon,'EdgeColor','Red','LineWidth',1);

pN = - Hcp\gcp';  % Newton Step
XN = cp(1) + pN(1); YN = cp(2) + pN(2);
MN = fcp + gcp(1)*(XN-x) + gcp(2)*(YN-y) + ...
     0.5*Hcp(1,1)*(XN-x).^2 + Hcp(1,2)*(XN-x).*(YN-y) + 0.5*Hcp(2,2)*(YN-y).^2;

nD = 50;
dD = Dmax/nD;
xtc = zeros(1,nD);
ytc = zeros(1,nD);

th = linspace(0,2*pi,361);
ct = cos(th);
st = sin(th);

for k=1:nD,
  Delta = k*dD;
  X = cp(1) + Delta*ct;
  Y = cp(2) + Delta*st;
  M = fcp + gcp(1)*(X-x) + gcp(2)*(Y-y) + ...
      0.5*Hcp(1,1)*(X-x).^2 + Hcp(1,2)*(X-x).*(Y-y) + 0.5*Hcp(2,2)*(Y-y).^2;
  lowpt = find(M==min(M));
  lowpt = lowpt(1);
  % disp(['M(lowpt)=',num2str(M(lowpt))]);
  if norm([XN;YN]-cp,2)>Delta | MN > M(lowpt),
    xtc(k) = X(lowpt);
    ytc(k) = Y(lowpt);
  else
    xtc(k) = XN;
    ytc(k) = YN;
  end;
end;

plot([cp(1),xtc],[cp(2),ytc],'g','LineWidth',1);

set(gca,'XLim',[cp(1)-gs,cp(1)+gs],'YLim',[cp(2)-gs,cp(2)+gs]);
set(gcf,'PaperPosition',[1,1,4,4]);
%print -deps trcurve1.eps

