%-----------------------------------------------------------------------
% 
% Driver script to compare how long it takes to create a 2D array of
% random values, using functions instead of scripts and using explicit
% loops versus Matlab's array notation.  This results in four
% different methods:
% 
% 1.    Using a script with internal loops: use a double-nested loop
% 2.    Using a script without internal loops:
% 3.    Using a function with internal loops:
% 4.    Using a function without internal loops:
% 
% The first and third ones just do something like
% 
%      for i = 1:n
%          for j = 1:n
%             H(j,i) = rand(1,1);
%          end
%      end
% 
% while the second and fourth do
% 
%     H  = rand(n,n);
% 
% According to hoary old wisdom, using loops is much slower than using
% array notation. Also, because a Matlab function is compiled the first
% time it is executed while scripts are interpreted line-by-line, using
% a script is much slower than writing things as a function. This script
% will test those hypotheses.
%
%       Created: 22 Jan 2008
% Last Modified: Wed 06 Feb 2008, 09:04 AM
%
%-----------------------------------------------------------------------

%-----------------------------
% How many timings to perform.
%-----------------------------
ntests = 100;

%---------------------------------------------------------
% Size of the matrix to be constructed. The time required
% will grow quadratically with this, so doubling this size
% will lead to a 4x increase in the amount of time.
%---------------------------------------------------------
n = 200;

%------------------------------------------------------------
% Names for indices, to make it easier to remember. That is,
% it is easier to figure out what the array times(fnoloops,:)
% references than times(4,:).  
%------------------------------------------------------------
sloops   = 1;
snoloops = 2;
floops   = 3;
fnoloops = 4;
%--------------------------------------------------------
% Array of character strings used for the legend command
%--------------------------------------------------------
methodnames = strvcat('sloops', 'snoloops', 'floops', 'fnoloops');

times = zeros(4,ntests);

for k = 1:ntests
    % Using a script with internal loops:
    tic;
    scriptloops;
    times(sloops, k) = toc;
    clear H i j;

    % Using a script without internal loops:
    tic;
    scriptnoloops;
    times(snoloops, k) = toc;
    clear H i j;

    % Using a function with internal loops:
    tic;
    H = functionloops(n);
    times(floops, k) = toc;
    clear H i j;

    % Using a function without internal loops:
    tic;
    H = functionnoloops(n);
    times(fnoloops, k) = toc;
    clear H i j;
end

%----------------------------
% Statistics for the timings
%----------------------------
avg_times = mean(times');  % Do "help mean" to see why I had to use the
                           % tranpose (times')  of the timing vector,
                           % instead of just doing mean(times).

if (ntests > 1)     % If a single timing, no point in doing std deviation
    std_times = std(times');
    errorbar([1:4], avg_times, -std_times, std_times, 'b*' )
    title('Average times with \pm error bars')
else
    plot([1:4], avg_times, 'b*' )
end

% Try the following with the more obvious abcissa [1:4] to see why
% I added 0.1 to each. This attaches given text to given coordinates.
text([1.1:1:4.1], avg_times, methodnames);


%---------------------------------------------------------------
% Look at all of the timing data. Gives a better feel for the 
% reliability of the data compared to just looking at averages 
% and the standard deviation, and is computationally reasonable
% given there's only 4*ntest ~ 400 observations to plot.
%---------------------------------------------------------------
figure; % Open a new graphics window to avoid losing the first one.
x = 1:ntests; % An array used for the x-coordinates in plotting. Notice
              % that you can change this to 2:ntests to avoid the first
              % measurement if it is a special case. Or use 2:2:ntests
              % to get every even-numbered timing. But in neither case
              % is it necessary to change the plot command below.

plot(x, times(sloops,x),   'b*', ...
     x, times(snoloops,x), 'r*', ...
     x, times(floops,x),   'bo', ... 
     x, times(fnoloops,x), 'ro');
legend(methodnames);
xlabel('Timing number');
ylabel('Time (sec)');
title('Data from testing scripts vs. functions and loops vs. vector ops')