Setting stims and time after SaveState.restore()

Anything that doesn't fit elsewhere.
Post Reply
wwlytton
Posts: 64
Joined: Wed May 18, 2005 10:37 pm

Setting stims and time after SaveState.restore()

Post by wwlytton » Wed Aug 13, 2014 2:16 pm

After a SaveState.restore() , modifications are difficult and seem to depend on doing things at the correct time in the finitialize() sequence. I don't yet fully understand the reasons at each step but will post further as I better understand it.

Here is the fairly simple sim I was working with:
mss = h.SaveState(); myfile = h.File()
lsy, lnc, lns, lfih = [], [], [], [] # list of synapse, netcon, netstim, finithandlers
dend = h.Section(); dend.insert('pas'); dend.insert('hh')
vvec= h.Vector(); vvec.record(dend(0.1)._ref_v) # record voltage
tvec=h.Vector(); tvec.record(h._ref_t) # record time
syn = h.Exp2Syn(0.5,sec=dend)
syn.tau1 = 0.05;syn.tau2 = 5.3;syn.e = 0 # AMPA
ns = h.NetStim();ns.number = 1;ns.start = 750;
nc = h.NetCon(ns,syn); nc.weight[0] = 0.1; nc.active(1)
lsy.append(syn); lns.append(ns); lnc.append(nc)

Note that when saving/restoring the model structure must be identical -- hence the netstim should be in place prior to the save. The save was done as a cvode.event() thru an finitializehandler as follows.

def sssave (): print "Saving at t=%g"%(h.t); mss.save()
def setsave (): h.cvode.event(509,sssave)
fihl=h.FInitializeHandler(0,setsave)

after writing to a file and restoring the following handlers were set up to manipulate the time and netstim before running a new simulation using fihset() (see at bottom)
fihset((3,begsets),(1,ssrestore),(2,finalsets))

def begsets (): # fih3 at very beginning
lns[0].start=600 # hidden -- will no longer see this when print out at fih3
def ssrestore (): # fih1
mss.restore(1) # restore takes place here
lns[0].number=3; lns[0].interval=50
def finalsets (): # fih2 at very end
h.cvode.re_init() # needed to get recording to continue from old t

