### Frogs, flies, and Julia

Here's the July PuzzlOR question, courtesy of John Toczek. Now that the competition is over, let's go to town. A frog is looking to catch his next meal just as a fly wanders into his pond. The frog jumps randomly from one lily pad to the next in hopes of catching the fly. The fly is unaware of the frog and is moving randomly from one red flower to another.

The frog can only move on the lily pads and the fly can only move on the flowers. The interval at which both the frog and the fly move to a new space is one second. They never sit still and always move away from the space they are currently on. Both the frog and the fly have an equal chance of moving to any nearby space including diagonals. For example, if the frog were on space A1, he would have a 1/3 chance each of moving to A2, B2, and B1.

The frog will capture the fly when he lands on the same space as the fly.

Question: Which space is the frog most likely to catch the fly?

This is another Markov chain puzzle, similiar in some ways to the post about modeling chutes and ladders as a Markov chain or the spy catcher puzzle. The twist here is that instead of being a static absorbing Markov chain or a single random process, the game ends when two discrete random processes converge on the same state.

Lately, I've been playing with the Julia language, so I'll do this one using Julia. Thanks to the awesome IJulia library, we don't have to give up the IPython notebook. For me, that's a huge selling point.

## Setting up the problem¶

Let's number these squares: Now that we have consistent way to refer to the states that the fly and the frog can be in, we will create adjacency matrices that define how they are able to move.

By convention, element $M_{ij}$ of an adjacency matrix tells you if state $j$ (the column) is reachable from state $i$ (the row). For example, if the fly is on square 9 then the only two adjacent nodes with flowers that it can reach next are squares 6 and 8. Therefore, in row 9 of the flower adjacency matrix, only the 6th and 8th elements are equal to one (indicating adjacency) while the rest are equal to zero (meaning they are unreachable from the current state).

In :
# the frog's possible transitions
adj_frog =  [0  1  0  1  1  0  0  0  0;
1  0  0  1  1  0  0  0  0;
0  0  0  0  0  0  0  0  0;
1  1  0  0  1  0  1  1  0;
1  1  0  1  0  0  1  1  0;
0  0  0  0  0  0  0  0  0;
0  0  0  1  1  0  0  1  0;
0  0  0  1  1  0  1  0  0;
0  0  0  0  0  0  0  0  0]

adj_fly = [0  0  0  0  0  0  0  0  0;
0  0  0  0  0  0  0  0  0;
0  0  0  0  0  1  0  0  0;
0  0  0  0  0  0  1  0  0;
0  0  0  0  0  0  0  0  0;
0  0  0  0  0  0  0  1  1;
0  0  0  1  0  0  0  1  0;
0  0  0  1  0  1  1  0  1;
0  0  0  0  0  1  0  1  0]

Out:
9x9 Array{Int64,2}:
0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0
0  0  0  0  0  1  0  0  0
0  0  0  0  0  0  1  0  0
0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  1  1
0  0  0  1  0  0  0  1  0
0  0  0  1  0  1  1  0  1
0  0  0  0  0  1  0  1  0

We'll also define a function to normalize the rows, and turn them into conditional probabilities. (Much more on this in the spy catcher post.)

In :
function stochastic_matrix_from_adjacency_matrix(M)

# initialize a Float64 matrix of same size as M
P = zeros(size(M))

# normalize each non-zero row by dividing by sum of row
for i in 1:size(P, 2)
if sum(M[i, :]) > 0
P[i, :] = M[i, :] / sum(M[i, :])
end
end

P
end

Out:
stochastic_matrix_from_adjacency_matrix (generic function with 1 method)

Using this function, we'll get the stochastic matrices for the fly and frog.

In :
P_fly = stochastic_matrix_from_adjacency_matrix(adj_fly)

Out:
9x9 Array{Float64,2}:
0.0  0.0  0.0  0.0   0.0  0.0   0.0   0.0  0.0
0.0  0.0  0.0  0.0   0.0  0.0   0.0   0.0  0.0
0.0  0.0  0.0  0.0   0.0  1.0   0.0   0.0  0.0
0.0  0.0  0.0  0.0   0.0  0.0   1.0   0.0  0.0
0.0  0.0  0.0  0.0   0.0  0.0   0.0   0.0  0.0
0.0  0.0  0.0  0.0   0.0  0.0   0.0   0.5  0.5
0.0  0.0  0.0  0.5   0.0  0.0   0.0   0.5  0.0
0.0  0.0  0.0  0.25  0.0  0.25  0.25  0.0  0.25
0.0  0.0  0.0  0.0   0.0  0.5   0.0   0.5  0.0 
In :
P_frog = stochastic_matrix_from_adjacency_matrix(adj_frog)

