#!/usr/bin/env python

import bisect
import math
import random
import sys

class Poisson_Process:
    def __init__(self, par, plimit = 1E-9):
        assert(par > 0 and par < 1)
        if par > 0.1:
            sys.stderr.write('Warning, poisson parameter dangerously high: %f\n' %(float(par)))
        self.par = float(par)
        assert(plimit > 0 and plimit < 1)
        self.plimit = plimit
        self.gen_prob_table()

    # public methods
    def random(self):
        r = random.uniform(0, 1)
        # Use binary search. self.probs[] is a cumulative discrete probability
        # table. The first entry is the probability of (self.n - 1)
        # events occuring. The kth entry (starting from 0) is the probability
        # of self.n-k-1 events occuring added with entry k-1. The last
        # value in table is fixed to 1.0 and r < 1.0 always, so
        # the maximum value returned by binary search is self.n - 1.
        i = bisect.bisect_left(self.probs, r)
        return self.n - i

    # private methods
    def poisson_probability(self, n):
        # See http://en.wikipedia.org/wiki/Poisson_distribution
        f = 1
        i = 2
        while i < n:
            f = f * i
            i += 1
        return pow(self.par, n) / f * math.exp(-self.par)

    def gen_prob_table(self):
        self.probs = []
        n = 0
        while True:
            p = self.poisson_probability(n)
            self.probs.append(p)
            if p < self.plimit:
                break
            n += 1
        assert(n > 1)
        # hack, reverse it for proper sorting order for binary search in
        # randomize_n()
        self.probs.reverse()
        self.n = n

        sys.stderr.write('Poisson cumulative probabilities: ')
        s = 0.0
        for i in xrange(len(self.probs[:-1])):
            self.probs[i] += s
            s = self.probs[i]
            sys.stderr.write('%f ' %(s))
        sys.stderr.write('\n')

        assert(s < 1)
        self.probs[-1] = 1.0

        sys.stderr.write('Poisson process: at most %d occurences / invocation\n' %(n - 1))
        sys.stderr.write('Poisson probability for 0 occurences: %f\n' %(1 - self.probs[-2]))


def poissontest():
    p = 0.05
    n = 1000000
    x = Poisson_Process(p)
    s = 0
    for i in xrange(n):
        y = x.random()
        print y
        s += y
    print "total", s
    print "error", float(s)/(p*n) - 1

poissontest()