%
%  math495/stat490 -- 12 nov 03 -- djm
%
%  w11snakes.m:

clear

%  set up markov matrix
P = zeros(101,101);

%  simple rolls
for k=1:6
    P = P + diag(ones(101-k,1),k)/6;
end

%  end conditions
P(101,101) = 1;

for k=1:5
    P(101-k,101) = (7-k)/6;
%  exact roll rule
%    P(101-k,101-k) = (6-k)/6;
end

%  10 snakes & 9 ladders
sinfo = 1 + [14  4; 31  9; 44 26; 56 53  ; 62 19; 64 42; ...
             84 28; 91 71; 95 75; 98 80 ];
linfo = 1 + [ 1 38;  6 16; 11 49; 21 60  ; 24 87; ...
             35 54; 51 67; 73 92; 78 100];
    
%  column changes
for k=1:size(sinfo,1)
    P(:,sinfo(k,2)) = P(:,sinfo(k,2)) + P(:,sinfo(k,1));
    P(:,sinfo(k,1)) = zeros(101,1);
end

for k=1:size(linfo,1)
    P(:,linfo(k,2)) = P(:,linfo(k,2)) + P(:,linfo(k,1));
    P(:,linfo(k,1)) = zeros(101,1);
end

%  row changes
for k=1:size(sinfo,1)
    P(sinfo(k,1),:) = P(sinfo(k,2),:);
end

for k=1:size(linfo,1)
    P(linfo(k,1),:) = P(linfo(k,2),:);
end

%  make square #78 an absorbing state
P(78+1,101) = 0;  P(78+1,78+1) = 1;

%  transpose & linear algebra
PT = transpose(P);

[evec,eval] = eig(PT);

figure(1);  clf;
image(1:101,1:101,P,'cdatamapping','scaled');  colorbar
title('\bf markov matrix')

figure(2);  clf
subplot(2,1,1)
plot(1:101,abs(diag(eval)),'x')
title('\bf absolute value of eigenvalues')
xlabel('\bf index')
ylabel('\bf |\lambda_k|')
axis([1 101 0 1])

subplot(2,1,2)
plot(0:100,evec(:,3),'x');  hold on
title('\bf eigenvector of slowest decay')
xlabel('\bf index (red circles are at snake tails!?)')
ylabel('\bf vector elements')
plot(sinfo(:,2)-1,evec(sinfo(:,2),3),'ro')
axis([0 100 -0.04 0.005])

%  mean # of steps to end
b = ones(99,1);
S = inv(eye(99,99)-P([1:78 80:100],[1:78 80:100]));
msteps = (eye(99,99)-P([1:78 80:100],[1:78 80:100]))\b;

endvar = (eye(99,99)-P([1:78 80:100],[1:78 80:100]))\(2*msteps - b);
endvar = endvar - msteps.*msteps;

msteps = [msteps(1:78) ; 0 ; msteps(79:end) ; 0];
endvar = [endvar(1:78) ; 0 ; endvar(79:end) ; 0];

mcheck = sum(transpose(S))';
mcheck = [mcheck(1:78) ; 0 ; mcheck(79:end) ; 0];

figure(3);  clf
plot([0:100],msteps,'x');  hold on
plot([78 100],[0 0],'o')
title('\bf expected time to finish')
xlabel('\bf initial square (0-100)')
ylabel('\bf expected number of rolls')

%  plot snakes and ladders
for k=1:size(sinfo,1)
    plot([sinfo(k,1) sinfo(k,2)]-1, ...
        [msteps(sinfo(k,1)) msteps(sinfo(k,2))],'r')
    plot([sinfo(k,1) sinfo(k,1)]-1,[msteps(sinfo(k,1)) 0],'r:')
end

for k=1:size(linfo,1)
    plot([linfo(k,1) linfo(k,2)]-1, ...
        [msteps(linfo(k,1)) msteps(linfo(k,2))],'g')
    plot([linfo(k,1) linfo(k,1)]-1,[msteps(linfo(k,1)) 0],'g:')
end

%  transience time
figure(4);  clf
subplot(2,1,1)
plot([0:99],[S(1,1:78) nan S(1,79:end)],'x');  hold on
plot(sinfo(:,2)-1,S(1,sinfo(:,2)-(sinfo(:,2)>79)),'ro')
plot(linfo(1:end-1,2)-1,S(1,linfo(1:end-1,2)-(linfo(1:end-1,2)>79)),'go')
title('\bf transience time:  how long am i there?')
xlabel('\bf square # (start square is 0, red are snake tails, green are ladder tops)')
ylabel('\bf mean occupancy time (in turns)')

figure(5);  clf;  hold on
msteps1 = [msteps ; 0];
sq = 0:99;
sq = 1 + sq + (((-1).^floor(sq/10))==-1).*(9 - 2*mod(sq,10)); 
axis([0 10 0 10])

for k1 = 0:9
    for k = 1:10
        hh = text(k-0.5,k1+0.3,[num2str(sq(10*k1+k),3)]);
        set(hh,'color','r')
        text(k-1.0,k1+0.7,['\bf ' num2str(msteps1(sq(10*k1+k)+1),4)])
    end
end

set(gca,'xticklabel',' ','yticklabel',' ')
title(['\bf expected rolls to finish (start = ' num2str(msteps(1),4) ')'])
grid on

figure(4);
subplot(2,1,2);  hold on
plot(0:100,msteps,'o')
plot(0:100,sqrt(endvar),'gx')
title(['\bf mean (o) & \surd{variance} (x)'])
xlabel(['\bf square #'])