function pi=maze(N,rho,alpha,M,start,finish)
% pi=maze(N,rho,alpha,M,start,finish)
% Outputs -
% pi: Output choice set
% Inputs -
% N: Number of samples to generate each round
% rho: fraction of best samples to take
% alpha: smoothing parameter
% M: maze (0 is free, 1 is blocked)
% start: starting co-ords
% finish: endinf co-ords

% works very poorly for large mazes!
% since we (on most its) "DIE" before we glean too
% much useful information from the particular chosen
% direction set.

z=5;

     if nargin<=3,

%			   start=[1,1];
%finish=[20,20];
%direc=2;
%M=[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%    1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 0];


			   start=[1,2];
finish=[20,19];
direc=2;
M=[ 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
    1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 0 1];



%			   start=[1,3];
%finish = [12,8];
% direc=3;
%M=[1 1 0 1 1 1 1 1 1 1 1 1;
%   1 0 0 0 0 0 1 1 0 0 0 1;
%   1 0 1 1 1 0 0 0 0 1 0 1;
%   1 0 0 1 0 0 1 1 1 1 0 1;
%   1 1 0 1 0 1 1 1 1 0 0 1;
%   1 1 0 1 0 1 1 1 1 0 1 1;
%   1 0 0 0 0 0 0 0 0 0 1 1;
%   1 1 1 1 0 1 1 1 0 1 1 1;
%   1 0 0 0 0 0 1 0 0 0 0 1;
%   1 0 1 1 1 0 1 0 1 1 0 1;
%   1 0 0 1 0 0 0 0 1 1 0 1;
%   1 1 1 1 1 1 1 0 1 1 1 1];
%start=[1,2];
%finish=[4,3];
% direc=3;
%M=[ 1 0 1 1;
%    1 0 0 1;
%    1 1 0 1;
%    1 1 0 1];

%start=[1,2];
%finish=[7,6];
% direc=3
%M=[1 0 1 1 1 1 1;
%   1 0 0 0 1 1 1;
%   1 0 1 0 1 1 1;
%   1 0 1 0 0 0 1;
%   1 0 1 1 1 0 1;
%   1 0 1 1 1 0 1;
%   1 1 1 1 1 0 1];

end

[m,n]=size(M);

dn=0;

for i=1:(m-2),
			   for j=1:(n-2),

			   if (M(i+1,j)==0)&&(M(i+1,j+1)==0)&&(M(i,j+1)==0)&&(M(i+2,j+1)==0)&&(M(i+1,j+2)==0),
			   dn=dn+1;
elseif(M(i+1,j)==0)&&(M(i+1,j+1)==0)&&(M(i,j+1)==0)&&(M(i+2,j+1)==0),
			   dn=dn+1;
elseif(M(i+1,j+1)==0)&&(M(i,j+1)==0)&&(M(i+2,j+1)==0)&&(M(i+1,j+2)==0),
			   dn=dn+1;
elseif(M(i+1,j)==0)&&(M(i+1,j+1)==0)&&(M(i+2,j+1)==0)&&(M(i+1,j+2)==0),
			   dn=dn+1;
elseif(M(i+1,j)==0)&&(M(i+1,j+1)==0)&&(M(i,j+1)==0)&&(M(i+1,j+2)==0),
			   dn=dn+1;
end
end
end

P=ones(dn,4);
P=P./4
Pold=zeros(dn,4);

if nargin==0,
N=500;
	       rho=.01;
alpha=.3;
end

tol=0.005;

kk=(1:n)';

% P(i,k) : prob that at desc. node i we move in direction "k"
%
%            1 
%          4 k 2
%            3
%

st=zeros(z,1);

count=1;

g=floor(rho*N);

while max(max(abs(P-Pold)))>tol,

for i=1:N,

X(:,i)=gmaze(P)';
S(i)=smaze(X(:,i),M,start,finish,direc);

end

[S,I]=sort(S);

w=zeros(dn,4);

for i=1:4,
w(:,i)=sum(X(:,I(1:g))==i,2);
end

w=w./g;

Pold=P;
P=alpha*w+(1-alpha)*P;
st(1:(z-1))=st(2:z);
st(z)=S(g);
[count,S(g),S(1)]

if (sum(st==st(z))==z),
break;
end
count=count+1;
end
pi=X(:,I(1));

