Monte Carlo modeling of dice-based treasure trails

trails.py

 #!/usr/bin/python3

# Monte Carlo modeling of dice-based "treasure trails"
# Ian McDougall, Jan 2023

import math
import random
import numpy as np
import numpy.ma as ma
import matplotlib as mpl
import matplotlib.pyplot as plt

def trial(runs, stop, dice):
	# Initialization
	step = ma.zeros(runs, dtype=int)
	clue = ma.zeros(runs, dtype=int)
	tally = np.zeros(stop, dtype=int)
	steps = np.zeros(stop, dtype=int)
	final = np.zeros(stop + dice, dtype=int)
	# hardmask stops assignment from unmasking, but does not stop assignment itself!
	step.harden_mask()
	clue.harden_mask()
	while ma.min(clue) < stop:
		roll = np.array([random.randint(1, dice) for i in clue]) + step
		clue[~ma.getmaskarray(clue)] = roll[~ma.getmaskarray(step)]
		clue = ma.masked_greater_equal(clue, stop)
		step.mask = ma.getmask(clue)
		for i in ma.compressed(clue): # "compressed" returns only non-masked data
			tally[int(i)] += 1
		step += 1 # masking seems to work for this kind of assignment
	step.soften_mask()
	clue.soften_mask()
	step.mask = ma.nomask
	clue.mask = ma.nomask
	for i in step:
		steps[int(i)] += 1
	for i in clue:
		final[int(i)] += 1
	return(tally, steps, final)

# the "pythonic" way to represent summaries is with "collections," but this seems like a hassle

def summary_mean(inarray):
	mean = sum((range(len(inarray)) * inarray)) / sum(inarray)
	return mean

def summary_var(inarray):
	mean = summay_mean(inarray)
	variance = sum(np.square(range(len(inarray)) - mean) * inarray)
	return variance

def summary_mode(inarray):
	mode = list(inarray).index(max(inarray))
	return mode

def summary_norm(inarray):
	outarray = inarray / sum(inarray)
	return outarray

size = 100000

(tally, steps, final) = trial(size, 7, 6)
mean_len = summary_mean(steps)
print("mean length of trail:", mean_len)
loot = final.copy()
np.put(loot, 7, 0)
print("there were ", sum(loot), "loot rolls, and ", final[7], "intel rolls.")
print("mean loot roll:", summary_mean(loot))
print("the most common clue given was type:", list(tally).index(max(tally)))
print("the least common clue type given was:", list(tally).index(min(tally[1:])))

print(steps)
print(final)
print(tally)

fig1 = plt.figure()
plt.bar(range(len(tally))[1:], tally[1:])
plt.bar(range(len(final))[1:], final[1:])
plt.xticks(range(len(final)), range(len(final)))
plt.ylabel('Result Frequency')
plt.title('Results of 1d6+n, Stopping on 7+, n={}'.format(size))
fig1.savefig('figure1.png', bbox_inches='tight')

stops = [5, 6, 7, 8, 9]
sides = [4, 6]
size = 50000

nrows = len(stops)
ncols = len(sides)
maxval = max(stops) + np.asarray(sides)

fig2, ax2 = plt.subplots(nrows=len(stops), ncols=len(sides), sharex='col', sharey=False )
for i in range(nrows):
	for j in range(ncols):
		#axn = plt.subplot(nrows, ncols, (i*ncols + j + 1))
		(tally, steps, final) = trial(size, stops[i], sides[j])
		ax2[i,j].bar(range(len(tally))[1:], tally[1:])
		ax2[i,j].bar(range(len(final))[1:], final[1:])
		ax2[i,j].set_xticks(range(maxval[j])[1:])
		ax2[i,j].set_xticklabels([])
		ax2[i,j].set_yticks([])
		nbar = summary_mean(steps)
		ax2[i,j].text(1, 1, '$\overline{n}$=' + "%.1f" % nbar, ha='right', va='top', transform=ax2[i,j].transAxes)
		ax2[i,j].set_frame_on(False)
#fig2.xticklabels(["Roll 1d{}+1".format(a) for a in sides])
for i in range(nrows):
	ax2[i,0].set_ylabel("{}+".format(stops[i]))
for j in range(ncols):
	ax2[0,j].set_xlabel("Roll 1d{}".format(sides[j]))
	ax2[0,j].xaxis.set_label_position('top')
fig2.tight_layout()
fig2.savefig('figure2.png',bbox_inches='tight')