from CoaSim import *
import sys

if len(sys.argv) != 5:
    print 'Usage:', sys.argv[0], 'beta noSeqs mergeTime noSamples'
    sys.exit(2)

beta      = float(sys.argv[1])
noSeqs    = int(sys.argv[2])
mergeTime = float(sys.argv[3])
noSamples = int(sys.argv[4])

class callback(object):
    def __init__(self):
        self.count = 0
    def recombinationEvent(self, *dummy):
        self.count += 1

# mapping seqs*beta*mergeTime -> rho
rhoMap = {(20,   0, 0.150): 0.19,
          (50,   0, 0.150): 0.15,
          (20, 100, 0.015): 2.22,
          (50, 100, 0.015): 1.20}
rho = rhoMap[(noSeqs,beta,mergeTime)]

counts = dict()
for i in xrange(noSamples):
    cb = callback()
    arg = simulate([], noSeqs, beta=beta, rho=rho,
                   callbacks=cb,
                   keepEmptyIntervals=True)
    try:
        counts[cb.count] += 1
    except KeyError:
        counts[cb.count] = 1


observedCounts = counts.keys()
observedCounts.sort()
total = sum(counts.values())
maxFrac = max([float(counts[oc])/total for oc in observedCounts])
for oc in observedCounts:
    print '%d: %6.2f%%' % (oc, 100.0*counts[oc]/total),
    print '*'*int(100*maxFrac*counts[oc]/total)