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.

Sign up for GitHub

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

Jump to bottom

Compartment condensation #1108

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Compartment condensation #1108

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e4cffdd
First test of compartment condensation
chbehrens Jun 14, 2019
42cdb39
Fixed unit error
chbehrens Jun 14, 2019
b700ae3
Condense now returns a map from old to new compartments
chbehrens Jun 27, 2019
15db034
Added tests for condensation
chbehrens Jul 5, 2019
46b5264
Updated description of condensation functions
chbehrens Jul 5, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
169 changes: 153 additions & 16 deletions brian2/spatialneuron/morphology.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1092,6 +1092,50 @@ def _replace_three_point_soma(compartment, all_compartments):
Morphology._replace_three_point_soma(child,
all_compartments)

def index_structure(self):
index_list = [self.indices[:]]
for c in self.children:
index_list += c.index_structure()
return index_list

def condense(self, lam=0.1):
'''
Condense the morphology section by section, controlled by parameter lam. The
recursive condensation process is taken care of by the function
`Section.condense_section`.

Once the condensation is done a map between old and new compartment indices
is assembled that can be used to set variables of a spatial neuron compartment-
wise without knowing the new indices.

Parameters
----------
lam : float
Parameter defining the degree of condensation. Small lam leaves
more compartments, larger lam results in stronger condensation.

Returns
-------
global_index_map : numpy.array
A 2 x self.n array containing the indices of all original compartments
(0 to n) and the ones of the corresponding new compartments.
'''
original_compartment_indices = self.index_structure()

if type(self) == Soma:
local_index_maps = [np.array([[0], [0]])]
for c in self._children:
local_index_maps += c.condense_section(lam)
else:
local_index_maps = self.condense_section(lam)
new_compartment_indices = self.index_structure()
global_index_map = np.empty([2, 0])
for s in range(self.total_sections):
section_index_map = np.vstack((original_compartment_indices[s],
new_compartment_indices[s][local_index_maps[s][1, :]]))
global_index_map = np.concatenate((global_index_map, section_index_map),axis=1)
return global_index_map.astype(int)

@staticmethod
def from_points(points, spherical_soma=True):
'''
Expand Down Expand Up @@ -1808,6 +1852,14 @@ def __init__(self, diameter, n=1, length=None, x=None, y=None, z=None,
(self.end_z - self.start_z) ** 2)
self._length = length

d_1 = self.start_diameter
d_2 = self.end_diameter
d_mid = 0.5 * (d_1 + d_2)
self._area = np.pi / 2 * (d_1 + d_2) * np.sqrt(((d_1 - d_2) ** 2) / 4 + self._length ** 2)
self._volume = np.pi * self._length * (d_1 ** 2 + d_1 * d_2 + d_2 ** 2) / 12
self._r_length_1 = np.pi / 2 * (d_1 * d_mid) / self._length
self._r_length_2 = np.pi / 2 * (d_mid * d_2) / self._length

def __repr__(self):
if all(np.abs(self.end_diameter - self.end_diameter[0]) < self.end_diameter[0]*1e-12):
# Constant diameter
Expand Down Expand Up @@ -1839,6 +1891,75 @@ def copy_section(self):
return Section(diameter=self._diameter, n=self.n, x=x, y=y, z=z,
length=length, type=self.type)

def condense_section(self, lam):
'''
Function to recursively condense compartments section by section.
It Condenses compartments in the current section and continues with
all child sections.

Parameters
----------
lam : float
Parameter defining the degree of condensation. Small lam leaves
more compartments, larger lam results in stronger condensation.
Returns
-------
index_maps : list
List that maps original local compartment indices (0, 1, ...)
to the new condensed compartments. It is used in the function
`condense` to assemble a global map of compartment indices.
'''
section_index_map = np.tile(np.arange(self.n), (2, 1))
for lamcrit in np.array([0.2, 0.4, 0.6, 0.8, 1.0]) * lam:
run_condensation = True
k = 0
while run_condensation:
k += 1
n = self.n
if n == 1:
break
for i in range(n - 1):
if (np.sqrt(self._area[i] / (self._r_length_2[i] + self._r_length_1[i + 1]) / meter) < lamcrit) \
or (np.sqrt(self._area[i + 1] / (self._r_length_2[i] + self._r_length_1[i + 1]) / meter) < lamcrit):
self.condensation_update(i)
section_index_map[1, section_index_map[1, :] > i] -= 1
break
if n == self.n:
run_condensation = False
index_maps = [section_index_map]
for c in self._children:
index_maps += c.condense_section(lam)
return index_maps

def condensation_update(self, i):
'''
Condense two compartments by merging compartment i and i+1 and
calculate parameters of the new compartment.

