[Question] cvode slows down model significantly

[xanderladd]

Hello NEURON team,

I am currently working with an L5 TTPC cell model that I have translated from HOC into python. When I turn on h.cvode.active(1) the model I am using becomes extremely slow (~ 45 seconds run time v.s. 5 seconds previously). Thus the results of my timing experiments are:

run model took: 6.373714923858643, cvode status: False
run model took: 49.04107451438904, cvode status: True

My question is -- why does Cvode result in such a dramatic slowdown? I have tried cvode using the same model and constant stimuli (set using IClamp amp) and it is indeed faster in that case. Ultimately, I am just looking to find the fastest model for these stimuli.

I thought the issue could be that I am using a high value for the DT (.09) but even when I use the regular DT (.02) I see the exact same results.
Here are the code solver stats

Cvode instance 0x55555591a6d0 cADpyr232_L5_TTPC1_0fb1ca4724[0].soma[0] statistics : 7878 CVode states
   13717 advance_tn, 9999 interpolate, 10000 init (0 due to at_time)
   85610 function evaluations, 21894 mx=b solves, 11475 jacobian setups

10000 init states seem like an issue.

Attached Working Example

cvode_issue.tar.gz

I attached a tarred version of the files needed for this code example. This setup might require recompile but should run without any external libraries besides NEURON, NumPy, os, and matplotlib. I also attached responses for each model, notably the cvode model has 20k timesteps instead of the expected 10k.

no cvode visualization

cvode visualization

It is attached but also here is the code for main.py

import os
os.chdir("neuron_files/")
from neuron import h
os.chdir("../")
import time
import numpy as np
from nrnpython_L5TTPC import Cell
import matplotlib.pyplot as plt


stim_num = '75'
ntimestep = 10000
h.load_file('neuron_files/run_model_cori.hoc') # in attached tar file
curr_stim = np.genfromtxt(f'stim_{stim_num}.csv',delimiter=',') #  in attached tar file


def run_model(param_set,dt,cell):

    volts_list, times_list = [], []
    # fill out model parameters
    cell.update_params(param_set)
    h.dt = dt
    if abs(dt*h.nstep_steprun*h.steps_per_ms - 1)  != 0:
            h.steps_per_ms = 1/(dt * h.nstep_steprun)
    # load stim and record variables
    v = h.Vector(curr_stim)
    v.play(cell.ic._ref_amp, dt)
    recording_vars = {k: h.Vector().record(v) \
               for k,v in zip(['t', 'v', 'stim_i', 'amp', 'dt'], \
                              [h._ref_t, cell.soma(.5)._ref_v,  \
                               cell.ic._ref_i, cell.ic._ref_amp, h._ref_dt])}

    h.v_init = -70
    h.tstop = ntimestep*dt 
    h.finitialize()
    # time simulation
    start = time.time()
    h.run()
    end = time.time()
    print(f"run model took: {end-start}, cvode status: {h.cvode.active()}")
    
    volts_list.append(recording_vars['v'].to_python()), times_list.append(recording_vars['t'].to_python())
    [recording_vars[key].play_remove() for key in recording_vars] # clear out recording vecs
    return np.array(volts_list), np.array(times_list)     
    
def plot_result(volts, postfix):
    fig = plt.figure()
    plt.plot(volts[0])
    plt.savefig(f'volts_run_{postfix}.png')


if __name__ == "__main__":
    dt = 0.09090909090909091
    # shows that using a large DT is not the main cause of this problem
    # dt =.02
    params = [2.0414857014281588e-05, -64.86061597380811, 7.526059454595134, 14.425443361666193, 0.029383427198445156, 3.390833660223704, 7.646600465203884, 0.5117614278381101, 0.0007967292525331661, 0.2212067846779962, 0.007515842711241828, 3.2673802922616044, 0.022771076489153656, 0.001231151499661452]
    
    
    l5_ttpc_cell = Cell()
    h.cvode.active(0)
    v1,t1 = run_model(params,dt, l5_ttpc_cell)
    h.cvode.active(1)
    v2, t2 = run_model(params,dt, l5_ttpc_cell)
    # plot and save reults
    plot_result(v1, 1), plot_result(v2, 2)
    print(h.cvode.statistics())    
   

