Using event delivery system with fixed time step method.

Anything that doesn't fit elsewhere.
Post Reply
ngreiner
Posts: 9
Joined: Tue Feb 28, 2017 8:12 am

Using event delivery system with fixed time step method.

Post by ngreiner » Fri Jun 29, 2018 5:56 am

Dear Ted,
I am encountering an issue when using the event delivery system while cvode_active() == 0.
Basically, the following code

Code: Select all

from neuron import h
h.load_file('stdrun.hoc')

sec = h.Section()
sec.insert('hh')

def hi():
	print 'hi'
	
h.dt = 0.025
h.tstop = 5.0

h.finitialize()
h.cvode.event(2.5, hi)
h.run()

print h.t
produces the following output

Code: Select all

NEURON -- VERSION 7.5 master (c693a84) 2017-12-07
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2016
See http://neuron.yale.edu/neuron/credits

5.0
indicating that the event that I try to set up is not taken into account during the simulation.
Instead, if I use the variable time step method, by modifying the previous code as follows

Code: Select all

from neuron import h
h.load_file('stdrun.hoc')

sec = h.Section()
sec.insert('hh')

def hi():
	print 'hi'
	
h.dt = 0.025
h.tstop = 5.0

h.cvode_active(1)

h.finitialize()
h.cvode.event(2.5, hi)
h.cvode.solve(h.tstop)

print h.t
I get the following output

Code: Select all

NEURON -- VERSION 7.5 master (c693a84) 2017-12-07
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2016
See http://neuron.yale.edu/neuron/credits

hi
5.0
which is what I expected to obtain with the first piece of code as well.
Finally, if I keep cvode_active()==0 but still try to use h.cvode.solve(h.tstop) instead of h.run(), as in

Code: Select all

from neuron import h
h.load_file('stdrun.hoc')

sec = h.Section()
sec.insert('hh')

def hi():
	print 'hi'
	
h.dt = 0.025
h.tstop = 5.0

h.finitialize()
h.cvode.event(2.5, hi)
h.cvode.solve(h.tstop)

print h.t
The output which I get is

Code: Select all

NEURON -- VERSION 7.5 master (c693a84) 2017-12-07
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2016
See http://neuron.yale.edu/neuron/credits

hi
1000000000.0
indicating that the integrator went much passed the required limit of 5.0.
Would you be able to help me understand these outcomes ?

Thank you very much in advance,

Nathan

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

Re: Using event delivery system with fixed time step method.

Post by ted » Fri Jun 29, 2018 1:21 pm

Good question. Has nothing to do with fixed step or adaptive integration. It's all about initialization.

run() is defined in nrn/share/nrn/lib/hoc/stdrun.hoc

Code: Select all