Parameters
----------
i : int
Index of the compartment within the section
'''
new_r_length_1 = 1 / (1 / self._r_length_1[i] + (1 / self._r_length_2[i] + 1 / self._r_length_1[i + 1])
* self._area[i] / (self._area[i] + self._area[i + 1]))
new_r_length_2 = 1 / (1 / self._r_length_2[i + 1] + (1 / self._r_length_2[i] + 1 / self._r_length_1[i + 1])
* self._area[i + 1] / (self._area[i] + self._area[i + 1]))
self._area = np.concatenate(
(self._area[:i], [self._area[i] + self._area[i + 1]], self._area[i + 2:])) * meter ** 2
self._volume = np.concatenate(
(self._volume[:i], [self._volume[i] + self._volume[i + 1]], self._volume[i + 2:])) * meter ** 3
self._r_length_1 = np.concatenate((self._r_length_1[:i], [new_r_length_1], self._r_length_1[i + 2:])) * meter
self._r_length_2 = np.concatenate((self._r_length_2[:i], [new_r_length_2], self._r_length_2[i + 2:])) * meter
if self._x is not None:
self._x = np.delete(self._x, i + 1)
self._y = np.delete(self._y, i + 1)
self._z = np.delete(self._z, i + 1)
self._diameter = np.delete(self._diameter, i + 1) * meter
self._length = np.concatenate((self._length[:i], [self._length[i] + self._length[i + 1]],
self._length[i + 2:])) * meter
self._n = self._n - 1

@property
def area(self):
r'''
Expand All @@ -1850,9 +1971,10 @@ def area(self):
respectively. Note that this surface area does not contain the area of
the two disks at the two sides of the truncated cone.
'''
d_1 = self.start_diameter
d_2 = self.end_diameter
return np.pi/2*(d_1 + d_2)*np.sqrt(((d_1 - d_2)**2)/4 + self._length**2)
# d_1 = self.start_diameter
# d_2 = self.end_diameter
# return np.pi / 2 * (d_1 + d_2) * np.sqrt(((d_1 - d_2) ** 2) / 4 + self._length ** 2)
return self._area

@property
def volume(self):
Expand All @@ -1864,9 +1986,10 @@ def volume(self):
:math:`d_2` are the diameter at the start and end of the compartment,
respectively.
'''
d_1 = self.start_diameter
d_2 = self.end_diameter
return np.pi * self._length * (d_1**2 + d_1*d_2 + d_2**2)/12
# d_1 = self.start_diameter
# d_2 = self.end_diameter
# return np.pi * self._length * (d_1 ** 2 + d_1 * d_2 + d_2 ** 2) / 12
return self._volume

@property
def length(self):
Expand Down Expand Up @@ -1922,9 +2045,10 @@ def r_length_1(self):
start and the midpoint of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
'''
d_1 = self.start_diameter
d_2 = (self.start_diameter + self.end_diameter)*0.5
return np.pi/2 * (d_1 * d_2)/self._length
# d_1 = self.start_diameter
# d_2 = (self.start_diameter + self.end_diameter) * 0.5
# return np.pi / 2 * (d_1 * d_2) / self._length
return self._r_length_1

@property
def r_length_2(self):
Expand All @@ -1933,9 +2057,10 @@ def r_length_2(self):
midpoint and the end of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
'''
d_1 = (self.start_diameter + self.end_diameter)*0.5
d_2 = self.end_diameter
return np.pi/2 * (d_1 * d_2)/self._length
# d_1 = (self.start_diameter + self.end_diameter) * 0.5
# d_2 = self.end_diameter
# return np.pi / 2 * (d_1 * d_2) / self._length
return self._r_length_2

@property
def start_x_(self):
Expand Down Expand Up @@ -2135,6 +2260,14 @@ def __init__(self, diameter, n=1, length=None, x=None, y=None, z=None,
(self.end_z - self.start_z) ** 2)
self._length = length