[xanderladd]

also I can move this issue to discussion if that is more appropriate.

Thank you in advance for all your help so far!!

[nrnhines]

Have not run the model yet but will do so. Until then, on a brief look at main.py, my first guess is that the performance problem is because of

v.play(cell.ic._ref_amp, dt)

This forces cvode to reinitialize itself after every dt interval and start from that time as though it is a different initial value problem (with a possibly different value for cell.ic.amp). I haven't looked yet at your cur_stim trajectory but from your graphs I expect it is just a long step that begins at 1000 and ends a little while after 8000. ie. no need to send an event every dt, an event is needed only when there is a discontinuity in cell.ic.amp.
CVode works best when events (discontinuities that force a reinitialization) are rare.

Ok. Looking back near the top of your post I see

13717 advance_tn, 9999 interpolate, 10000 init

So CVode reinitialized itself 10000 times and each time interpolating itself back to the exact i*dt time to take into account a new value for cell.ic.amp .

I expect that the abscissa on your second plot is not your recorded time.

I'll get back after I run the model and see if my speculations are confirmed. If so a simple modification to v.play should suffice to fix the performance problem.

[nrnhines]

With

$ git diff main.py
diff --git a/main.py b/main.py
index eea91f7..4fda438 100644
--- a/main.py
+++ b/main.py
@@ -1,3 +1,4 @@
+import discons
 import os
 os.chdir("neuron_files/")
 from neuron import h
@@ -23,8 +24,11 @@ def run_model(param_set,dt,cell):
     if abs(dt*h.nstep_steprun*h.steps_per_ms - 1)  != 0:
             h.steps_per_ms = 1/(dt * h.nstep_steprun)
     # load stim and record variables
-    v = h.Vector(curr_stim)
-    v.play(cell.ic._ref_amp, dt)
+    v, tv = discons.discon(curr_stim, dt)
+    print(v, tv)
+    v = h.Vector(v)
+    tv = h.Vector(tv)
+    v.play(cell.ic._ref_amp, tv)
     recording_vars = {k: h.Vector().record(v) \
                for k,v in zip(['t', 'v', 'stim_i', 'amp', 'dt'], \
                               [h._ref_t, cell.soma(.5)._ref_v, 

and

$ cat discons.py
def discon(vin, dt):
    vout = []
    tout = []
    vprev = -1e9
    for i, v in enumerate(vin):
        if v != vprev:
            if i > 0:
                tout.append(i * dt)
                vout.append(vprev)
            tout.append(i * dt)
            vout.append(v)
        vprev = v
    return vout, tout

I get

run model took: 10.246886491775513, cvode status: False
[0.0, 0.0, 0.25999999375692084, 0.25999999375692084, 0.0] [0.0, 89.18181818181819, 89.18181818181819, 750.4545454545455, 750.4545454545455]
run model took: 8.526818037033081, cvode status: True

Cvode instance 0x557b9d225850 cADpyr232_L5_TTPC1_0fb1ca4724[0].soma[0] statistics : 4946 CVode states
   2882 advance_tn, 2 interpolate, 3 init (0 due to at_time)
   6975 function evaluations, 4079 mx=b solves, 463 jacobian setups

But when you plot recording_vars['v'] don't forget to plot against recording_vars['t']

[xanderladd]

Confirming that specifying dt caused all of those init calls. As you specified, discons is now working to specify change points that actually require a call to init. So now I am also getting

run model took: 6.169771432876587, cvode status: False
run model took: 5.358117580413818, cvode status: True

Cvode instance 0x55555591a6d0 cADpyr232_L5_TTPC1_0fb1ca4724[0].soma[0] statistics : 4946 CVode states
   3011 advance_tn, 2 interpolate, 3 init (0 due to at_time)
   7297 function evaluations, 4272 mx=b solves, 504 jacobian setups

which shows the expected behavior of cvode. Thank you for your help it is greatly appreciated and I will close this issue.