%-----------------------------------------------------
%
% Matlab script corresponding to a hill height problem.
% The data involved had:  
%  
% From ground level,
%   A is measured 1236 meters high
%   B is measured 1941 meters high
%   C is measured 2417 meters high
%  
% From A, B is 711 meters higher 
% From A, C is 1177 meters higher 
% From B, C is 475 meters higher 
%  
% Compared to test_hills1.m, this has more output and
% plots the residual vector and its absolute value, as
% a measure of the accuracy of the observations.
% test_hills1.m and experiment_least_squares1.m had the
% least squares function doing the plotting and output,
% but this version shifts it to the driver test_hills2.m
% so that the least squares function can be used for more
% general problems.
%
%-------------------------
%
% Started : 21 Jan 2009
% Last Modified: Tue 24 Feb 2009, 05:46 PM
%
%-----------------------------------------------------

A = [ 1 0 0
      0 1 0
      0 0 1
      -1 1 0
      -1 0 1
      0 -1 1];

b = [1236; 1941; 2417; 711; 1177; 475];

[h, outlier, outlier_index] = experiment_least_squares2(A, b);

[number_of_observations, number_of_parameters] = size(A);
h = A\b;

s = ['Best least squares fit for data is A = ', num2str(h(1)), ...
     ', B = ', num2str(h(2)), ', C = ', num2str(h(3))];
disp(s);
    
s = ['Worst outlier is observation %s', num2str(outlier_index)];
disp(s);

residual = b - A*h;
display(['Residual vector is ', num2str(residual')]);

figure;
plot(1:number_of_observations, residual, 'r-',...
     1:number_of_observations, residual, 'b*');
title('Residual Vector');
xlabel('Residual Vector Component');
ylabel('Residual Vector Values');
figure;

figure;
plot(1:number_of_observations, abs(residual), 'r-',...
     1:number_of_observations, abs(residual), 'b*');
title('Absolute Values of Residual Vector');
xlabel('Residual Vector Component');
ylabel('Residual Vector Values');

[outlier, outlier_index] = max(abs(residual));

s = sprintf('Worst outlier is observation number %s', num2str(outlier_index));
disp(s);

return