d_1 = self.start_diameter
d_2 = self.end_diameter
d_mid = 0.5 * (d_1 + d_2)
self._area = np.pi / 2 * (d_1 + d_2) * np.sqrt(((d_1 - d_2) ** 2) / 4 + self._length ** 2)
self._volume = np.pi * self._length * (d_1 ** 2 + d_1 * d_2 + d_2 ** 2) / 12
self._r_length_1 = np.pi / 2 * (d_1 * d_mid) / self._length
self._r_length_2 = np.pi / 2 * (d_mid * d_2) / self._length

def __repr__(self):
s = '{klass}(diameter={diam!r}'.format(klass=self.__class__.__name__,
diam=self.diameter[0])
Expand Down Expand Up @@ -2170,7 +2303,8 @@ def area(self):
diameter. Note that this surface area does not contain the area of
the two disks at the two sides of the cylinder.
'''
return np.pi * self._diameter * self.length
# return np.pi * self._diameter * self.length
return self._area

@property
def start_diameter(self):
Expand Down Expand Up @@ -2202,7 +2336,8 @@ def volume(self):
where :math:`l` is the length of the compartment, and :math:`d` is its
diameter.
'''
return np.pi * (self._diameter/2)**2 * self.length
# return np.pi * (self._diameter/2)**2 * self.length
return self._volume

@property
def r_length_1(self):
Expand All @@ -2211,7 +2346,8 @@ def r_length_1(self):
start and the midpoint of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
'''
return np.pi/2 * (self._diameter**2)/self.length
# return np.pi/2 * (self._diameter**2)/self.length
return self._r_length_1

@property
def r_length_2(self):
Expand All @@ -2220,4 +2356,5 @@ def r_length_2(self):
midpoint and the end of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
'''
return np.pi/2 * (self._diameter**2)/self.length
# return np.pi/2 * (self._diameter**2)/self.length
return self._r_length_2
110 changes: 110 additions & 0 deletions brian2/tests/test_morphology.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -702,6 +702,114 @@ def test_tree_soma_from_swc_3_point_soma():
_check_tree_soma(soma, coordinates=True, use_cylinders=False)


def _check_condensation(morphology, compartment_map, coordinates=False, use_cylinders=True):

# check condensed morphology
# number of compartments per section
assert morphology.n == 1
assert morphology['1'].n == 4
assert morphology['2'].n == 2

# number of compartments per subtree
assert morphology.total_compartments == 7
assert morphology['1'].total_compartments == 4
assert morphology['2'].total_compartments == 2

# number of sections per subtree
assert morphology.total_sections == 3
assert morphology['1'].total_sections == 1
assert morphology['2'].total_sections == 1

assert_allclose(morphology.diameter, [30]*um)

# Check that distances (= distance to root at midpoint)
# correctly follow the tree structure
# Note that the soma does add nothing to the distance
assert_equal(morphology.distance, 0 * um)
assert_allclose(morphology['1'].distance, [20, 50, 70, 90]*um)
assert_allclose(morphology['2'].distance, [10, 35]*um)
assert_allclose(morphology.end_distance, 0 * um)
assert_allclose(morphology['1'].end_distance, 100 * um)
assert_allclose(morphology['2'].end_distance, 50 * um)

assert_allclose(morphology.diameter, 30*um)
assert_allclose(morphology['1'].start_diameter, [8, 6, 4, 2]*um)
assert_allclose(morphology['1'].diameter, [7, 5, 3, 1]*um)
assert_allclose(morphology['1'].end_diameter, [6, 4, 2, 0]*um)
assert_allclose(morphology['2'].start_diameter, np.ones(2) * 4*um)
assert_allclose(morphology['2'].diameter, np.ones(2) * 4*um)
assert_allclose(morphology['2'].end_diameter, np.ones(2) * 4*um)

