Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

change volume method #2781

Merged
merged 17 commits into from May 6, 2024
Merged
20 changes: 18 additions & 2 deletions share/lib/python/neuron/rxd/region.py
Expand Up @@ -344,8 +344,14 @@
def _short_repr(self):
return "Extracellular"

def volume(self, index):
def volume(self, index=None):
"""Returns the volume of the voxel at a given index"""
if index is None:
if numpy.isscalar(self.alpha):
vol = self._nx * self._ny * self._nz * numpy.prod(self._dx) * self.alpha

Check warning on line 351 in share/lib/python/neuron/rxd/region.py

View check run for this annotation

Codecov / codecov/patch

share/lib/python/neuron/rxd/region.py#L350-L351

Added lines #L350 - L351 were not covered by tests
else:
vol = numpy.sum(self.alpha) * numpy.prod(self._dx)
return vol

Check warning on line 354 in share/lib/python/neuron/rxd/region.py

View check run for this annotation

Codecov / codecov/patch

share/lib/python/neuron/rxd/region.py#L353-L354

Added lines #L353 - L354 were not covered by tests
if numpy.isscalar(self.alpha):
return numpy.prod(self._dx) * self.alpha
return numpy.prod(self._dx) * self.alpha[index]
Expand Down Expand Up @@ -1000,6 +1006,16 @@
else:
raise RxDException("Cannot set secs now; model already instantiated")

def volume(self, index):
def volume(self, index=None):
"""Returns the volume of the voxel at a given index"""
initializer._do_init()
if index is None:
vol = 0
if hasattr(self, "_vol") and any(self._secs3d):
vol += numpy.sum(self._vol)
if hasattr(self, "_geometry") and any(self._secs1d):
vol += numpy.sum(
[self._geometry.volumes1d(sec) for sec in self._secs1d]
)
return vol
return self._vol[index]
26 changes: 14 additions & 12 deletions share/lib/python/neuron/rxd/rxd.py
Expand Up @@ -905,25 +905,27 @@ def _setup_matrices():
num_3d_indices_per_1d_seg = numpy.asarray(
num_3d_indices_per_1d_seg, dtype=numpy.int64
)
hybrid_grid_ids.append(-1)
hybrid_grid_ids = numpy.asarray(hybrid_grid_ids, dtype=numpy.int64)

hybrid_indices3d = numpy.asarray(hybrid_indices3d, dtype=numpy.int64)
rates = numpy.asarray(rates, dtype=float)
volumes1d = numpy.asarray(volumes1d, dtype=float)
volumes3d = numpy.asarray(volumes3d, dtype=float)
dxs = numpy.asarray(grids_dx, dtype=float)
set_hybrid_data(
num_1d_indices_per_grid,
num_3d_indices_per_grid,
hybrid_indices1d,
hybrid_indices3d,
num_3d_indices_per_1d_seg,
hybrid_grid_ids,
rates,
volumes1d,
volumes3d,
dxs,
)
if hybrid_grid_ids.size > 1:
set_hybrid_data(
num_1d_indices_per_grid,
num_3d_indices_per_grid,
hybrid_indices1d,
hybrid_indices3d,
num_3d_indices_per_1d_seg,
hybrid_grid_ids,
rates,
volumes1d,
volumes3d,
dxs,
)

# TODO: Replace this this to handle 1d/3d hybrid models
"""
Expand Down
20 changes: 20 additions & 0 deletions test/rxd/test_region_volume.py
@@ -0,0 +1,20 @@
def test_region_volume(neuron_nosave_instance):
h, rxd, _ = neuron_nosave_instance
dend1 = h.Section("dend1")
dend2 = h.Section("dend2")
diff = 1e-10
dend1.nseg = 4
dend2.nseg = 4
cyt1 = rxd.Region(dend1.wholetree(), "i", dx=0.25)
cyt2 = rxd.Region(dend2.wholetree(), "i", dx=0.25)
ca1 = rxd.Species(cyt1, name="ca1", charge=2, initial=1e-12)
ca2 = rxd.Species(cyt2, name="ca2", charge=2, initial=1e-12)
dend1.L = 10
dend1.diam = 2
dend2.L = 10
dend2.diam = 2
rxd.set_solve_type(domain=[dend1], dimension=1)
rxd.set_solve_type(domain=[dend2], dimension=3)
h.finitialize()
assert abs(cyt1.volume() - sum(ca1.nodes(cyt1).volume)) < diff
assert abs(cyt2.volume() - sum(ca2.nodes(cyt2).volume)) < diff