# Written by Robert Richardson for Stat 624
# Python script for simulating Chutes and Ladders

import numpy as np              # numerical library
import pandas as pd             # data wrangling library
import matplotlib.pyplot as plt # plotting library
from scipy.stats import norm 
import random

random.seed(12)

# Set up 100 spaces on a board
board = np.linspace(1,100,100)

# Randomly sample 40 locations for chutes and ladders
events = np.random.choice(board, 40, replace=False)
chutes = events[0:20]
ladders = events[20:40]

chutes_start = chutes.reshape(10,2).max(axis=1)
chutes_end = chutes.reshape(10,2).min(axis=1)

ladders_start = ladders.reshape(10,2).min(axis=1)
ladders_end = ladders.reshape(10,2).max(axis=1)

# Simulating one turn
def turn(Player):
  tempPos = Player + np.random.choice(6,1)+1
  if tempPos in chutes_start:
    ch = np.where(chutes_start == tempPos)
    tempPos = chutes_end[ch]
  if tempPos in ladders_start:
    ch = np.where(ladders_start == tempPos)
    tempPos = ladders_end[ch]
  return tempPos

# Simulate 1000 games with one player
nrun = 1000
nturn = [0]*nrun
for i in range(nrun):
  k=0
  Player_pos = 0
  while Player_pos < 100:
    Player_pos = turn(Player_pos)
    k = k+1
  nturn[i] = k

## Print out the monte carlo estimate of the number of turns with monte carlo error
est = np.mean(nturn)
print("The Monte Carlo Estimate of the number of turns is %4.3f" % est)
ci = est + np.array([-1.,1.])*norm.ppf(.975)*np.std(nturn)/np.sqrt(nrun)
print("The 95%% confidence interval is (%4.3f,%4.3f)." % (ci[0],ci[1]))

# Plot number of turns
n, bins, patches = plt.hist(x=nturn, bins='auto', color='#0504aa',alpha=0.7, rwidth=0.85)
# plt.show()
plt.savefig("hist.pdf")

# 2 player game
def play2():
  win = 0
  Pl1 = 0
  Pl2 = 0
  while Pl1 < 100 and Pl2 < 100:
    Pl1 = turn(Pl1)
    Pl2 = turn(Pl2)
  if Pl1 >= 100:
    win = 1
  return win

n1_win = [play2() for i in range(nrun)]

est2 = np.mean(n1_win)
print("The Monte Carlo Estimate of the number of times \n the first player wins is %4.3f" % est2)

ci_a = est2 + np.array([-1.,1.])*norm.ppf(.975)*np.std(n1_win)/np.sqrt(nrun)
print("The 95%% confidence interval is (%4.3f,%4.3f)." % (ci_a[0],ci_a[1]))

ci_b = est2 + np.array([-1.,1.])*norm.ppf(.975)*np.sqrt(est2*(1-est2))/np.sqrt(nrun)
print("The 95%% confidence interval is (%4.3f,%4.3f)." % (ci_b[0],ci_b[1]))