# Check additional parameters that change during condensation
assert_allclose(morphology['1'].length, [40, 20, 20, 20]*um)
assert_allclose(morphology['2'].length, [20, 30]*um)
assert_allclose(morphology['1'].area, [943.02723161, 314.55171931, 188.73103159, 62.91034386]*um**2)
assert_allclose(morphology['2'].area, [251.32741229, 376.99111843]*um**2)
assert_allclose(morphology['1'].volume, [1780.23583703, 397.93506945, 146.60765717, 20.94395102]*um**3)
assert_allclose(morphology['2'].volume, [251.32741229, 376.99111843]*um**3)
assert_allclose(morphology['1'].r_length_1, [2.34645164, 2.35619449, 0.9424778 , 0.15707963]*um)
assert_allclose(morphology['2'].r_length_1, [1.25663706, 0.62831853]*um)
assert_allclose(morphology['1'].r_length_2, [1.99112587, 1.57079633, 0.4712389, 0.]*um)
assert_allclose(morphology['2'].r_length_2, [1.25663706, 1.25663706]*um)

# Check compartment mapping
assert_equal(compartment_map, np.array([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [0, 1, 1, 2, 3, 4, 5, 5, 6, 6, 6]]))

if coordinates:
# Coordinates should be absolute
# section: soma
assert_allclose(morphology.start_x, 100*um)
assert_allclose(morphology.x, 100*um)
assert_allclose(morphology.end_x, 100*um)
assert_allclose(morphology.y, 0*um)
assert_allclose(morphology.z, 0*um)
# section: cable['1']
step = 20 / np.sqrt(2) * um
assert_allclose(morphology['1'].start_x, 100 * um + [0, 2, 3, 4] * step)
assert_allclose(morphology['1'].x, 100 * um + [0.5, 2, 3, 4] * step + step/2)
assert_allclose(morphology['1'].end_x, 100 * um + [1, 2, 3, 4] * step + step)
assert_allclose(morphology['1'].start_y, [0, 2, 3, 4] * step)
assert_allclose(morphology['1'].y, [0.5, 2, 3, 4] * step + step/2)
assert_allclose(morphology['1'].end_y, [1, 2, 3, 4] * step + step)
assert_allclose(morphology['1'].z, np.zeros(4) * um)
# section: cable['2']
step = 10 / np.sqrt(2) * um
assert_allclose(morphology['2'].start_x, 100 * um + [0, 2] * step)
if use_cylinders:
assert_allclose(morphology['2'].x, 100 * um + [0.5, 3] * step + step / 2)
assert_allclose(morphology['2'].end_x, 100 * um + [1, 4] * step + step)
assert_allclose(morphology['2'].start_y, -([0, 2] * step))
if use_cylinders:
assert_allclose(morphology['2'].y, -([0.5, 3] * step + step / 2))
assert_allclose(morphology['2'].end_y, -([1, 4] * step + step))
if use_cylinders:
assert_allclose(morphology['2'].z, np.zeros(2) * um)


@attr('codegen-independent')
def test_condensation_schematic():
soma = Soma(diameter=30*um)
soma.L = Section(n=5, diameter=[8, 8, 6, 4, 2, 0]*um,
length=np.ones(5)*20*um) # tapering truncated cones
soma.R = Cylinder(n=5, diameter=4*um, length=50*um)
compartment_map = soma.condense(0.007)

_check_condensation(soma, compartment_map)


@attr('codegen-independent')
def test_condensation_coordinates():
soma = Soma(diameter=30*um, x=100*um)
soma.L = Section(n=5, diameter=[8, 8, 6, 4, 2, 0]*um,
x=np.linspace(0, 100, 6)/np.sqrt(2)*um,
y=np.linspace(0, 100, 6)/np.sqrt(2)*um) # tapering truncated cones
soma.R = Cylinder(n=5, diameter=4*um,
x=[0, 50]*um/np.sqrt(2), y=[0, -50]*um/np.sqrt(2))
compartment_map = soma.condense(0.007)

_check_condensation(soma, compartment_map, coordinates=True)


@attr('codegen-independent')
def test_construction_incorrect_arguments():
### Morphology
Expand Down Expand Up @@ -1312,6 +1420,8 @@ def test_str_repr():
test_tree_soma_from_points_3_point_soma_incorrect()
test_tree_soma_from_swc()
test_tree_soma_from_swc_3_point_soma()
test_condensation_schematic()
test_condensation_coordinates()
test_construction_incorrect_arguments()
test_from_points_minimal()
test_from_points_incorrect()
Expand Down