Cannot get Vector.play() and NetCon/NetStim work with CVODE and SaveState.restore()

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

Post Reply
JustasB
Posts: 28
Joined: Tue Dec 08, 2015 8:17 pm
Location: Tempe, AZ, USA
Contact:

Cannot get Vector.play() and NetCon/NetStim work with CVODE and SaveState.restore()

Post by JustasB »

I am having trouble using Vector.play() and NetCon/NetStim together with CVODE and SaveState.restore(). I run the simulation until my model reaches steady state, then I save with SaveState.save(), and then I inject an arbitrary current into the cell after SaveState.restore(). A NetCon/NetStim combination is used to collect membrane potential values at some interval.

However, I cannot get both Vector.play() AND NetCon/NetStim to work. If I use SaveState.restore(), Vector.play() does not work, if I use SaveState.restore(1) the simulation crashes on h.continuerun() when the NetCon/NetStim combination is present. Is there a way to make this work?

Below is a minimal code to reproduce the problem:

Open python session and run:

Code: Select all

from neuron import h,gui

class Collector:
    def __init__(self, collection_period_ms, expr):
        from neuron import h
        self.collection_period_ms = collection_period_ms
        self.stim = h.NetStim(0.5)
        self.stim.start = 0
        self.stim.number = 1e9
        self.stim.noise = 0
        self.stim.interval = self.collection_period_ms
        self.con = h.NetCon(self.stim, None)
        self.con.record((self.collect, expr))
        self.fih = h.FInitializeHandler(self.clear)
        self.clear()
    def clear(self):
        self.values = []
    def collect(self, new_value):
        self.values.append(new_value())

class Tester:
    def build(self):
        self.soma = h.Section()
        self.soma.insert('pas')
        self.soma.insert('hh')
        self.ic = h.IClamp(0.5, sec=self.soma)
        self.ic.delay = 0
        self.ic.dur = 1e9
        self.ic.amp = 0
        self.col = Collector(1.0, lambda: h.t) # Including this line and using .restore(1), crashes continuerun()
    def save(self):
        self.ss = h.SaveState()
        self.ss.save()
        self.sf = h.File('state.bin')
        self.ss.fwrite(self.sf)
    def restore(self,arg):
        self.ns = h.SaveState()
        self.sf = h.File('state.bin')
        self.ns.fread(self.sf)
        self.ns.restore(arg)


t = Tester()
t.build()
h.stdinit()
h.tstop = 500
h.cvode_active(1)
h.nrncontrolmenu()
h.newPlotV()
h.run()
t.save()
Then in another python session, run:

Code: Select all

from neuron import h,gui

class Collector:
    def __init__(self, collection_period_ms, expr):
        from neuron import h
        self.collection_period_ms = collection_period_ms
        self.stim = h.NetStim(0.5)
        self.stim.start = 0
        self.stim.number = 1e9
        self.stim.noise = 0
        self.stim.interval = self.collection_period_ms
        self.con = h.NetCon(self.stim, None)
        self.con.record((self.collect, expr))
        self.fih = h.FInitializeHandler(self.clear)
        self.clear()
    def clear(self):
        self.values = []
    def collect(self, new_value):
        self.values.append(new_value())

class Tester:
    def build(self):
        self.soma = h.Section()
        self.soma.insert('pas')
        self.soma.insert('hh')
        self.ic = h.IClamp(0.5, sec=self.soma)
        self.ic.delay = 0
        self.ic.dur = 1e9
        self.ic.amp = 0
        self.col = Collector(1.0, lambda: h.t) # Including this line and using .restore(1), crashes continuerun()
    def save(self):
        self.ss = h.SaveState()
        self.ss.save()
        self.sf = h.File('state.bin')
        self.ss.fwrite(self.sf)
    def restore(self,arg):
        self.ns = h.SaveState()
        self.sf = h.File('state.bin')
        self.ns.fread(self.sf)
        self.ns.restore(arg)


t = Tester()
t.build()

# Create a current ramp
rv = h.Vector([0,   0,  10,   0])
tv = h.Vector([0, 500, 600, 600])
rv.play(t.ic._ref_amp, tv, 1)

# Restore steady state and run ramp
h.stdinit()
t.restore(0) # 0 does not play() vector - no ramp - as expected per: https://www.neuron.yale.edu/neuron/static/new_doc/simctrl/savstate.html#SaveState.restore
h.cvode_active(1)
h.nrncontrolmenu()
h.newPlotV()
h.continuerun(700)
t.col.values # Shows NetStim/Con collected values for 500-700ms - as expected


# Restore, but this time preserve the event queue
h.stdinit()
t.restore(1) # 1 causes continuerun() crash when netstim/Collector is present
h.cvode_active(1)
h.nrncontrolmenu()
h.newPlotV()
h.continuerun(700) # This line crashes NEURON

The crash shows this message:

Code: Select all

Assertion failed!

Program: [...]\python.exe
File: ../../../nrn/src/nrniv/../nrncvode/netcvode.cpp, Line 2233

Expression: tdiff == 0.0 || ( gcv_->tstop_begin_ <= tt && tt <= gcv_->tstop_end_)
What am I missing? Is it possible to make Vector.play() and NetStim/NetCon work with CVODE and SaveState?
JustasB
Posts: 28
Joined: Tue Dec 08, 2015 8:17 pm
Location: Tempe, AZ, USA
Contact:

