Multiple invocations of psolve when using direct memory interface

[iraikov]

Hello,
Thanks again for your work on direct memory integration with NEURON. One of my use cases involves repeated invocations of psolve as part of an optimization loop. This works in standard NEURON, but when attempting to do the same in CoreNEURON via the direct memory interface with NEURON, either psolve returns immediately and does not conduct the simulation, or causes a segmentation fault. For example, adding a second call to psolve in the test_spikes.py script causes a segmentation fault. Should it be expected the repeated calls to psolve in the same context is supported in CoreNEURON? And if so, is there a specific initialization that has to be performed in between calls to psolve? My current code only invokes finitialize before the next psolve invocation. Thanks a lot for your help.

[pramodk]

@iraikov : we have here very simple test to run multiple invocations but it may not be capturing all model aspects and there might be a bug.

If it's possible to provide a test to reproduce the issue then it would help us to debug and fix it faster.

@nrnhines : Did you encounter similar issue recently?

[iraikov]

Hi @pramodk thanks, this test script causes a segmentation fault on my computer as well. Perhaps there is something peculiar to my setup (MPICH+gcc on Debian). Let me try on Frontera with intel to see if the same issue occurs.

[pramodk]

Interesting! I can tomorrow try it on Ubuntu/Debian system. Could you tell me your cmake arguments used for neuron / coreneuron build?

[iraikov]

Hi @pramodk the test script worked with Intel MPI on Frontera. Another difference was that I built without Interviews. The compilation options for MPICH+gcc were:

cmake ../nrn -DCMAKE_INSTALL_PREFIX=$HOME/bin/nrnpython3 -DNRN_ENABLE_CORENEURON=ON -DPYTHON_EXECUTABLE=which python3 -DNRN_ENABLE_MPI=ON -DCMAKE_CXX_FLAGS="-D MPICH_SKIP_MPICXX=1" -DNRN_ENABLE_RX3D=ON

[nrnhines]

Fixed a segmentation fault on ubuntu 18.4, gcc7.5.0 with "here very simple test". The #377 will be merged soon.
However there are a number of other issues in this domain that need some work. I.e. a sequence of psolve with or without intervening finitialize, and return of all data when psolve ends.

[pramodk]

@iraikov : could you update coreneuron to latest master and try your test again? If you are building via NEURON then:

cd external/coreneuron/
git checkout master && git pull

and usual CMake & make install command to build.

[iraikov]

Thanks @nrnhines @pramodk I will try out the latest coreneuron. On a related note, my specific use case for the "sequence of psolve" functionality is an optimization loop, where the parameters of some density mechanisms are modified before the next psolve is run. Is it correct to expect that the modified model configuration will be transferred to coreneuron when the next psolve is executed? Thank you very much for the rapid fixes.

@iraikov : could you update coreneuron to latest master and try your test again? If you are building via NEURON then:

cd external/coreneuron/
git checkout master && git pull

and usual CMake & make install command to build.

[nrnhines]

Is it correct to expect that the modified model configuration will be transferred to coreneuron when the next psolve

Yes. This is the case where you will have intervening calls to finitialize. (by default a psolve directed to CoreNEURON will in fact
do the equivalent of finitialize() without a v_init arg. ie. voltages are left as is but the INITIAL blocks of all mechanisms are called, the event queue is cleared, vector record and play restarts, etc.)

[iraikov]

Hi @nrnhines @pramodk while the second psolve does not cause a segmentation fault, something does not look right with the recorded spike vectors. I am attaching a modified direct transfer spike script, which attempts to run psolve twice and compare the CoreNeuron spike train after each run. It does not look like the CoreNeuron spike vector is resized to 0 when finitialize gets called. Let me know if this is actually expected behavior.

import numpy as np
from neuron import h, gui
import sys

h('''create soma''')
h.soma.L=5.6419
h.soma.diam=5.6419
h.soma.insert("hh")
ic = h.IClamp(h.soma(.5))
ic.delay = 50.
ic.dur = 100.0
ic.amp = 0.01
h.tstop = 250

# for testing external mod file
#h.soma.insert("CaDynamics_E2")

