vector element differs for constant or variable index

[tedballou]

I also am having a problem with inconsistent data values, in my case in a hoc program designed to calculate motor unit force following a train of spikes. The program seems to run fine, but when I use a very long spike train, the force exhibits anomalous blips. While debugging this I found NEURON to be in a very strange state in which the value of a vector was inconsistent, depending on whether the index was specified by a constant or by a variable assigned the value of the constant:

oc>ii
        82669
oc>trainforce.x[82669]
        0
oc>trainforce.x[ii]
        9.0028106

The problem occurs on a WIndows 7 64-bit system running NEURON 7.2, and also on a WIndows XP pro system running NEURON 7.2; it does not occur on my macbook running OSX 10.5.8, on NEURON 7.0.

[ted]

Strange indeed. For each of your machines can you please tell me
1. its OS and whether it the 64 or 32 bit version of that OS
2. the detailed version number of NEURON (from the first line of NEURON's banner)
3. the value hoc returns when calculating
ii - 82669

[tedballou]

Hi Ted, Thanks for the quick reply. I think this is pretty spooky....

1. its OS and whether it the 64 or 32 bit version of that OS
2. the detailed version number of NEURON (from the first line of NEURON's banner)

Windows 7 64 bit. NEURON -- VERSION 7.2 (445:d8e244719b0d) 2010-06-10
Windows XP 32 bit. NEURON -- VERSION 7.2 (422:615063ede00a) 2010-02-26
Mac Leopard 32 bit. NEURON -- Release 7.0 (281:80827e3cd201) 2009-01-16

3. the value hoc returns when calculating
ii - 82669

oc>ii - 82669
-1.4551915e-11
oc>

[ted]

oc>ii - 82669
-1.4551915e-11
oc>

I see the error under 32 bit CentOS 5.5
Well, that's the clue--the cause is roundoff error. ii is just a tad too small to be treated as 82669. Instead it is being treated as if it were 82668, which is the last nonzero element in the Vector.

A workaround is to increase the value of float_epsilon slightly. See the documentation of float_epsilon
http://www.neuron.yale.edu/neuron/stati ... at_epsilon
Note that the default value of float_epsilon is 1e-11. Consequently ii falls a bit shy of what it needs to be in order to be treated as equal to 82669.

If I increase float_epsilon just a bit
float_epsilon = 2e-11
then trainforce.x[ii] returns 0

One might well ask why roundoff error is more of a problem on one machine with one OS, but not on another with a different OS. And why ii--obviously an index--is being computed in a way that gives it a nonzero fractional part. Perhaps here
for ii=t1.x[j]/dltaT, upperlimit {
or here
for ii=0,maxTidx-1 {
since maxTidx is computed as
maxTidx = (t1.x[t1.size-1] + TwtchDur)/dltaT

And incidentally what protects this
upperlimit = t1.x[j+1]/dltaT-1
from having a nonzero fractional part?

(all of these are in twitch.hoc)

Suggestion: wrap ratios inside int() (or int(ratio + 0.5)) as necessary to force true integer results
for ii=int(t1.x[j]/dltaT), upperlimit {
maxTidx = int((t1.x[t1.size-1] + TwtchDur)/dltaT)
upperlimit = int(t1.x[j+1]/dltaT-1)
and see if that doesn't fix the problem. I just tried it (including the + 0.5 part) and the bug vanished. Whether the resulting twitch waveform is correct, or whether I have only introduced some other bug (by adding 0.5 before applying int()) is up to you to judge.

[tedballou]

Thanks Ted,

Forcing the indices to be integral does the trick!

Cheers, Ted

[ted]

That's a relief . . . a successful WAG (well, mayby only half-WAG) means I don't have to try to actually understand how your program works.