Out:
9x9 Array{Float64,2}:
0.0       0.333333  0.0  0.333333  0.333333  0.0  0.0       0.0       0.0
0.333333  0.0       0.0  0.333333  0.333333  0.0  0.0       0.0       0.0
0.0       0.0       0.0  0.0       0.0       0.0  0.0       0.0       0.0
0.2       0.2       0.0  0.0       0.2       0.0  0.2       0.2       0.0
0.2       0.2       0.0  0.2       0.0       0.0  0.2       0.2       0.0
0.0       0.0       0.0  0.0       0.0       0.0  0.0       0.0       0.0
0.0       0.0       0.0  0.333333  0.333333  0.0  0.0       0.333333  0.0
0.0       0.0       0.0  0.333333  0.333333  0.0  0.333333  0.0       0.0
0.0       0.0       0.0  0.0       0.0       0.0  0.0       0.0       0.0

Now let's write some functions to simulate the frog and fly game.

In :
function move(state, P):
"""
Given a current state and a transition matrix, execute one
transition.

state -- an int representing the current state
P -- the stochastic matrix
"""

# find stochastic vector given current state
probs = P[state, :]

# randomly pick one of these possibilities proportional
# to its probability of being selected
findfirst(cumsum(probs[1:end]) .>= rand())
end

function simulate(n_simulations)
"""
Run a given number of fly/frog game simulations.

n_simulations -- number of simulations to conduct
"""

# initialize the output array
result = zeros(Int64, n_simulations, 2)

for i in 1:n_simulations

# set initial states from the problem definition
fly = 3
frog = 5

# initialize a counter for the number of moves
n_moves = 0

# while the fly hasn't been eaten ...
while fly != frog

# execute one discrete tick
fly = move(fly, P_fly)
frog = move(frog, P_frog)

# increment the move counter
n_moves += 1

end

# fill in the results from the game
result[i, :] = [n_moves, fly]

end

# return the number of moves and location
result
end

Out:
simulate (generic function with 1 method)

Having set up all the necessary machinery, we'll run 50,000 iterations of the simulation.

In :
using DataFrames

n_simulations = 50000
result = simulate(n_simulations)
df = DataFrame(moves=result[:, 1], end_square=result[:, 2])

Out:
movesend_square
1224
247
3104
497
5138
634
737
847
928
10194
1157
12138
13198
14127
1567
16374
1777
18128
1948
2074
21178
2234
2344
2448
25214
2628
27124
28134
2968
3058

## How many moves does the game last?¶

Here are some descriptive statistics for the number of moves:

In :
describe(df[:moves])

Min      2.0
1st Qu.  5.0
Median   9.0
Mean     12.22894
3rd Qu.  16.0
Max      108.0
NAs      0
NA%      0.0%


So on average, it took around 12 moves for the fly to get eaten. The fact that it takes so long may be a little less surprising when we consider that the fly has to random-walk itself out of the top-right corner down onto the lily pads before even being reachable by the frog.

We can also look at a histogram of the number of moves:

In :
using Gadfly

plot(df, x="moves", Geom.histogram(bincount=30),
xintercept=[mean(df[:moves])], Geom.vline(color="red"))

Out:

Some of the games lasted a surprisingly long time for such a tiny 3x3 board.

## What was the most common end square?¶

In :
_, count = hist(df[:end_square], 0:9)
end_square_counts = DataFrame(count=count, frequency=count./sum(count))

Out:
countfrequency
100.0
200.0
300.0
4174800.3496
500.0
600.0
7111980.22396
8213220.42644
900.0
In :
plot(end_square_counts, x=1:9, y="frequency", color=1:9, Geom.bar,
Guide.xticks(ticks=[0:1:10]),
Guide.xlabel("end square"))

Out:

So we have our answer, which is that square number 8 was where the frog most often catches the fly. Previous: Using Python for linear optimization

Any comments or suggestions? Let me know.