Re: Cannot get Vector.play() and NetCon/NetStim work with CVODE and SaveState.restore()

Post by JustasB »

Strangely, writing a post here has made it easier to solve my own problem!

The problem appears to be with the NetStim.start time, before the restore, being set to 0. If I set the NetStim.start equal to the h.t value after the restore, the crash dissapears. So I believe I've found a workaround. It is:

Code: Select all

t = Tester()
t.build()

rv = h.Vector([0,   0,  10,   0])
tv = h.Vector([0, 500, 600, 600])
rv.play(t.ic._ref_amp, tv, 1)

t.restore(1)               # Will set the h.t to the time it was saved
t.col.stim.start = h.t # Ensures NetStim.start is not before h.t

h.stdinit()
t.restore(1)
h.cvode_active(1)
h.nrncontrolmenu()
h.newPlotV()
h.continuerun(700) # Runs without crashing with the expected ramp
However, I cannot tell if the requirement to make sure NetStim.start is >= h.t of the saved state by design or a bug. NetStim docs do not mention restrictions on .start and a somewhat related post appears to have been a bug that has been fixed.

Thoughts?
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Cannot get Vector.play() and NetCon/NetStim work with CVODE and SaveState.restore()

Post by ted »

One of the nice things about initialization to a saved state is that, after executing

savestate.restore()

you should be able to do

h.t = 0 # why allow t and all graphs to start at a completely accidental time?

Also, since your new simulation would start at t = 0, the question of what happens when NetStim.start is less than the h.t of the saved state would never have arisen.

Note that since savestate.restore() changes states,
1. the adaptive integrator must be notified that it must update its internal copy of the initial states
and
2. changing states after calling h.stdinit() (well, actually after calling h.finitialize(), which is one of the things that h.stdinit() does) can result in the storage of incorrect values in the first element of recorded Vectors

The Python equivalent of a common, effective hoc idiom for initialization that involves savestate.restore() is

Code: Select all

h.finitialize(h.v_init)
h.savestate.restore()
h.t = 0
if (h.cvode.active()):
  h.cvode.re_init() # make adaptive integrator update its copy of initial states
h.frecord_init()# re-initialize Vector recording
A much simpler way to capture a variable at fixed intervalls is to use the Vector class's
vdest.record(var_reference, Dt)
syntax--see the documentation at https://www.neuron.yale.edu/neuron/stat ... tor.record. Easier to read and maintain, and it employs internal compiled code which is much faster than interpreter-based data capture.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: Cannot get Vector.play() and NetCon/NetStim work with CVODE and SaveState.restore()

Post by hines »

I'll need to look into this a bit more but for now I just want to caution that if SaveState is chosen so that events on the queue are part of the
state, that one cannot in general set h.t = 0 because all the outstanding events will be delivered at a time after the save time. Also many
ArtificialCell (such as NetStim) store the previous event time and presume that event arrival involves monotonically increasing time.
For the case that events are not part of the saved state, I do believe it is allowed to set h.t = 0 as long as no existing object relies on a stored value of
t as some earlier time. e.g. no NetStim exists. But I am not certain of that and a test would
be in order. The relevant documentation fragments (which are ambiguous in this matter) are:

class SaveState
...
The outstanding event delivery times are absolute. When restored, all outstanding events will be cleared and the restored event times and NetCon info will take their place. Note that it is not in general possible to change the value of t in a network simulation since most NET_RECEIVE blocks keep t0 (the last event time) as part of their state

SaveState.restore
...
If the arg is 1, then the event queue is not cleared and no saved events are put back on the queue. Therefore any Vector.play and/or FInitializeHandler events on the queue after finitialize() are not disturbed.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: Cannot get Vector.play() and NetCon/NetStim work with CVODE and SaveState.restore()

Post by hines »

I've reproduced the results of your code. An alternative to your fix using t.restore(0), which restores the event queue to the time of save.
is to realize that the save phase did not have any knowledge of rv.play(t.ic._ref_amp, tv, 1) and in particular that on restore there needs to
be an event for the vector play which comes from the save phase. So during the save phase I use:

Code: Select all

  t = Tester()
  t.build()
  # Create a current ramp
  rv = h.Vector([0,   0]) 
  tv = h.Vector([0, 510])
  rv.play(t.ic._ref_amp, tv, 1)
  h.stdinit()
and for the restore phase I used

Code: Select all

  t = Tester()
  t.build()
  # Create a current ramp
  rv = h.Vector([0,   0,  20,   0])
  tv = h.Vector([0, 510, 600, 600])
  rv.play(t.ic._ref_amp, tv, 1)
  # Restore steady state and run ramp
  h.tstop = 700
  h.stdinit()
Note that I started the ramp at 510 instead of 500 to emphasize that there needs to be an event saved. The vectors themselves do not have to be
the same size or values as they are not part of the saved state. Ending the Vector at 500 for the save state at t=500, could only work if the
event was not already delivered. Of course, if the y, t points of the vector play are the same for save and restore, there should be no problem
Post Reply