# Root finding and integration [closed]

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.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

• I tried running the code, and it failed with a NameError. – Peter Taylor Aug 2 '19 at 10:05

• Spacing Between variables and operators, there should be spaces to improve readability. Any code like prop_func1,prop_func2=funcs should instead be prop_func1, prop_func2 = funcs.
• Lots of Parameters: When you have tons of parameters you are passing/receiving to/in a method, having them all on one line can decrease readability significantly. You should put each parameter on its own line, to make it clear what is being passed/received.
• Docstrings: You should have a docstring at the beginning of every method, class, and module you write. This helps documentation determine what your code is supposed to accomplish. You had the right idea with comments above the methods explaining what they do, now just move them inside the function between """ docstring here """.
• Import Statements: All import statements should be at the top of the program, regardless of when you use them in your code. My personal preference is to have them ordered alphabetically, but it's up to you. The convention is to have them at the top of the program. (PEP-8 import statement info here)
• Main Guard: Let's say you want to reuse this system for another program. You import the file, write your code, and WHAT??? Everything is messed up. Having a main guard is a way to prevent this from happening. Just include everything outside of methods in a if __name__ == '__main__' statement, and once you import the file, since it's not the main file being run, that code won't run. An answer from StackOverflow better explains this.
• Constants: Any and all constants you have in your program should be UPPERCASE, to make it clear they are constants.

Updated Code

"""
Program implements the modified Next Reaction Method with time varying propensities
"""

import matplotlib.pyplot as plt
import numpy as np

from functools import partial, update_wrapper
from scipy import integrate, optimize

# 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):
""" Calculation Method """
partial_func = partial(func, *args, **kwargs)
update_wrapper(partial_func, func)
return partial_func

#@jit
def sinputt(t, amplitude=5.0, frequency=0.2):
""" Calculation Method """
return amplitude * (1 + np.sin(frequency * t)) + 1

#@jit
def propensity_1(state, t, k):
""" Calculation Method """
return k * sinputt(t=t)

#@jit
def propensity_2(state, t, k):
""" Calculation Method """
return k * state

def functionieren(t, state, func):
""" Returns the function to be integrated """
return func(t=t, state=state)

def integrieren(t, dt, state, func):
""" Returns the integral """
return np.abs(
func=functionieren,
a=t,
b=dt,
args=(state, func),
epsabs=1.e-4,
epsrel=1.e-4,
limit=50)[0]
)

def optimizieren(dt, t, state, func, P_T):
""" Function for which root is to be found """
return P_T - integrieren(t, dt, state, func)

def solve(t, dt, state, func, P_T):
""" Function that finds root, thus t+dt """
return optimize.newton(
func=optimizieren,
x0=dt,
args=(t, state, func, P_T)
)

def run_step(P, T, state, t, propensities, funcs):
""" One step in the algorithm """
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

def run_sim(P, T, state, t, propensities, funcs):
""" The algorithm, final time=100 """
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)

if __name__ == '__main__':

# 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=t0+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
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.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