h.cvode.use_fast_imem(1)
h.cvode.cache_efficient(1)

pc = h.ParallelContext()
pc.set_gid2node(pc.id()+1, pc.id())
mync = h.NetCon(h.soma(0.5)._ref_v, None, sec=h.soma)
pc.cell(pc.id()+1, mync)


nrn_spike_t = h.Vector()
nrn_spike_gids = h.Vector()
pc.spike_record(-1, nrn_spike_t, nrn_spike_gids)

v = h.Vector()
v.record(h.soma(.5)._ref_v, sec = h.soma)
i_mem = h.Vector()
i_mem.record(h.soma(.5)._ref_i_membrane_, sec = h.soma)
tv = h.Vector()
tv.record(h._ref_t, sec=h.soma)
h.run()
vstd = v.cl()
tvstd = tv.cl()
i_memstd = i_mem.cl()

from neuron import coreneuron
coreneuron.enable = True
coreneuron.verbose = 0

corenrn_all_spike_t = h.Vector()
corenrn_all_spike_gids = h.Vector()
    
pc.spike_record( -1 if pc.id() == 0 else (pc.id()),
                     corenrn_all_spike_t,
                     corenrn_all_spike_gids )

h.stdinit()
pc.psolve(h.tstop)


print('nrn_spike_t: %s' % str(nrn_spike_t.to_python()))
print('nrn_spike_gids: %s' % str(nrn_spike_gids.to_python()))
print('corenrn_all_spike_t: %s' % str(corenrn_all_spike_t.to_python()))
print('corenrn_all_spike_gids: %s' % str(corenrn_all_spike_gids.to_python()))
sys.stdout.flush()

# check spikes match
assert(len(nrn_spike_t)) # check we've actually got spikes
assert(len(nrn_spike_t) == len(nrn_spike_gids)); # matching no. of gids
assert(np.all(np.isclose(nrn_spike_t.to_python(), corenrn_all_spike_t.to_python())))
assert(nrn_spike_gids.to_python() == corenrn_all_spike_gids.to_python())

# running second time for testing multiple executions
h.stdinit()
pc.psolve(h.tstop)

print('nrn_spike_t: %s' % str(nrn_spike_t.to_python()))
print('nrn_spike_gids: %s' % str(nrn_spike_gids.to_python()))
print('corenrn_all_spike_t: %s' % str(corenrn_all_spike_t.to_python()))
print('corenrn_all_spike_gids: %s' % str(corenrn_all_spike_gids.to_python()))
sys.stdout.flush()

# check spikes match
assert(len(nrn_spike_t)) # check we've actually got spikes
assert(len(nrn_spike_t) == len(nrn_spike_gids)); # matching no. of gids
assert(np.all(np.isclose(nrn_spike_t.to_python(), corenrn_all_spike_t.to_python())))
assert(nrn_spike_gids.to_python() == corenrn_all_spike_gids.to_python())

print("v min: %g v max: %f" % (np.min(v.to_python()), np.max(v.to_python())))
if not bool(tv.eq(tvstd)):
    print("Voltage times are different")
    sys.exit(-1)
if v.cl().sub(vstd).abs().max() >= 1e-10:
    print(v.cl().sub(vstd).abs())
    print("Voltage difference greater than or equal to 1e-10")
    sys.exit(-1)
if i_mem.cl().sub(i_memstd).abs().max() >= 1e-10:
    print("i_mem difference greater than or equal to 1e-10")
    sys.exit(-1)

print("Voltage times and i_membrane_ are same and difference is less than 1e-10")
h.quit()

[nrnhines]

Thanks Ivan. #380 fixes this and will likely be merged to master after review by @pramodk

[iraikov]

Thanks @nrnhines !

[pramodk]

@iraikov : neuron and coreneuron is now updated. From main neuron repository you can do git pull && git submodule update --init --recursive and it should update all submodules.

We have also updated README significantly and one can now do nrnivmodl -coreneuron mod_dir to compile mod files for neuron as well as coreneuron. If you see any issues, please let us know!

[iraikov]

Thanks @pramodk the multiple invocations test now works as expected.

[pramodk]

Perfect! Thanks for confirmation !