euler213

euler213v1.py

import numpy as np
# build transition matrix
transition_matrix = np.zeros((900,900))

# corner cases
corner_transition_idx = [(0, 1), (0, 30), (29, 28), (29, 59), (870, 871), (870, 840), (899, 898), (899, 869)]

# edge cases
top_idx = [(x, x+30) for x in range(1, 29)] + [(x, x-1) for x in range(1, 29)] + [(x, x+1) for x in range(1, 29)]
left_idx = [(x, x+30) for x in range(30, 870, 30)] + [(x, x-30) for x in range(30, 870, 30)] + [(x, x+1) for x in range(30, 870, 30)]
right_idx = [(x, x+30) for x in range(59, 899, 30)] + [(x, x-30) for x in range(59, 899, 30)] + [(x, x-1) for x in range(59, 899, 30)]
bottom_idx = [(x, x-30) for x in range(871, 899)] + [(x, x-1) for x in range(871, 899)] + [(x, x+1) for x in range(871, 899)]
edge_transition_idx = top_idx + left_idx + right_idx + bottom_idx

# everything else
interior_transition_idx = []
for x in range(31,869):
    if x % 30 in (0, 29):
        continue
    interior_transition_idx += [(x, x+1), (x, x-1), (x, x+30), (x, x-30)]

# combine
for x, y in corner_transition_idx:
    transition_matrix[x, y] = 0.5

for x, y in edge_transition_idx:
    transition_matrix[x, y] = 1/3

for x, y in interior_transition_idx:
    transition_matrix[x, y] = 0.25

# compute A^50
after_fifty = np.linalg.matrix_power(transition_matrix, 50)

# get probability of no flea
no_flea_A = 1 - after_fifty

# multiply across rows
per_cell_no_flea = np.prod(no_flea_A, axis=0)

# sum all to final expectation value
np.sum(per_cell_no_flea)

euler213v2.py

import numpy as np

# build transition matrix
k1_diag = np.ones(899) * 0.25
k1_diag[np.arange(1, 30) * 30 - 1] = 0
k1_diag[np.arange(1, 30) * 30] = 1/3
k1_diag[[0, 870]] = 0.5
k1_diag[1:29] = 1/3
k1_diag[871:] = 1/3

k30_diag = np.ones(870) * 0.25
k30_diag[np.arange(1, 30) * 30 - 1] = 1/3
k30_diag[np.arange(1, 29) * 30] = 1/3
k30_diag[[0, 29]] = 0.5
k30_diag[1:29] = 1/3

tm = np.diag(k1_diag, k=1) + np.diag(k1_diag[::-1], k=-1) + np.diag(k30_diag, k=30) + np.diag(k30_diag[::-1], k=-30)

# compute expectation value
np.sum(np.prod(1 - np.linalg.matrix_power(tm, 50), axis=0))