Precision of t in fadvance

[neuromau]

Hello, I'm using a sinusoidally varying netstim in my model and so I added code to fadvance() to update the interval of the netstim. I only want to update the interval every so often and so I added code to check the value of t before resetting the interval. However, the check doesn't work because calculations involving t seem to be less precise than I expected.

For example below, I checked t to print its value every 100 ms. The same code outside of fadvance for a loop variable (r) produces the expected result whether I use the mod operator (%) or explicitly calculate it myself. But in fadvance(), it's not as precise and so I run into trouble even if I increase float_epsilon.

This code:

Code: Select all

{load_file("nrngui.hoc")}				// Standard definitions - NEURON library file

tstop = 1000
float_epsilon = 1e-5

print "Before run:"
for r=0, 999 {
	if ((r - int(r/100)*100)==0) {print "calc  r = ", r, "          r%100 = ", r%100}

	if (r%100 == 0) {print "%     r = ", r}
}
print " "
print "During run:"

proc advance() {
	if ((t - int(t/100)*100)==0) {print "calc  t = ", t, "          t%100 = ", t%100}

	if (t%100 == 0) {print "%     t = ", t}
	fadvance()
}

run()
quit()

produces this output:

Code: Select all

Before run:
calc  r = 0           r%100 = 0 
%     r = 0 
calc  r = 100           r%100 = 0 
%     r = 100 
calc  r = 200           r%100 = 0 
%     r = 200 
calc  r = 300           r%100 = 0 
%     r = 300 
calc  r = 400           r%100 = 0 
%     r = 400 
calc  r = 500           r%100 = 0 
%     r = 500 
calc  r = 600           r%100 = 0 
%     r = 600 
calc  r = 700           r%100 = 0 
%     r = 700 
calc  r = 800           r%100 = 0 
%     r = 800 
calc  r = 900           r%100 = 0 
%     r = 900 
 
During run:
calc  t = 0           t%100 = 0 
%     t = 0 
calc  t = 100           t%100 = 1.4168222e-11 
%     t = 100 
calc  t = 200           t%100 = 100 
calc  t = 300           t%100 = 100 
calc  t = 400           t%100 = 100 
calc  t = 500           t%100 = 100 
calc  t = 600           t%100 = 100 
calc  t = 700           t%100 = 3.5515768e-10 
%     t = 700 
calc  t = 800           t%100 = 7.1895556e-10 
%     t = 800 
calc  t = 900           t%100 = 1.0827534e-09 
%     t = 900 

And if I don't change float_epsilon, during fadvance() I only get:

Code: Select all

During run:
calc  t = 0           t%100 = 0 
%     t = 0 
calc  t = 600           t%100 = 100 

I don't understand why the calculation is less reliable within fadvance()?

[ted]

Roundoff error with finite precision arithmetic can do a lot of strange-looking things. Before delving into that, I'd first like to address the issue that led you down this rocky path. Do you want to generate a train of spike events whose ISIs are directly or inversely proportional to sin(2*PI*f*t)? Were you thinking of doing this by altering a NetStim's interval parameter in the course of a simulation, e.g. at some user-specified times t0, t1 . . . ti ?

[neuromau]

Do you want to generate a train of spike events whose ISIs are directly or inversely proportional to sin(2*PI*f*t)? Were you thinking of doing this by altering a NetStim's interval parameter in the course of a simulation, e.g. at some user-specified times t0, t1 . . . ti ?

Yes and Yes.

[neuromau]

Right now I have something to the effect of:

Code: Select all

period = 125
baserate = 10 // maximum firing rate (ms, interval)
amp = 45 // half of the range between min and max firing rate (ms, interval)
tstop = 1000

freq = 1000/period // Hz

numinrns = 3
double phasepref[numinrns]

phasepref[0] = .1
phasepref[1] = .3
phasepref[2] = .8

objref sinvecs[numinrns] // numinrns gives number of distinct input types (each type has a characteristic sin wave pattern to the intervals)
for r = 0, numinrns-1 {
	sinvecs[r] = new Vector(tstop/dt)
	sinvecs[r].sin(freq,2*phasepref[r],dt) // phasepref is an array of the preferred oscillation phase for each distinct input type
	sinvecs[r].mul(amp)
	sinvecs[r].add(amp+baserate)
}

... // create netstims and append to cells list

objref myrate
proc advance() {
	if (t%1 == 0) { // only update the intervals every 1 ms
		for r = 1, cells.count()-1 { // first cell in cells list is not a netstim
			myrate = cells.o(r) // myrate.pretype gives the index into the array of sin wave vectors for distinct input types
			myrate.pp.interval = sinvecs[myrate.pretype].x[t/dt]
		}
	}
	fadvance()
}

run()

[ted]

I think the best way to do this is with NMODL code. Two more questions before I have a go at it:
First, just to make sure I don't misunderstand, you want something that generates spike events at deterministic intervals (at least for starts), right?
Second, is the parameter "period" in ms or seconds?

[neuromau]

generates spike events at deterministic intervals (at least for starts), right?

Yes

parameter "period" in ms or seconds?

It's in ms.

Thanks for your help!