Incremental simulation in direct memory mode

[iraikov]

Hello,
In standard NEURON, I use a checkpointing strategy where the overall simulation time is divided into equal segments, psolve is run sequentially for each segment, and the spike times and state traces are saved at the end of the psolve run for that segment. My understanding is that in standard NEURON, h.t is updated to the stop time for the segment, and when psolve is run again, the simulation starts from that time. However, when using CoreNeuron, it appears that the simulation always starts at t=0, even when h.t is set to a positive number. I'd like to ask whether the outlined behavior of psolve in standard NEURON is intended by design, and whether it would be possible to have consistent psolve start time semantics in both standard NEURON and CoreNEURON? Many thanks.

[pramodk]

I'd like to ask whether the outlined behavior of psolve in standard NEURON is intended by design, and whether it would be possible to have consistent psolve start time semantics in both standard NEURON and CoreNEURON?

Without looking into all details, I think it should be simple to simple to fix this behaviour.

I would like to understand bit more about the use case though - do you change any parameters or any model properties when you stop for each "segment" of the simulation?

[iraikov]

Thanks @pramodk In this specific case, no parameters are changed. The only purpose is to save the simulation results up to the current tstop.

[nrnhines]

That is a current limitation of the implementation which we intend to remove soon.

[nrnhines]

@iraikov It has been a long road but I think we finally have this resolved with PR neuronsimulator/nrn#1192 . There may be a bit of cleaning up to do, and likely a lot more testing, before it is merged into the master. But if it is not too much trouble for you, it would be nice to get your experience of it as a drop in substitute for your current "checkpoint" strategy with NEURON. With

from neuron import coreneuron
coreneuron.enable=True

Then pc.psolve(tstop) should execute on CoreNEURON and give exactly the same results (except for round off error differences due to high performance optimizations). I guess if you are using a GPU it might be better to hold off til all our CI
passes on that architecture. I've only tested on x86_64. Also, I don't believe you are using extracellular, LinearMechanism, CVode, HocEvents (that use HOC or Python during pc.psolve), or RxD, but , if you are, then I have to admit that CoreNEURON will not work for you. Let me know if you have any questions.

[iraikov]

Hi @nrnhines thank you so much for your hard work! I will test it out as soon as possible. I do have a question related to the calculation of LFP: you have advised in the past to use use_fast_imem to calculate i_membrane with CoreNEURON, but if sampling many cells, there might be too much memory overhead to wait until the end of the simulation to compute an LFP signal from the individual i_membrane recordings. My question whether now it might be possible to have a strategy where CoreNEURON is run for a short time (e.g. 1 ms), then LFP is computed for that time period, and then the simulation is resumed. Or would it be too much overhead to keep transferring data from CoreNEURON so frequently? Thank you!

[nrnhines]

Hi @iraikov . I think the overhead is too much for that strategy. In the future, the plan is to have NEURON and CoreNEURON share the same memory representation (CoreNeuron's) of the model but that is a long way off. Sadly, I've gotten into the habit of using Python callbacks at each time step to do the lfp[i] = M[i][j]*imem[j] per rank calculation where M has only one or a few rows and then reduce the ranks at tstop ( PtrVector is used to gather the imem on a rank into a Vector). For CoreNEURON, the best way I can think of at present is to construct a mod file (one instance per thread per lfp electrode) which has a constant Vector of scale factors computed by NEURON (A row of the above M matrix) and looks into the NrnThread to access the i_membrane_ array. Then at each time step the i_membrane_ value is multiplied by its scale factor and added to the mechanism instance's accumulation RANGE variable for that component of the (single) lfp electrode.
With those variables being recorded in a Vector, they are available to NEURON when CoreNEURON returns. That way, you can do long interval integration on CoreNEURON. I'm pretty confident that will work for cpu and if it works properly on NEURON then it should also work on CoreNEURON. I'm not sure what happens to performance on GPU.

[nrnhines]

@pramodk I've been thinking further about a mod file for calculating lfp that would work in CoreNEURON. It would seem to be amenable to the shared structure approach just added to patstim.mod. In the lfp case, the shared structure is just an array of doubles that contain the constant position factors computed by NEURON. The array should be in one-to-one corespondence with the i_mem array so that the lfp contribution is computed as

lfp = 0.0;
for (i = 0; i < nt.end; ++i) {
  lfp += factor[i]*nt.nrn_fast_imem->nrn_sav_rhs[i]
}

Of course if the nodes are permuted, then I guess the factor array should be copied and permuted instead of shared.
Anyway, one question to you is whether this would be better accomplished with a special gpu version that could do the sum of products in log time in parallel. Another speculative question, since this has now arisen twice in a short time span, is whether it is justified to introduce a new nmodl keyword analogous to BBCOREPOINTER that might be named NRNCORE_SHARED_DATA.
With that, the nmodl translator would emit something similar to the recently added VERBATIM blocks in the NEURON and CoreNEURON pattern.mod files for Info** pattern_stim_info_ref(...). I hesitate going that far however since it all may be moot given the latest ideas about closer integration of NEURON and CoreNEURON where the only low level memory is that of the layout currently in CoreNEURON.

[pramodk]

@pramodk I've been thinking further about a mod file for calculating lfp that would work in CoreNEURON

@nrnhines : could you take a look at what is done as part of #492? (see the doc). It was added for BBP use case with neurodamus+coreneuron interaction but may be the same/similar could be adapted for generic implementation?

[nrnhines]

Yes!
It's great not to have to do this again and not so well :)
After a quick skimming, only one detail stands out for me. The test is tests/unit/lfp/lfp.cpp so it is not clear at the moment what a python test would look like in nrn/tests/coreneuron using NEURON to setup the model, request lfp calculation, and verify that the various direct interaction modes all produce the same result.