needed documentation - behaivour of ShapePlot affecting cell geometry

[lucasg52]

Context

I was working on a project that involved preforming deep searches for propagation thresholds in cells with randomly generated geometries (via python).
As such, tiny changes in cell geometry could completely throw off the data.

Overview of the issue

After losing my marbles over some extremely inconsistent simulation behaviors, I realized that using ShapePlot (in order to view channel activity) was causing small changes in the diameters of axon sections near branch points, and was causing smaller changes everywhere else.
Eventually, I accepted that I would no longer the GUI to observe my simulations, since I couldn’t find any information on how to remedy this problem.

Obviously, I am now aware that ShapePlot is intended to adjust diameters near branch points via tapering (though it appears to cast float64 diameters down to float32 somewhere along the way--maybe this is unintended), but I think this should be documented, especially considering that this was the recommended method of viewing simulations when I attended the NEURON course.

I am also aware that pt3d data must be generated for this to be possible, but there's no reason this should destroy previously defined diameter information. One simply does not expect their model to literally change when they try to view it in the GUI.

NEURON setup

• Version: 8.2.2 HEAD (93d41fa) 2022-12-15
• Installation method: pip install
• OS + Version: Linux Mint 21.2 x86_64
• Compiler + Version: probably gcc

Minimal working example - MWE

from neuron import h

# We begin by creating two sections with unequal, fractional diameters:

sec_a, sec_b = h.Section("a"), h.Section("b")
sec_a.diam = 0.6
sec_b.diam = 0.2
sec_a.nseg = sec_b.nseg = 9

# connect them to make a 'T'
sec_b.connect(sec_a(0.5))
h.topology()

# To see the effects of PlotShape, we can easily scrape diameters using an h.RangeVarPlot, along
# with h.Vector-s:
rvp = h.RangeVarPlot("diam", sec_a(0), sec_a(1))
x, diams1, diams2, diams3 = [h.Vector() for i in range(4)]
# save the initial diameters, then compare to new diameters:
rvp.to_vector(diams1,x)

h.define_shape()    # define_shape, despite the fact that it generates pt3d data, does not affect
                    # diameter.

rvp.to_vector(diams2,x)
diff = list(h.Vector().copy(diams1).sub(diams2))
print(f"\n{diff=}")

h.PlotShape(False)  # PlotShape does affect diameter.

rvp.to_vector(diams3,x)
diff2 = list(h.Vector().copy(diams3).sub(diams2))
print(f"\n{diff2=}")

[nrnhines]

Thank you for the issue. Chatgpt more or less misses the point about why this happens when asked "With the NEURON simulator, ShapePlot seems to change diameters. Why is that?" But gets it, with a subsequent, "it looks like some of the differences are due to 3d points are stored in 32 bit floats"

The diameter change with your above example is entirely due to define_shape() storing the the 3d point information
it calculates from diam,L double precision data into 32bit single precision floating point structure of x,y,z,d arrays. (PlotShape calls define_shape internally). Many functions, such as segment.area() and finitialize(), check the diam_changed flag and update diam based on the 3d point data. Thus, when I modified your code to insert a fragment:

$ git diff
diff --git a/test.py b/test.py
index 0073773..c731983 100644
--- a/test.py
+++ b/test.py
@@ -25,6 +25,13 @@ rvp.to_vector(diams2,x)
 diff = list(h.Vector().copy(diams1).sub(diams2))
 print(f"\n{diff=}")
 
+sec_a(.5).area()
+diams2a = h.Vector()
+rvp.to_vector(diams2a,x)
+diff2a = list(h.Vector().copy(diams2a).sub(diams1))
+print(f"\n{diff2a=}")
+
+
 h.PlotShape(False)  # PlotShape does affect diameter.
 
 rvp.to_vector(diams3,x)

The output becomes

$ python test.py

|---------|       a(0-1)
      `--------|       b(0-1)


diff=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

diff2a=[2.3841857821338408e-08, 2.3841857821338408e-08, 2.3841857821338408e-08, 2.3841858154405315e-08, 2.3841857710316106e-08, 2.3841858154405315e-08, 2.3841858043383013e-08, 2.3841857266226896e-08, 2.3841858154405315e-08, 2.3841858043383013e-08, 2.3841858043383013e-08]

diff2=[2.3841857821338408e-08, 2.3841857821338408e-08, 2.3841857821338408e-08, 2.3841858154405315e-08, 2.3841857710316106e-08, 2.3841858154405315e-08, 2.3841858043383013e-08, 2.3841857266226896e-08, 2.3841858154405315e-08, 2.3841858043383013e-08, 2.3841858043383013e-08]

I propose to mention the 32bit floating point truncation in https://nrn.readthedocs.io/en/latest/python/modelspec/programmatic/topology/geometry.html#define_shape
and also refer from there to more serious interpolation issues mentioned in https://nrn.readthedocs.io/en/latest/python/modelspec/programmatic/topology/geometry.html#pt3dconst
Also add references to define_shape in Shape and PlotShape

Lastly, in case your 3-d shapes involve contours, you should be aware of the pull request #2981 bug fix. That is currently languishing due to CI failures which are likely to require re-creation of some of the standard comparison data.