%Solves a small maze, similar to the one designed in MP2
%by various methods.
clear
% the maze is 4 x 4, with food along the left side only.
ndim = 4;
P = zeros(ndim,ndim);
%initialize all P's and then set the constant border
% use uniform initial guess set below.
P(:,:) = 0.50;
%corners are inaccessible
P(1,1) = NaN;
P(4,4) = NaN;
P(1,4) = NaN;
P(4,1) = NaN;
%left side has food
P(2,1) = 1;
P(3,1) = 1;
%other sides do not have food
P(1,2) = 0;
P(1,3) = 0;
P(4,2) = 0;
P (4,3) = 0;
P(2,4)=0;
P(3,4) = 0;
PInit = P; %store copy of initial P matrix
%
% run iteratove loop over the interior variable cells
%this is a Gauss Seidel iteration, done explicitly over the
%probabilities to be determined.
tol = 1.E-6;
%iterate to convergence
converged = false;
nIter = 0;
while (~converged)
    Pold = P;
    for i = 2:3
        for j = 2:3
            %each cell probability is average of the neighbors
            P(i,j) = (0.25)*(P(i,j+1) + P(i,j-1) + P(i+1,j) + P(i-1,j));
        end
    end
    %check tolerance
    nIter = nIter + 1;
    error = max(max (abs(P-Pold)));
    if (error < tol)
        converged = true;
        disp ('Gauss Seidel Method')
        nIter
        P     %print out the maze of probabilities
    end
end
% run iteratove loop over the interior variable cells
%this is a Jacoby iteration, done explicitly over the
%probabilities to be determined. Uses only values from
%previous iteration in formulas
%iterate to convergence
%reset P to initial values
P = PInit;
converged = false;
nIter = 0;
while (~converged)
    Pold = P;
    for i = 2:3
        for j = 2:3
            %each cell probability is average of the neighbors
            P(i,j) = (0.25)*(Pold(i,j+1) + Pold(i,j-1) + ...
                                  Pold(i+1,j) + Pold(i-1,j));
        end
    end
    %check tolerance
    nIter = nIter + 1;
    error = max(max (abs(P-Pold)));
    if (error < tol)
        converged = true;
        disp ('Jacobi Explicit Method')
        nIter
        P     %print out the maze of probabilities
    end
end


% solve by matrix methods
% contains coefficients of variable terms in the equations above
A = [1 -.25 -.25 0; -.25 1 0 -.25; -.25 0 1 -.25; 0 -.25 -.25 1];
b = [0.25; 0; 0.25; 0];   %contains constant terms from equations
%Solve for x
%X is the set of probabilities for the cells
% Built in matrix method
disp ('Built-in Matrix Method')
X = A\b
%
%solve by iterative Jacobi method (p. 236 of turner)
%A and b from above 
D = diag(A);
n = size(A,1);
A_D = A - diag(D);
x = zeros(n,1);
nits = 15;
for k = 1:nits
    x = (b - A_D*x)./D;
    S(:,k)=x;
end
%each column of S is an iterative set of the variables.
%Iterative Jacoby method in matrix form
S