coupons

#!/usr/bin/python

# Coupon Simulator 2015
# let's calculamate perbability! Simulates the "balls in boxes"
# problem with 50 boxes: how many throws, on average, until each box has at
# least 1 ball in it? Another form of the "coupon collector's problem".
# usage: coupons.py <# of trials> <-y> (to skip interactive approval of > 2500 trials)
# - lardlung

import random
import sys
import argparse
import math

results = list()

parser = argparse.ArgumentParser(description = "coupons.py: a probability simulator of the coupon collector's dilemma, or, balls and boxes.")
parser.add_argument("trialcount", help = "<#> number of trials to run", type = int)
parser.add_argument("-y", "--yes", action = 'store_true', help = "optional; skip interactive approval of trialcount > 2500")
parser.add_argument("-b","--bins", help = "optional; number of boxes/bins to use. Default is 50", default = 50, action='store')



# do one trial.

def one_run(bins):

    myBins     = list()
    throwCount = 0
    poop       = 1

    for i in range(0,bins):
        myBins.insert(i,0)

    while poop == 1:

        myThrow     = random.randint(0, bins - 1)
        myBins[myThrow] = myBins[myThrow] + 1
        throwCount  = throwCount + 1

        if 0 not in myBins:
            poop = 2
            return throwCount
        else:
            continue

# average results of multiple trials.

def do_results(trialCount):

    for i in range( 0, trialCount ):
        results.insert(i,0)
        results[i] = int(one_run(bins))

    trialsAverage = sum( results ) / float(trialCount)
    return trialsAverage

# find a mathematical approximation of what the results should be.
# em_gamma is an approximation of the euler-mascheroni constant!
# http://mathworld.wolfram.com/Euler-MascheroniConstantApproximations.html

def do_calc(bins):

    em_gamma = math.pi/(2*math.e)
    calcResult = bins * ( em_gamma + math.log(bins) )
    return calcResult

# calculate the percent difference between the approximation and the trial results.

def do_percent(calcResults,avgResults):
    x = ( abs(calcResults - avgResults) ) / ((calcResults + avgResults)/2)
    pDiff = 100 * float(x)
    return pDiff


# go main loop, go!

try:
    args = parser.parse_args()
    trialCount = args.trialcount

    if args.bins == False:
        bins = 50
    else:
        bins = int(args.bins)

    if trialCount >= 1 and trialCount <= 2500:

        avgResults = do_results(trialCount)
        calcResults = do_calc(bins)
        pDiff = do_percent(avgResults,calcResults)

        print "trials:", trialCount, "bins:", bins
        print "expected result(calculated): {0:.3f}".format(calcResults)
        print "averaged result:", avgResults
        print "percent difference: {0:.3f}%".format(pDiff)

    elif trialCount > 2500:

        if not args.yes:
            verify = raw_input("Are you SURE? y/n ")


            if verify in ('Y','y'):

                avgResults = do_results(trialCount)
                        calcResults = do_calc(bins)
                        pDiff = do_percent(avgResults,calcResults)

                        print "trials:", trialCount, "bins:", bins
                        print "expected result(calculated): {0:.3f}".format(calcResults)
                        print "averaged result:", avgResults
                        print "percent difference: {0:.3f}%".format(pDiff)


            elif verify in ('N','n'):
                print "aborted, but have 1 quick trial anyways:", one_run(bins)
                sys.exit(0)

            else:
                print "invalid input: aborted."
                sys.exit(0)

        elif args.yes:

            avgResults = do_results(trialCount)
                calcResults = do_calc(bins)
                pDiff = do_percent(avgResults,calcResults)

                    print "trials:", trialCount, "bins:", bins
                print "expected result(calculated): {0:.3f}".format(calcResults)
                    print "averaged result:", avgResults
                    print "percent difference: {0:.3f}%".format(pDiff)


        else:
            print "Usage: coupons.py <# of trials> <-y> (to skip interactive approval of > 2500 trials)"
            sys.exit(1)


    else:
        print "bombElse", results


except:
    print "bombExcept"
    sys.exit(0)