Algorithm Alley
by Shehrzad Qureshi


Listing One
function fixed_pt_iteration_1
% FIXED_PT_ITERATION_1   un-optimized fixed-point iteration
num_it = 100; % # of iterations
num_vars = 2;
sigma_term = .125;

load variables.mat; % get F, x & y axes

% get starting point for fixed-pt iteration
fixed_pt_solution = get_initial_guess(x_axis, y_axis, F);

length_x_axis = length(x_axis);
length_y_axis = length(y_axis);

for n=1:num_it
  denominator = 0; % denominator of expression given by equations 2(a) & 2(b) 
  numerator = [0 0]; % numerator of expression given by equations 2(a) & 2(b)
    
  for ii = 1:length(x_axis)
        for jj = 1:length(y_axis)
            xy = fixed_pt_solution; % [ xi yi ]
            axis_pts = [x_axis(ii) y_axis(jj)];
            term = F(ii,jj)*gaussian(xy, axis_pts, sigma_term);
            denominator = denominator + term;
            numerator(1) = numerator(1) + x_axis(ii) * term;
            numerator(2) = numerator(2) + y_axis(jj) * term;
        end
    end
    % update refined answer
    fixed_pt_solution = numerator./denominator;    
end
disp(sprintf('Answer is (%f,%f)',fixed_pt_solution(1),fixed_pt_solution(2)));
function initial_guess = get_initial_guess(x_axis, y_axis, F)
%
% GET_INITIAL_GUESS - returns starting point for fixed point iteration
%                     evaluates expression given by 3(a) & 3(b)
x_initial = 0;
y_initial = 0;
[M N] = size(F);

for ii=1:M
    for jj=1:N
        x_initial = x_initial + x_axis(ii) * F(ii,jj);
        y_initial = y_initial + y_axis(jj) * F(ii,jj);
    end
end
initial_guess = [ x_initial y_initial ] ./ sum(F(:));
function term = gaussian(xy, axis_pts, sigma_term)
%
% GAUSSIAN - evaluates 2-dimensional gaussian given by equation in Fig. 4

term = 1.0;
for dim=1:2
    term = term * exp(-(xy(dim)-axis_pts(dim)) * 
                                   (xy(dim)-axis_pts(dim)) / (sigma_term));
end


Listing Two
function fixed_pt_iteration_2
% FIXED_PT_ITERATION_2   optimized fixed-point iteration

num_it = 100; % # of iterations
num_vars = 2;
sigma_term = .125;

load variables.mat; % get F, x & y axes
[M N] = size(F);

% form x_grid by tiling x' along column dimension
x_grid = repmat(x_axis', [1 N]);
% same for y dimension by except tile along the row dimension
y_grid = repmat(y_axis, [9 1]);

% get starting point for fixed-pt iteration
fixed_pt_solution = get_initial_guess(x_axis, y_axis, F);

for n=1:num_it
    x_approx = repmat(fixed_pt_solution(1),[9 9]);
    y_approx = repmat(fixed_pt_solution(2),[9 9]);
    
    x_exp = exp( -((x_approx-x_grid).*(x_approx-x_grid))./sigma_term );
    y_exp = exp( -((y_approx-y_grid).*(y_approx-y_grid))./sigma_term );
    
    term1 = F .* x_exp .* y_exp;
    denominator = sum(term1(:));
    
    term2_x = x_grid .* term1;
    term2_y = y_grid .* term1;
    
    numerator = [ sum(term2_x(:)) sum(term2_y(:)) ];
    
    fixed_pt_solution = numerator./denominator;
end

disp(sprintf('Answer is (%f,%f)',fixed_pt_solution(1),fixed_pt_solution(2)));
function initial_guess = get_initial_guess(x_axis, y_axis, F)
%
% GET_INITIAL_GUESS - returns starting point for fixed point iteration
%                     evaluates expression given by 3(a) & 3(b)
%
x_initial = 0;
y_initial = 0;
[M N] = size(F);

for ii=1:M
    for jj=1:N
        x_initial = x_initial + x_axis(ii) * F(ii,jj);
        y_initial = y_initial + y_axis(jj) * F(ii,jj);
    end
end
initial_guess = [ x_initial y_initial ] ./ sum(F(:));

Listing Three
function profile_fixed_pt(which_implementation, N)
% PROFILE_FIXED_PT   demonstrates use of Matlab profiler
%    LOAD_DRR(which_implementation, N) runs the function indicated
%    by which_implementation N times.
%    Use '1' for fixed_pt_iteration_1 (slow implementation)
%        '2' for fixed_pt_iteration_2 (optimized implementation)

if which_implementation == 1
    fun_str = 'fixed_pt_iteration_1';
elseif which_implementation == 2
    fun_str = 'fixed_pt_iteration_2';
else
    error('Must specify 1 or 2');
end
profile on;
for ii=1:N
    eval(fun_str);
end
profile off;
profile report;






3