proc run() {
        running_ = 1
        stdinit()
        continuerun(tstop)
}
The same file contains the definition of stdinit(). One of the things that stdinit() does is call init(), and init() in turn calls finitialize(). And finitialize() clears the event queue, discarding the event that had been planted by h.cvode.event(2.5, hi).
h.cvode.solve() advances the solution in time, without calling finitialize() (see the documentation of CVode.solve in the Programmer's Reference).

To force a particular action at a time that is known before a simulation is launched, it is generally best to use an FInitializeHandler because this ensures that the event is placed on the queue during the process of initialization (see the Programmer's Reference documentation of FInitializeHandler).

To fix your first example, change
h.cvode.event(2.5, hi)
to
# h.cvode.event(2.5, hi)
fih = h.FInitializeHandler(2.5, hi)

Then every time you call
h.run()
your hi() procedure will be called when t==2.5 (you might want to add a statement that reports the value of t to hi(), just to verify).

By the way, since run() (eventually) calls finitialize(), there's no need for your code to call h.finitialize() directly.

ngreiner
Posts: 9
Joined: Tue Feb 28, 2017 8:12 am

Re: Using event delivery system with fixed time step method.

Post by ngreiner » Fri Jun 29, 2018 2:13 pm

I see. Thank you for the quick reply.

Would you be kind enough to also explain to me why using h.cvode.solve(h.tstop) leads to h.t=10000000 (despite the fact that h.tstop==5.0) when cvode.active()==0 ?

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

Re: Using event delivery system with fixed time step method.

Post by ted » Thu Jul 05, 2018 12:12 pm

Would you be kind enough to also explain to me why using h.cvode.solve(h.tstop) leads to h.t=10000000 (despite the fact that h.tstop==5.0) when cvode.active()==0 ?
Sounds as reasonable as any other value for t would be. If the adaptive integrator isn't active, why would trying to exercise adaptive integration cause t to have any particular value at all? In a less imperfect world, an appropriate "now, now, mustn't do that" message would be generated. In a completely perfect world, h.cvode.solve() would never be called unless the adaptive integrator is active.

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

Re: Using event delivery system with fixed time step method.

Post by hines » Sat Jul 07, 2018 6:14 pm

cvode.active(1) is primarily for using cvode with fadvance() and ParallelContext.psolve(tstop). It is merely inadvertent that cvode.solve(tstop)
does its work independently of the value of cvode.active(). Note, though, that cvode.solve does not make sense without a precursor call
to cvode.re_init() which intializes the cvode integrator to the current state of the model (and that does not do its work without cvode.active() == 1).
So cvode.solve is not much practical use without cvode.active(1) as the latter is necessary for cvode.re_init() to work.

breakwave922
Posts: 79
Joined: Mon May 19, 2014 10:39 pm

Re: Using event delivery system with fixed time step method.

Post by breakwave922 » Sun Oct 14, 2018 6:15 pm

Dear Ted and Hines,

Since my question is similar to the thread here, I guess it's better to ask it here rather than having a new thread. I've been mired in the problem for several days...

I have a network model which I already implemented calculating LFP at each time step. But to simplify, I don't have to calculate it at every time step, but just need to calculate it at time steps with less resolution. For example, my network model simulates at dt=0.05, but for LFP, I only need to calculate it at 1ms step interval. Having checked several threads relating to LFP calculating in real-time, I found I can use cvode.event and FInitializeHandler to realize. It works well initially, calling calculation of LFP and record it at time points I need. But after I increased my simulation time, e.g. tstop=9,000ms, when I checked recorded time points for LFP, they were not accurately at 1ms step beyond some time point.

I created a simple version of code similar to the one in this thread, which can also reproduce my problems, where I called a proc time_rec() at every 1ms to record time progress.

Code: Select all

load_file("nrngui.hoc")
dt = 0.05
steps_per_ms= 20
tstop = 9000

dt_LFP=1
objref testvec
 testvec=new Vector()

proc time_rec(){
testvec.append(t)
//print "hi at: ",t
cvode.event(t + dt_LFP, "time_rec()")	
}

objref fih2
fih2 = new FInitializeHandler("time_rec()")

run()

strdef cmd
proc writeLFP() { localobj f
    sprint(cmd,"test_time.dat")
    f = new File(cmd)
    f.wopen()
    testvec.printf(f)
    f.close()
 }
 writeLFP()
However, when checking saved test_time.dat file with time points recorded, you will find that till 8,255ms time points were saved correctly, but beyond it, the intervals increased to 1.05ms, until the end of simulation length.

Code: Select all

8249
8250
8251
8252
8253
8254
8255
8256.05
8257.1
8258.15
8259.2
8260.25
...
8990
8991.05
8992.1
8993.15
8994.2
8995.25
8996.3
8997.35
8998.4
8999.45
Not sure what causes this.
Any help would be appreciated.

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

Re: Using event delivery system with fixed time step method.

Post by ted » Tue Oct 16, 2018 12:16 am

Looks like a roundoff error problem. You need a dt that satisfies two criteria:

1. A whole multiple of dt must equal your sampling interval. That's pretty obvious, and 0.05 does satisfy this criterion.

2. The binary representation of dt must have a finite number of bits, AND that representation must be short enough to be handled by double precision arithmetic without error. In retrospect this is also obvious; if dt does not satisfy this criterion, sooner or later roundoff error will rear its ugly head. 0.05 does NOT satisfy this criterion, because it equals the product 1/2 * 1/2 * 1/5. The binary representation of the product of the first two terms is just 0.01; so far so good. But the binary representation of 1/5 is 0.0011001100110011... so 0.05 cannot be represented exactly with a finite number of bits.

Try a dt that is a whole multiple of an integer power of 1/2, e.g. 0.0625 (1/16) or 0.03125 (1/32).
Does the problem go away?

breakwave922
Posts: 79
Joined: Mon May 19, 2014 10:39 pm

Re: Using event delivery system with fixed time step method.

Post by breakwave922 » Tue Oct 16, 2018 12:25 pm

ted wrote:
Tue Oct 16, 2018 12:16 am

2. The binary representation of dt must have a finite number of bits, AND that representation must be short enough to be handled by double precision arithmetic without error. In retrospect this is also obvious; if dt does not satisfy this criterion, sooner or later roundoff error will rear its ugly head. 0.05 does NOT satisfy this criterion, because it equals the product 1/2 * 1/2 * 1/5. The binary representation of the product of the first two terms is just 0.01; so far so good. But the binary representation of 1/5 is 0.0011001100110011... so 0.05 cannot be represented exactly with a finite number of bits.

Try a dt that is a whole multiple of an integer power of 1/2, e.g. 0.0625 (1/16) or 0.03125 (1/32).
Does the problem go away?
Ted,

It makes sense. After I change dt= 0.0625 (1/16) to meet the second criteria, the problem goes away.
Thanks a lot.

Post Reply