"""
craps.py

A craps statistical simulation.

-- aside to my python students --

 An example for the intro programming with python course.

 Note that this file includes :
  * a descripription at the top of the file,
    including author and date
  * functions, each with a docstring and clearly defined inputs/outputs
  * tests for some functions as appropriate
  * implentation details in comments within functions
  * an example session of what it looks like when it runs

 Your midterm project should also have these pieces, perhaps in
 several files.  For example, if your program is a graphical one,
 including screen captures of what it looks like would be helpful.

 There are a few python syntax constructions here which
 may be new to you; if so, check the python documentation
 when reviewing this. In particular, these may be new :
  i += 1           # add one to i; same as "i = i + 1"
  value in [7,11]  # boolean expression, same as "value==7 or value==11"

-- about this program --

For the rules of craps, see for example http://en.wikipedia.org/wiki/Craps

This code implements
 (a) one crap game,
 (b) a simulation of m craps games, and a count of wins/m_games (i.e. 45/100)
 (c) a series of n such simulations (i.e. [45/100, 48/100, ...])
 (d) a statistical analysis of the series which finds
     an estimate of the odds of winning at craps and its uncertainty.

Two globals, debug_print and interactive, control whether
the user is asked for input, and how much output is generated.

Running the program with debug_print=interactive=False
on a 2.66GHz intel core i7 mac laptop (and using the
unix "time" function to see how long it takes) gives :

  $ time python craps.py 
  --- Estimated odds of winning at craps ---
    Running 100 series of 10000 game simulations ...
    done; winning probabilities are [0.487, 0.489, ..., 0.499]
    Odds of winning = 0.49313 +- 0.00109  with 95% confidence.
  real	0m13.844s
  user	0m13.794s
  sys	0m0.041s

The true odds (see ./craps_analysis.txt) are 244/495 = 0.49292929...
which puts this result of 0.493 +- 0.001 right on target.

 Jim Mahoney | Oct 2011 | GPL
"""

from random import randrange 
from doctest import testmod
from math import sqrt

debug_print = False
interactive = False


def d6(n):  
  """ Roll n 6-sided dice and return the sum. """
  sum = 0
  for i in range(n):
    sum = sum + randrange(1,7)
  return sum


def test_dice(dice, n=1, low=1, high=6):
  """ Return true if many repititions of d6(n) is between low and high.
      >>> test_dice(d6, n=1, low=1, high=6)  # one die, sum is from 1 to 6
      True
      >>> test_dice(d6, n=2, low=2, high=12) # two dice, sum is 2 to 12
      True
  """
  for i in range(1000):
    roll = dice(n)
    if roll < low:
      return False
    if roll > high:
      return False
  return True


def craps():
  """ Play one game of craps. Return True if won, False if lost. """
  nrolls = 1
  first_roll = d6(2)
  if debug_print:
    print "  playing craps: "
    print "    roll %i : %i " % (nrolls, first_roll)
  if first_roll in [7, 11]:
    if debug_print:
      print "    win on roll %i" % nrolls
    return True
  elif first_roll in [2, 3, 12]:
    if debug_print:
      print "    lose on roll %i" % nrolls
    return False
  else:
    while True:
      nrolls += 1
      next_roll = d6(2)
      if debug_print:
        print "    roll %i : %i " % (nrolls, next_roll)
      if next_roll == first_roll:
        if debug_print: print "    win on roll %i" % nrolls
        return True
      elif next_roll == 7:
        if debug_print: print "    lose on roll %i" % nrolls
        return False


def count_wins(n_games):
  """ Return number of wins when playing n_games of craps. """
  wins = 0
  for i in range(n_games):
    if craps():
      wins = wins + 1
  return wins


def craps_series(n_series, n_games):
  """ Return a list of n_series of probabilities from count_wins(n_games). """
  result = [0] * n_series     # Create array for results, e.g. [0, 0, ... ]
  for i in range(n_series):
    result[i] = float(count_wins(n_games))/n_games
  return result


def average(numbers):
  """ Return the mean of a list of numbers.
      >>> average([1,2,3,4])
      2.5
  """
  return float(sum(numbers))/len(numbers)


def standard_deviation(numbers, s=True):
  """ Return one of two common definitions of the
      standard deviation (S.D.) of a given list of numbers :
        s=False : return sigma = sqrt(average((numbers - mean)**2))
        s=True  : return s = sigma * sqrt(n)/sqrt(n-1)
      The difference is that sigma is the S.D. of those numbers
      as a complete population, while s estimates the S.D. of
      a parent population that the numbers were drawn from.
      >>> "%.4f" % standard_deviation([1,2,3,4])
      '1.2910'
      >>> "%.4f" % standard_deviation([1,2,3,4], s=False)
      '1.1180'
  """
  mean = average(numbers)
  n = len(numbers)
  squared_diffs = [0] * n
  for i in range(n):
    squared_diffs[i] = (numbers[i] - mean)**2
  sigma2 = average(squared_diffs)
  if s:
    return sqrt(sigma2 * n / (n-1))
  else:
    return sqrt(sigma2)


def string_summary(numbers):
  """ Return string summary of list of numbers.
      >>> string_summary([2.3345, 0.11113, 3, 4, 5, 16.22289])
      '[2.334, 0.111, ..., 16.223]'
  """
  if len(numbers) < 3:
    return str(numbers)
  else:
    return "[%.3f, %.3f, ..., %.3f]" % (numbers[0], numbers[1], numbers[-1])


def main():
  print "--- Estimated odds of winning at craps ---"

  # Set number of games one simulation, and number of series to simulate.
  if interactive:
    try:
      n_games = input("  Number of games? (10) ")
    except:
      n_games = 10
    try:
      n_series = input("  Number of series? (5) ")
    except:
      n_series = 5
  else:
    n_games = 10000
    n_series = 100

  # Show results from one simulated series of games.
  if interactive:
    print "  Simulating %i games..." % n_games
    n_wins = count_wins(n_games)
    print "  done; number of simulated wins = %i" % n_wins
    print "  Odds of winning = %.2f" % (float(n_wins)/n_games)
    print

  # Compile statistics from many simulated series.
  print "  Running %i series of %i game simulations ..." % (n_series, n_games)
  results = craps_series(n_series, n_games)
  print "  done; winning probabilities are " + string_summary(results)

  mean = average(results)
  sd_population = standard_deviation(results)
  sd_mean = sd_population/sqrt(n_series)
  # sd_population is the standard deviation of any one series, 
  # that is, the standard deviation of the numbers in the results list.
  # But to estimate the uncertainty in our mean, we want the standard
  # deviation of the mean, namely how much that would likely change
  # if we were to do run the whole program again. A fundamental
  # result from statistics is that this sd_mean is
  # 1/sqrt(len(results)) smaller.

  print "  Odds of winning = %.5f +- %.5f  with 95%% confidence." \
      % (mean, 2*sd_mean)


if __name__ == "__main__":
  testmod()                # run the tests in the docstrings
  main()