% Define a quadratic functional by Q(x) = || Ax-b ||^2. It arises as a
% discrepancy term for the linear equation Ax=b. Here x and b are
% 2-dimensional vectors.
%
% Let w be a nonzero 2-dimensional vector. Further, denote the location of
% the minimum of Q by z. In other words,
%
% Q(z) = min{Q(x)},
%
% where the minimum is taken over all points x in the plane. We assume
% here that the minimum z is unique.
%
% We wish to show visually why the t-derivative of Q(z+tw) vanishes at t=0.
%
% Samuli Siltanen October 2015

% Parameters for controlling the appearance of the plots
thinline  = 1;
thickline = 2;
fsize     = 16;

% Choose direction vector w
w = [-.2;1];
w = w/norm(w);

% Define A and b
A = [[1 0.2];[0 1]];
b = [0;0];

% construct evaluation points for t
Nt = 12;
t  = linspace(-1,1,Nt);

% Construct the parameterized path in the x1-x2 plane
path1 = w(1)*t;
path2 = w(2)*t;

% Create 2D evaluation points
[X1,X2] = meshgrid(linspace(-1,1,128));

% Evaluate Q
Q = Qfunctional(X1,X2,A,b);
Q(abs(X1+1i*X2)>1) = NaN;

% Plot the result
for iii = 2:Nt
    figure(1)
    clf
    subplot(1,2,1)
    mesh(X1,X2,Q)
    hold on
    plot3([path1(1),path1(iii)],[path2(1),path2(iii)],[0,0],'b','linewidth',thinline)
    plot3(...
        [path1(iii),path1(iii)],...
        [path2(iii),path2(iii)],...
        [0,Qfunctional(path1(iii),path2(iii),A,b)],'k','linewidth',thinline)
    axis square
    axis([-1 1 -1 1 0 1])
    axis square
    title(['w = [', num2str(w(1),'%1.2f'), ',', num2str(w(2),'%1.2f'),']'],'fontsize',fsize)
    set(gca,'xtick',[-1,0,1],'fontsize',fsize)
    set(gca,'ytick',[-1,0,1],'fontsize',fsize)
    set(gca,'ztick',[0,.5,1],'fontsize',fsize)
    
    % Record camera parameters
    cva = get(gca,'CameraViewAngle');
    cp  = get(gca,'CameraPosition');
    ct  = get(gca,'CameraTarget');
    
    subplot(1,2,2)
    surf(X1,X2,Q)
    shading interp
    alpha 0.2
    hold on
    plot3([path1(1),path1(iii)],[path2(1),path2(iii)],[0,0],'b','linewidth',thinline)
    hold on
    plot3(...
        [path1(iii),path1(iii)],...
        [path2(iii),path2(iii)],...
        [0,Qfunctional(path1(iii),path2(iii),A,b)],'b','linewidth',thinline)
    plot3(path1(1:iii),path2(1:iii),Qfunctional(path1(1:iii),path2(1:iii),A,b),...
        'r','linewidth',thickline)
    axis square
    axis([-1 1 -1 1 0 1])
    axis square
    set(gca,'CameraViewAngle',cva,'CameraPosition',cp,'CameraTarget',ct);
    title(['t = ', num2str(t(iii),'%1.2f')],'fontsize',fsize)
    set(gca,'xtick',[-1,0,1],'fontsize',fsize)
    set(gca,'ytick',[-1,0,1],'fontsize',fsize)
    set(gca,'ztick',[0,.5,1],'fontsize',fsize)

    
    pause
    
    
    
    % Save image to file
    %eval(['print -dpng frames/Qpic', num2str(iii,'%03d'), '.png'])
end

% % Crop the images
% for iii = 2:Nt
%     % Read in the image
%     filename = ['frames/Qpic', num2str(iii,'%03d'), '.png'];
%     im = imread(filename,'png');
%     
%     % Crop the image
%     imR = im(186:708,100:1120,1);
%     imG = im(186:708,100:1120,2);
%     imB = im(186:708,100:1120,3);
%     im2 = zeros(size(imR,1),size(imR,2),3);
%     im2(:,:,1) = imR;
%     im2(:,:,2) = imG;
%     im2(:,:,3) = imB;
%     
%     % Write the image to file
%     imwrite(im2,filename,'png');
%     
%     % Monitor the run
%     if mod(iii,10)==0
%         disp(['Cropped ',num2str(iii),' out of ',num2str(Nt)])
%     end
%     
% end