for (x)

[HPCGuy]

Ted mentioned that the code I was working with had an error:

It used the code for (x): --

1) There is a lot of code that uses this sytnax
2) How can I make this code correct? What would I use in its place?

[ted]

There is nothing wrong with for (x) itself. It does what it is supposed to do. The problem is that the ability to write a program gives one the ability to write code that does something other than what one expects. Consequently, it is always necessary to verify that code works as intended. Also, it is often necessary to reread documentation because all implications of even the simplest statements may not be immediately apparent.

For example, consider this excerpt from the Programmer's Reference entry about the keyword for:

for (var) stmt

Loops over all segments of the currently accessed section. var begins at 0 and ends at 1. In between var is set to the center position of each segment. Ie. stmt is executed nseg+2 times.

Let's try this out.
oc>create axon
oc>access axon
oc>nseg = 1
oc>diam = 2

Code: Select all

oc>for (x) print x, diam(x)
0 2 
0.5 2 
1 2 

Which is what we expected, right?

Flush with understanding, we plow ahead with this
oc>for (x) diam(x) = x
and expect what? (be honest, now . . . ) But what really happens?

Code: Select all

oc>for (x) print x, diam(x)
0 1 
0.5 1 
1 1 

Which is what we should have expected, given what the Programmer's Reference says.

======== Start of rhetorical aside ========

A hand in the back raises and a voice asks "I don't get it."

OK, let's pretend to be the computer and execute these for loops ourselves, starting with the first pass through the first for loop (the loop that assigns values to diam(x)) and continuing to the last pass through the second for loop (the loop that reports the values of diam(x)).

The first for loop, which assigns values to diam:
1. On the first pass, x is 0, so diam(0) is set to 0. But the node at 0 really belongs to axon's first segment, so diam(0.5) is set to 0.
2. On the second pass, x is 0.5, so diam(0.5) is set to 0.5.
3. On the third and final pass, x is 1, so diam(1) is set to 1. But the node at 1 really belongs to axon's last segment, so diam(0.5) is set to 1.

The second for loop, which reports the values of diam:
1. On the first pass, reports diam(0), but the node at 0 really belongs to axon's first segment, so diam(0) is 1 i.e. the diameter at the middle of the first segment.
2. On the second pass, reports diam(0.5), which is 1.
3. On the third pass, reports diam(1), but the node at 1 really belongs to axon's last segment, so diam(0) is 1 i.e. the diameter at the middle of the last segment.

======== End of rhetorical aside ========

What would have happened if nseg were 3? The first two segments would have the "correct" (expected) diameters of 1/6 and 1/2, but the last segment's diameter would be 1 instead of the expected 5/6. To generalize, if f(x) is a function of x and foo is a range variable,
for (x) foo(x) = f(x)
will assign a value to foo in the last segment of the currently accessed section that is too large or too small, depending on whether f() increases or decreases over the interval 1-1/(2*nseg) . . . 1

The fix is simple. Don't use for (x) to assign values to range variables. Instead, use for (x,0). Quoting from the Programmer's Reference again,

for (var, expr) stmt

If the expression evaluates to a non-zero value, it is exactly equivalent to for (var) stmt If it evaluates to 0 (within float_epsilon ) then the iteration does not include the 0 or 1 points. Thus for(x, 0) { print x } is exactly equivalent to for (x) if (x > 0 && x < 1) { print x }

Let's verify that this works.
oc>for (x,0) diam(x) = x
Now what do we expect? (again, be honest . . . ). Here's what we get:

Code: Select all

oc>for (x) print x, diam(x)
0 0.5 
0.5 0.5 
1 0.5 

======== Start of rhetorical aside ========

A voice in the back mutters "This is terrible. I'm avoiding for (x) like the plague."

Not so fast, friend. for (x) is perfectly OK for accessing the values of range variables. It's particularly useful for membrane potential, because NEURON really does calculate the value of v at the 0 and 1 nodes. If the 0 or 1 end of a section is "floating" i.e. not attached to any other section or signal source (e.g. synapse, current clamp, voltage clamp), the potential at that point will be the same as the potential at the adjacent internal node. But if it is attached to a signal source or to the 0 or 1 end of another section, the potential at the point of attachment is calculated algebraically from the currents delivered by the signal sources, the membrane potentials at the adjacent internal nodes, and the values of the axial resistances between those nodes and the point of attachment (think "voltage divider").

So to wrap it all up, for (x) is fine for accessing the values of range variables, but for (x,0) is what should be used when assigning values to range variables.

======== End of rhetorical aside ========

And yes, there are many published models that use for (x) to assign values to range variables, with the result that certain parameters are too large or too small in the last segment of one or more sections. Which only goes to show that the modelers did not bother to verify the correctness of their model specifications. Are these errors fatal? Probably not in most cases. Remember that for many, if not most models, there is a fairly substantial volume of parameter space that will generate simulations that reproduce the particular phenomenon of interest or satisfy the criterion that concerned the modelers. But that doesn't mean anyone who is developing new code should deliberately embed errors in it, or neglect the responsibility of verifying that their model setup code actually does what they think it does. Go ye now and sin no more. Or someone else might put it, "Don't be evil, or at least don't be any more evil than is absolutely necessary to win."