8
\$\begingroup\$

I am trying to implement the modified Next Reaction Method with time varying propensities as mentioned in sections IV and V of this paper. At one step in the process this expression must be evaluated: \$\int_t^{t+\Delta t\mu}a(X(t,s), s) \textrm{d}s - (P_k - T_k) = 0\$. In the end this evaluation should tell us \$t + \Delta t\mu\$.

The code turns out to be quite slow, and since I need to run a huge system (many \$a(X(t,s))\$'s), this slow code is not an option. Below is the code, and comments on what it does, and what I've tried to optimize so far.

Maybe I should try Cython, but I do not know how, and for which parts it would be most beneficial; tips are welcome. Maybe other packages for integration, and root finding would help. But again, I am not aware of any such packages, or whether it would help me greatly, and tips are very welcome.

In the real project, this code is in a subclass of some more general class of simulation methods. This means that much less arguments are passed around, due to making objects attributes of this class.

import numpy as np
from scipy import integrate, optimize
from functools import partial, update_wrapper

# In my original code, the propensity functions (propensity_1, propensity_2)
# are generated automatically from a matrix and a bunch of parameter vectors,
# that is why setup might seem redundant in this code snippet.
# In one futile effort I tried @jit as shown below, but this made the code 
# drastically slower. 
def wrapped_partial(func, *args, **kwargs):
    partial_func = partial(func, *args, **kwargs)
    update_wrapper(partial_func, func)
    return partial_func

#@jit
def sinputt(t, amplitude=5.0, frequency=0.2):
    return amplitude * (1 + np.sin(frequency * t)) + 1
#@jit
def propensity_1(state, t, k):
    return k * sinputt(t=t)
#@jit
def propensity_2(state, t, k):
    return k * state

# Function to be integrated
def functionieren(t, state, func):
    return func(t=t, state=state)
# Returns integral
def integrieren(t, dt, state, func):
    return np.abs(
        integrate.quad(func=functionieren, a=t, b=dt, args=(state, func), epsabs=1.e-4, epsrel=1.e-4, limit=50)[0])
# Function for which root is to be found
def optimizieren(dt, t, state, func, P_T):
    return P_T - integrieren(t, dt, state, func)
# Function that finds root, thus t+dt
def solve(t, dt, state, func, P_T):
    return optimize.newton(func=optimizieren,x0=dt,args=(t,state,func,P_T))
# One step in the algorithm
def run_step(P, T, state, t, propensities, funcs):
    dtk = np.divide((P - T), propensities)
    r_i = np.argmin(dtk)
    dt = dtk[r_i, 0]
    t += dt
    if r_i == 0:
        state += 1
    elif r_i == 1:
        state -= 1

    T += propensities * dt
    P[r_i, 0] -= np.log(np.random.random(1))

    prop_func1,prop_func2=funcs
    #HERE IS THE BUGGER
    dtp=solve(t=t,dt=t+1.0,state=state,func=prop_func1,P_T=P[0,0]-T[0,0])
    propensity1=prop_func1(state=state, t=dtp)
    propensity2=prop_func2(state=state,t=t)

    propensities = np.array([[propensity1],
                              [propensity2]])
    return P, T, state, t, propensities

#The algorithm, final time=100
def run_sim(P, T, state, t, propensities, funcs):
    state_track = [state0]
    t_track = [t0]
    while t < 100:
        P, T, state, t, propensities = run_step(P, T, state, t, propensities,funcs)
        state_track.append(np.copy(state))
        t_track.append(np.copy(t))
    return np.array(state_track), np.array(t_track)

# Parameters
k1 = 2.0
k2 = 2.0

# Initial conditions
t0 = 0
state0 = 5

# Initial conditions internal clocks
T0 = np.zeros((2, 1))
P0 = -np.log(np.random.random((2, 1)))

# Generation of propensity functions
prop_func1 = wrapped_partial(propensity_1, k=k1)
prop_func2 = wrapped_partial(propensity_2, k=k2)

# Initial propensities
propensity1=solve(t=t0,dt=t+1.0,state=state0,func=prop_func1,P_T=P0[0,0]-T0[0,0])
propensity2=prop_func2(state=state0,t=t0)
propensities0=np.array([[propensity1],
                       [propensity2]])

# If dividing by a 0 propensity inf is returned
np.seterr(divide="ignore")
state_track,t_track=run_sim(P0,T0,state0,t0,propensities0,(prop_func1,prop_func2))

# Plot of the results as they are intended
import matplotlib.pyplot as plt
t_sin = np.arange(0, 100, 0.1)
input_sin = sinputt(t=t_sin)
plt.plot(t_sin, input_sin,label='input')
plt.step(t_track, state_track,label='state')
plt.legend(bbox_to_anchor=(0.65, 1), loc=2, borderaxespad=0.)
plt.show()

# Uncomment the %timeit below to see which steps are slow
# %timeit run_step(P, T, state, t, propensities,(prop_func1,prop_func2))
# %timeit solve(t=t,dt=t+1.0,state=state,func=prop_func1,P_T=P[0,0]-T[0,0])
# %timeit integrieren(t, t+1.0, state, prop_func1)
# %timeit functionieren(t, state, prop_func1)

# Output timeit
# 1000 loops, best of 3: 294 µs per loop
# 1000 loops, best of 3: 265 µs per loop
# 10000 loops, best of 3: 42.7 µs per loop
# The slowest run took 7.98 times longer than the fastest. This could mean 
# that an intermediate result is being cached.
# 1000000 loops, best of 3: 1.74 µs per loop
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.