Attempting to run multiple simulations of the Gillespie algorithm for a set of stochastic chemical reactions, in Python, in less than 10 minutes
Question: I have written Python code that generates a plot. In the code below, when I
set `maxtime = 0.1`, the program takes ~50s and when I set `maxtime = 1`, it
takes ~420s. I need to run it for `maxtime = 1000`.

I am familiar with Python syntax and writing "Matlabic" code, but am lost in
writing natively "Pythonic" code. As a result, I need help in **optimizing
this code for runtime, specifically in the two outer`for` loops and inner
`while` loop.**

  1. How can I make the code suitable for use with Numba or Cython? 
  2. If that's not possible, do I need to use functions, or `map`, or `lambda` statements?

Unfortunately, my Spyder IDE for Python is freezing up everytime I try to
profile the code. I would include those details if I could!

The code is below. Thanks.

    
    
    import numpy as np
    import matplotlib.pyplot as plt
    import pylab as pl
    import math
    
    maxtime = 1
    Qpvec = np.logspace(-2,1,101,10)
    lenQp = len(Qpvec)
    delta = 1e-3
    P0vec = Qpvec/delta
    SimNum = 100
    PTotStoch = np.zeros((SimNum,lenQp))
    k_minus = 1e-3
    k_cat = 1-1e-3
    k_plus = 1e-3
    zeta = 1e-4
    D0 = 10000
    kappa_M = (k_cat+k_minus)/zeta
    QpDeg = np.true_divide(1000*D0*Qpvec*k_cat,1000*Qpvec + kappa_M)
    
    for lenQpInd in range(0,lenQp):
        for SimNumInd in range(0,SimNum):
            Qp = Qpvec[lenQpInd]
            P0 = P0vec[lenQpInd]
            DP0 = 0
            P = math.floor(P0)
            DP = DP0
            D = D0
            time = 0
    
            while time < maxtime:
                u_time = pl.rand()
                u_event = pl.rand()
                rates=np.array([Qp,zeta*P*D,k_minus*DP,k_cat*DP])
                PTot = P + DP
                kT = np.sum(rates)
                tot = np.cumsum(rates)
                deltaT = -np.log(1-u_time)/kT
                time += deltaT
                if time > maxtime:
                    PTotStoch[SimNumInd,lenQpInd] = PTot
                    break
                elif u_event*kT < tot[0]:
                    P += 1
                elif u_event*kT < tot[1]:
                    P -= 1
                    DP += 1
                    D -= 1
                elif u_event*kT < tot[2]:
                    P += 1
                    DP -= 1
                    D += 1
                elif u_event*kT < tot[3]:
                    DP -= 1
                    D += 1   
    
    PMean = PTotStoch.mean(axis=0)
    PStd = PTotStoch.std(axis=0)
    
    plt.figure(0)
    plt.plot(Qpvec,PMean,marker=".") 
    plt.errorbar(Qpvec,PMean,yerr = 2*PStd, xerr = None)
    plt.show()
    


Answer: You can try compiling with [PyPy](http://pypy.org/). It doesn't look like
matplotlib is supported so you may have to include this part in a separate
file.

