using LinearAlgebra, PyPlot

M = zeros(101,101)
for i = 2:100
    M[i,max(1,i-6):(i-1)] .= 1/6
end
# last row
for i = 1:6
    M[101,101-i] = (7-i)/6 # = 6/6, 5/6, 4/6, ..., 1/6
end
M[101,101] = 1 # once we get to the last space, we stay there
M

e₁ = zeros(101); e₁[1] = 1
M*e₁

M^2 * e₁

sum(M, dims=1) # sum M along the first dimension, i.e. sum Mᵢⱼ over i, i.e. sum each column

eigvals(M)

λ, X = eigen(M)
det(X)

X[:,end]

# Plot the probabilities on something like a chutes and ladders board.  We won't show the starting position (0)
# since that is not on the board.
function plotchutes(p)
    P = transpose(reshape(p[2:101], 10,10)) # reshape to a 10×10 array and transpose to row-major
    # every other row should be reversed, corresponding to how players "zig-zag" across the board:
    for i = 2:2:10
        P[i,:] = reverse(P[i,:])
    end
    imshow(P, aspect="equal", cmap="Reds", norm=PyPlot.matplotlib["colors"]["LogNorm"](vmin=1e-3, vmax=1), origin="lower", interpolation="none")
    colorbar(label="probability")
    xticks(-0.5:1:10, visible=false)
    yticks(-0.5:1:10, visible=false)
    grid()
    for i = 1:10, j = 1:10
        n = (i-1)*10 + j
        x = iseven(i) ? 10-j : j-1
        y = i-1
        text(x,y, "$n", horizontalalignment="center", verticalalignment="center")
    end
end

fig = figure()
# using Interact
# @manipulate for n in slider(1:100, value=1)  
for n in (1,2,3,4,5,10,20,30,40,50,100)
    display(
    withfig(fig) do
        plotchutes(M^n*e₁)
        title("distribution after $n moves")
    end
    )
end

M^100*e₁

plot(0:100, [(M^n * e₁)[end] for n = 0:100], "b.-")
xlabel("number of moves n")
ylabel("probability of finishing in ≤ n moves")
grid()
title("boring chutes & ladders (no chutes or ladders)")

plot(1:100, diff([(M^n * e₁)[end] for n = 0:100]), "b.-")
xlabel("number of moves n")
ylabel("probability of finishing in n moves")
grid()
title("boring chutes & ladders (no chutes or ladders)")

sum((1:100) .* diff([(M^n * e₁)[end] for n = 0:100]))

T = zeros(101,101)

for t in (1=>39, 4=>14, 9=>31, 28=>84, 36=>44, 51=>67, 80=>100, 71=>91, # ladders
          16=>6, 47=>26, 49=>11, 56=>53, 64=>60, 92=>73, 95=>75, 98=>78) # chutes
    T[t[2]+1,t[1]+1] = 1
end

# Set T[j,j] = 1 in spaces with no chute/ladder
for j = 1:101
    if all(T[:,j] .== 0)
        T[j,j] = 1
    end
end

sum(T,dims=1)

sum(T*M, dims=1)

plotchutes(T*M*e₁)
title("probability distribution after 1 move")

fig = figure()
# using Interact
# @manipulate for n in slider(1:100, value=1)
for n in (1,2,3,4,5,10,20,30,40,50,100)
    display(
    withfig(fig) do
        plotchutes((T*M)^n*e₁)
        title("distribution after $n moves")
    end
    )
end

plot(1:100, diff([((T*M)^n * e₁)[end] for n = 0:100]), "b.-")
plot(1:100, diff([(M^n * e₁)[end] for n = 0:100]), "r.--")
xlabel("number of moves n")
ylabel("probability of finishing in n moves")
grid()
title("number of moves to finish chutes & ladders")
legend(["chutes & ladders", "no chutes/ladders"])

sum((1:1000) .* diff([((T*M)^n * e₁)[end] for n = 0:1000]))

semilogy(0:350, [1-((T*M)^n * e₁)[end] for n = 0:350], "b.-")
xlabel("number of moves n")
ylabel("probability of NOT finishing in n moves")
grid()
title("chance of long chutes & ladders game")

A = (T*M)'[1:100,1:100]  # the 100x100 upper-left corner of (TM)ᵀ
N = inv(I - A)           # N[i,j] = expected number of visits to i starting at j
(N * ones(100))[1]       # expected number of moves to finish starting at 1