% A comparision of gradient descent and conjugate gradient (guess who wins)

function main()

   % data
   A=[17, 2; 2, 7];  % the matrix 
   b=[2, 2]';        % right-hand side
   x0=[0, 0]';       % the initial guess  

   % linewidth and font size
   lw= 2; 
   fs = 25;

% colors
   red=[0.867 0.06 0.14];
   blue = [0, 129, 205]/256;
   green = [0, 200,  70]/256;
   black = [0, 0, 0];
   white = 0.99*[1, 1, 1];

   % Set up the plotting window
   figure(1); clf; set(gca, 'fontsize', fs); hold on; axis equal; axis off; 
   s = 0.16; x = A\b;
   Ax = x(1)-s; Bx = x(1)+s; Ay = x(2)-2.0*s; By =  x(2)+s;
   plot([Ax Bx Bx Ax Ax], [Ay Ay By By Ay], 'color', blue, 'linewidth', lw/2); % plot a blue box
   s=0.005; plot(Ax-s, Ay-s, '*', 'color', white); plot(Bx+0.5*s, By+0.5*s, '*', 'color', white); %markers
   Box = [Ax Bx Ay By];
   axis (Box);
   
   % plot the contours of the quadratic form associated with A and b
   plot_contours(A, b, Box, lw, blue);

% Do conjugate gradient and gradient descent.
% For the first one, start a bit shifted so that the two graphs don't overlap.
   shift = 0.0015*[1, -1];
   small_rad=0.002;
   tol = eps;
   x = conj_gradient(A, b, x0, tol, lw, red,   small_rad,  shift);
   x = grad_descent (A, b, x0, tol, lw, green, small_rad);

% text
   small = 0.015;
   text(x0(1)-2*small, x0(2)-1.6*small, 'x', 'fontsize', fs); 
   text(x0(1)-0.5*small, x0(2)-3*small, '0', 'fontsize', floor(0.7*fs)); 
   text(x(1)+small, x(2)+small, 'x', 'fontsize', fs);

% some balls for beauty
   small_rad = 0.003;
   ball(x0(1)+shift(1)/2, x0(2)+shift(2)/2, small_rad, blue)
   ball(x(1), x(2), small_rad, blue)

% save to disk as eps and svg
   saveas(gcf, 'Conjugate_gradient_illustration.eps', 'psc2');
   plot2svg('Conjugate_gradient_illustration.svg');
   
   
function x = conj_gradient(A, b, x, tol, lw, color, small_rad, shift)

   r=A*x - b;
   d=-r;

   while norm(r) > tol

      % a pretty ball for beauty, to cover imperfections when two segments are joined
      ball(x(1)+shift(1), x(2)+shift(2), small_rad, color);

      alpha = -dot(r, d)/dot(A*d, d);
      x0 = x;
      x = x + alpha*d;
      
      r=A*x - b;
      beta = dot(A*r, d)/dot(A*d, d);
      d0 = d;
      d = -r + beta*d;

      plot([x0(1), x(1)]+shift(1), [x0(2), x(2)]+shift(2), 'color', color, 'linewidth', lw)

   end

function x = grad_descent(A, b, x, tol, lw, color, small_rad)

   r=A*x - b;
   d=-r;

   while norm(r) > tol

      % a pretty ball for beauty, to cover imperfections when two segments are joined
      ball(x(1), x(2), small_rad, color);
      
      alpha = -dot(r, d)/dot(A*d, d);
      x0 = x;
      x = x + alpha*d;
      
      r=A*x - b;
      beta = 0; %beta = dot(A*r, d)/dot(A*d, d);
      d0 = d;
      d = -r + beta*d;

      plot([x0(1), x(1)], [x0(2), x(2)], 'color', color, 'linewidth', lw)

   end

   

function plot_contours (A, b, Box, lw, color);

   N=200;  % number of points (don't make it big, code will be slow)

   E = A\b; % the exact solution, around which we will draw the contours
   B = 0.12;
   [X, Y]=meshgrid(linspace(Box(1)-B, Box(2)+B, N), linspace(Box(3)-B, Box(4)+B, N));     % X and Y coordinates

% the quadratic form f= (1/2)*x'*A*X-b'*x;
   f = inline('0.5*A(1, 1)*X.*X + A(1, 2)*X.*Y+0.5*A(2, 2)*Y.*Y-b(1)*X-b(2)*Y', 'X', 'Y', 'A', 'b');
   Z = 0.5*A(1, 1)*X.*X + A(1, 2)*X.*Y+0.5*A(2, 2)*Y.*Y-b(1)*X-b(2)*Y;

   
   % prepare to draw the contours
   x0 = A\b; f0 = f(x0(1), x0(2), A, b);
   No = 25; % number of contours
   Levels = (linspace(f0, 1, No)-f0).^2+f0;

% Plot the contours with 'contour' in figure(2), and then with 'plot' in figure(1).
% This is to avoid a bug in plot2svg, it can't save output of 'contour'.
   figure(2); clf; hold on;
   for i=1:length(Levels)

      figure(2);
      [c, stuff] = contour(X, Y, Z, [Levels(i), Levels(i)]);

      [m, n]=size(c);
      if m > 1 & n > 0

	 % extract the contour from the contour matrix and plot in figure(1)
	 l=c(2, 1);
	 x=c(1,2:(l+1));  y=c(2,2:(l+1)); 
	 figure(1); plot(x, y, 'color', color, 'linewidth', lw/2);

      end
   end
   figure(1);

function ball(x, y, r, color)
   Theta=0:0.1:2*pi;
   X=r*cos(Theta)+x;
   Y=r*sin(Theta)+y;
   H=fill(X, Y, color);
   set(H, 'EdgeColor', 'none');