It seems odd that the NetStim had to be modified at 2 different locations.
The use of re_init() to get t correct also was mysterious and was something I attempted to examine further.
Without this, when recording under cvode (tvec.record(h._ref_t) # record time), the time in tvec[0] was 0
In order to fix this it appears necessary to do cvode.re_init() at the very end of finitialize()

This was effected by using flag 2 to an FInitializeHandler() -- this takes place then just before the return
doing this at flag 1 or using frecord_init() does not seem to work

It is not clear to me why this would be the case since it appears that re_init is already being called after handler1 in the finitialize() sequence. The 't' is preserved at the end of finitialize() regardless of whether this additional call is made.
nrncvode_set_t(t);
,,,,
nrn_fihexec(1); /* after INITIAL blocks, before fcurrent*/
,,,,
cvode_finitialize(); # calls re_init (/usr/site/nrniv/nrn/src/nrncvode/cvodestb.cpp:117:577)
nrn_record_init();
,,,,
nrn_fihexec(2); /* just before return */

def fihset (*args):
'set using tuple of number and routine'
global fihl
try: fihl
except: fihl=[0,0,0,0]
for i,rtn in args:
fihl=h.FInitializeHandler(i,rtn)

wwlytton
Posts: 64
Joined: Wed May 18, 2005 10:37 pm

Re: Setting stims and time after SaveState.restore()

Post by wwlytton » Wed Aug 13, 2014 4:47 pm

Here's full runnable code demonstrating the points made above
run python and then
execfile('file.py') # file.py is where you saved it

Code: Select all

from neuron import h,gui
import sys,os,time
h.load_file("stdrun.hoc")
from matplotlib import pyplot as mp
import numpy as np
mp.ion()
mss = h.SaveState(); myfile = h.File()
lsy, lnc, lns, lfih = [], [], [], [] # list of synapse, netcon, netstim, finithandlers
h.cvode.active(1)
h.tstop = 510

# setup sim: mkdend() insertsyn()
def mkdend ():
  global dend
  dend = h.Section(); dend.insert('pas'); dend.insert('hh')

def rec ():
  global vvec,tvec
  vvec= h.Vector(); vvec.record(dend(0.1)._ref_v) # record voltage
  tvec=h.Vector(); tvec.record(h._ref_t) # record time

## create a synapse with NetStim,NetCon
def insertSyn ():
  syn = h.Exp2Syn(0.5,sec=dend)
  syn.tau1 = 0.05;syn.tau2 = 5.3;syn.e = 0 # AMPA
  ns = h.NetStim();ns.number = 1;ns.start = 750;
  nc = h.NetCon(ns,syn); nc.weight[0] = 0.1; nc.active(1)
  lsy.append(syn); lns.append(ns); lnc.append(nc)


# savestate() and restorestate()
def fihset (*args):
  'set using tuple of number and routine'
  global fihl
  try: fihl
  except: fihl=[0,0,0,0]
  for i,rtn in args:
    fihl[i]=h.FInitializeHandler(i,rtn)

def fihclr (*args):
  global fihl
  fihl[:]=[0,0,0,0]
  print 'Confirming no FInitializeHandlers: {0:d}'.format(int(h.List('FInitializeHandler').count()))

def savestate (filestr):
  mss.save()
  myfile.wopen(filestr)
  mss.fwrite(myfile)
  myfile.close()
  
def loadstate (filestr):
  try: 
    myfile.ropen(filestr)
    mss.fread(myfile)
    myfile.close()
  except: print sys.exc_info()[0],' :',sys.exc_info()[1] 

def sssave (): print "Saving at t=%g"%(h.t); mss.save()
def setsave (): h.cvode.event(509,sssave)

def begsets (): # fih3 at very beginning
  prstiminfo('begin init')
  lns[0].start=600 # hidden -- will no longer see this when print out at fih3
  prstiminfo('end of fih3')

def ssrestore (): # fih1
  mss.restore(1)
  lns[0].number=3
  lns[0].interval=50
  st=prst()

def finalsets (): # fih2 at very end
  h.cvode.re_init()
  prstiminfo('final after fih2')

# full sequence
def saverun (f='aa.st'):
  mkdend(); rec(); insertSyn()
  fihset((0,setsave))
  h.run()
  savestate(f) # optionally save to a file


def restorerun (f='aa.st'):
  h.tstop=1e3
  # mkdend(); rec(); insertSyn() # don't need if resuming without restarting
  # loadstate(f)
  fihset((3,begsets),(1,ssrestore),(2,finalsets))
  h.run()
  
# plot and misc
def plot ():
  global f1
  f1=mp.plot(tvec, vvec)
  mp.xlabel('Time (ms)'); mp.ylabel('Vm (mV)')
  mp.show()

def rmnetstims ():
  lsy[:], lnc[:], lns[:] = [], [], [] # list of synapse, netcon, netstim
  chknetstim()

def chknetstim (): print '%d NetCons; %d NetStims; %d Exp2Syns'%(h.List('NetCon').count(),h.List('NetStim').count(),h.List('Exp2Syn').count())

def prst (st1=False):
  st0=[x for x in [h.t,dend.v,dend.m_hh,dend.h_hh,dend.n_hh]]
  print ['%.3f'%(x) for x in st0]
  if st1: print "Diffs:",[x-y for x,y in zip(st0,st1)]
  return st0

def prstiminfo (str=""):
    ns=lns[0]; 
    print str,':t={} ns.start={} ns.number={} ns.interval={}'.format(h.t,ns.start,ns.number,ns.interval)

# run sequence
saverun()
fihclr()
restorerun()
plot()

hines
Site Admin
Posts: 1510
Joined: Wed May 18, 2005 3:32 pm

Re: Setting stims and time after SaveState.restore()

Post by hines » Sun Aug 17, 2014 10:26 am

It seems odd that the NetStim had to be modified at 2 different locations.
At SaveState.save at t=509 The information that was saved included t=509, NetStim.start=750, number=1, and interval=10 and also the fact that there was an event on the queue to be delivered at 750.
The last is irrelevant since you are using SaveState.restore(1) which does not restore the queue, however it does restore the start, number, and interval. So the trick you were able to accomplish
was to set the Netstim.start to 600 from a handler just before the calls to INITIAL which ends up putting 600 on the queue. Then after the restore (which was after the calls to INITIAL) you could
reset the values of number and interval which were used when the event was delivered to the NET_RECEIVE block at 600.
In diagnosing what is going on, it is helpful to print the event queue (note that printing the queue from a type 3 FInitializeHandler will show the information from the end of the previous run), using
h.cvode.print_event_queue()
The use of re_init() to get t correct also was mysterious...
I added some print statements and see
ssrestore leave t=555.6 pc.t(0)=555.6
finalsets enter t=555.6 pc.t(0)=0
finalsets leave t=555.6 pc.t(0)=555.6
So the fact that the statement in the nrn_finitialize() implmentation
cvode_finitialize();
set the thread time to 0 but not the hoc time is the source of the puzzle. (cvode.re_init() uses the hoc t to initialize and also set the thread times). Normally
finitialize() sets the hoc t to 0 and SaveState.restore is called after the return from finitialize(). I failed to imagine that the restore would be done from an FInitializeHandler so that
hoc t would be changed after the c statement
t = 0.;
and before the call to cvode_finitialize().
Your work around is ok, but I should modify cvode_finitialize() to use the hoc t to initialize cvode

wwlytton
Posts: 64
Joined: Wed May 18, 2005 10:37 pm

Re: Setting stims and time after SaveState.restore()

Post by wwlytton » Tue Aug 19, 2014 10:36 am

if i understand correctly, there's a single thread here and, as in the case of multiple threads, the individual thread values have to be explicitly set to correspond to the time given by the interpreters h.t

It appears that in order to see these threads that one needs to set up an explicit ParallelContext
pc=h.ParallelContext()

then I was able to print out and compare the 2 different t's

print str,':h.t={} pc.t={} ns.start={} ns.number={} ns.interval={}'.format(h.t,pc.t(0),ns.start,ns.number,ns.interval)

hines
Site Admin
Posts: 1510
Joined: Wed May 18, 2005 3:32 pm

Re: Setting stims and time after SaveState.restore()

Post by hines » Tue Aug 19, 2014 11:10 am

Hoc h.t and all the thread times should generally agree when the interpreter has control, especially on return from finitialize() and fadvance(). Your's was a case that I consider a bug and
the fix is in http://www.neuron.yale.edu/hg/neuron/nr ... a0859594a8
The circumstances in which thread times may differ from each other and from h.t is while hoc/python is executing a callback function during simulation during use of the local variable time step and/or pc.psolve(...). In the latter each thread may integrate for a 'minimum netcon spike delay interval' before synchronizing.

Post Reply