How to sample Poisson processes in Python faster if the rate is non-constant?

I'm sampling a Poisson process at a millisecond time scale where the rate is not fixed. I discretise the sampling process by checking in each interval of size delta whether there is an event there or not based on the average rate in that interval. Since I'm using Python it's running a bit slower than I would hope it to be. The code I'm currently using is the following:

import numpy
def generate_times(rate_function,max_t,delta):
    times = []
    for t in numpy.arange(delta,max_t,delta):
        avg_rate = (rate_function(t)+rate_function(t+delta))/2.0
        if numpy.random.binomial(1,1-math.exp(-avg_rate*delta/1000.0))>0:
            spikes.extend([t])
    return times

The rate function can be arbitrary, I'm not looking for a closed form solution given a rate function.

If you want some parameters to play with you can try:

max_t = 1000.0
delta = 0.1
high_rate = 100.0
low_rate = 0.0
phase_length = 25.0
rate_function = (lambda x: low_rate + (high_rate-low_rate)*0.5*(1+math.sin(2*math.pi*x/phase_length)))