From e4cffdd882b64e18b42250a587e68bdc61c392d5 Mon Sep 17 00:00:00 2001 From: Christian Behrens Date: Fri, 14 Jun 2019 16:00:31 +0200 Subject: [PATCH 1/5] First test of compartment condensation --- brian2/spatialneuron/morphology.py | 105 +++++++++++++++++++++++++---- 1 file changed, 93 insertions(+), 12 deletions(-) diff --git a/brian2/spatialneuron/morphology.py b/brian2/spatialneuron/morphology.py index e45cf7ce5..a5fd0ea6c 100644 --- a/brian2/spatialneuron/morphology.py +++ b/brian2/spatialneuron/morphology.py @@ -1092,6 +1092,25 @@ def _replace_three_point_soma(compartment, all_compartments): Morphology._replace_three_point_soma(child, all_compartments) + def condense(self, lam=0.1): + ''' + Condense the morphology section by section, controlled by parameter lambda + Has to be a recursive process, probably call a condense function in each section? + problems to solve: + - get mapper between old/new compartment indices + - create new morphology or change old one? + :param lam: + :return: + ''' + + if type(self) == Soma: + for c in self._children: + c.condense_section(lam) + else: + self.condense_section(lam) + + return self + @staticmethod def from_points(points, spherical_soma=True): ''' @@ -1808,6 +1827,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 @@ -1839,6 +1866,56 @@ 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): + 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): + # print(np.sqrt(self._area[i] / (self._r_length_2[i] + self._r_length_1[i + 1]) / meter)) + 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): + # print('condensing \n', self.n) + self.condensation_update(i) + # print('condensed \n', self.n) + break + if n == self.n: + run_condensation = False + for c in self._children: + c.condense_section(lam) + + def condensation_update(self, i): + ''' + condenses compartment i+1 into compartment i, it's not important which one is the origin and which the target + :param i: + :return: + ''' + # + # attention: need to catch cases where i-1,i+2 do not exist! + # calculate new r_length_1: distribute internal resistance between + # proportionally to r_m_2/(r_m_1 + r_m_2) + # check formulas! + 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 + 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:])) + self._n = self._n - 1 + @property def area(self): r''' @@ -1850,9 +1927,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): @@ -1864,9 +1942,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): @@ -1922,9 +2001,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): @@ -1933,9 +2013,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): From 42cdb39eb3088849d27e1734c7b0582feea8fba6 Mon Sep 17 00:00:00 2001 From: Christian Behrens Date: Fri, 14 Jun 2019 16:40:03 +0200 Subject: [PATCH 2/5] Fixed unit error --- brian2/spatialneuron/morphology.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/brian2/spatialneuron/morphology.py b/brian2/spatialneuron/morphology.py index a5fd0ea6c..c4e095d47 100644 --- a/brian2/spatialneuron/morphology.py +++ b/brian2/spatialneuron/morphology.py @@ -1913,7 +1913,8 @@ def condensation_update(self, i): 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:])) + 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 From b700ae35519071cf83ff085e15f6a251400675d4 Mon Sep 17 00:00:00 2001 From: chbehrens Date: Thu, 27 Jun 2019 16:30:21 +0200 Subject: [PATCH 3/5] Condense now returns a map from old to new compartments --- brian2/spatialneuron/morphology.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/brian2/spatialneuron/morphology.py b/brian2/spatialneuron/morphology.py index c4e095d47..60223dd45 100644 --- a/brian2/spatialneuron/morphology.py +++ b/brian2/spatialneuron/morphology.py @@ -1092,6 +1092,12 @@ 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 lambda @@ -1102,14 +1108,21 @@ def condense(self, lam=0.1): :param lam: :return: ''' + original_compartment_indices = self.index_structure() if type(self) == Soma: + local_index_maps = [np.array([[0], [0]])] for c in self._children: - c.condense_section(lam) + local_index_maps += c.condense_section(lam) else: - self.condense_section(lam) - - return self + local_index_maps = self.condense_section(lam) + new_compartment_indices = self.index_structure() + final_index_list = np.empty([2, 0]) + for s in range(self.total_sections): + section_index_list = np.vstack((original_compartment_indices[s], + new_compartment_indices[s][local_index_maps[s][1, :]])) + final_index_list = np.concatenate((final_index_list, section_index_list),axis=1) + return final_index_list.astype(int) @staticmethod def from_points(points, spherical_soma=True): @@ -1867,6 +1880,7 @@ def copy_section(self): length=length, type=self.type) def condense_section(self, lam): + 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 @@ -1876,17 +1890,17 @@ def condense_section(self, lam): if n == 1: break for i in range(n - 1): - # print(np.sqrt(self._area[i] / (self._r_length_2[i] + self._r_length_1[i + 1]) / meter)) 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): - # print('condensing \n', self.n) self.condensation_update(i) - # print('condensed \n', self.n) + 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: - c.condense_section(lam) + index_maps += c.condense_section(lam) + return index_maps def condensation_update(self, i): ''' From 15db0342b79947958532b5ed1a32b89ad71b779a Mon Sep 17 00:00:00 2001 From: Christian Behrens Date: Fri, 5 Jul 2019 13:27:44 +0200 Subject: [PATCH 4/5] Added tests for condensation --- brian2/spatialneuron/morphology.py | 27 +++++-- brian2/tests/test_morphology.py | 110 +++++++++++++++++++++++++++++ 2 files changed, 130 insertions(+), 7 deletions(-) diff --git a/brian2/spatialneuron/morphology.py b/brian2/spatialneuron/morphology.py index 60223dd45..f79bec5ea 100644 --- a/brian2/spatialneuron/morphology.py +++ b/brian2/spatialneuron/morphology.py @@ -1923,9 +1923,10 @@ def condensation_update(self, i): (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 - self._x = np.delete(self._x, i + 1) - self._y = np.delete(self._y, i + 1) - self._z = np.delete(self._z, i + 1) + 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 @@ -2231,6 +2232,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]) @@ -2266,7 +2275,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): @@ -2298,7 +2308,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): @@ -2307,7 +2318,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): @@ -2316,4 +2328,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 diff --git a/brian2/tests/test_morphology.py b/brian2/tests/test_morphology.py index 6ed78a887..cef5a2f2f 100644 --- a/brian2/tests/test_morphology.py +++ b/brian2/tests/test_morphology.py @@ -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 @@ -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() From 46b5264caa5375ae13084508b5c1049e0c9e019e Mon Sep 17 00:00:00 2001 From: Christian Behrens Date: Fri, 5 Jul 2019 14:52:29 +0200 Subject: [PATCH 5/5] Updated description of condensation functions --- brian2/spatialneuron/morphology.py | 68 +++++++++++++++++++++--------- 1 file changed, 48 insertions(+), 20 deletions(-) diff --git a/brian2/spatialneuron/morphology.py b/brian2/spatialneuron/morphology.py index f79bec5ea..513220e8a 100644 --- a/brian2/spatialneuron/morphology.py +++ b/brian2/spatialneuron/morphology.py @@ -1100,13 +1100,25 @@ def index_structure(self): def condense(self, lam=0.1): ''' - Condense the morphology section by section, controlled by parameter lambda - Has to be a recursive process, probably call a condense function in each section? - problems to solve: - - get mapper between old/new compartment indices - - create new morphology or change old one? - :param lam: - :return: + 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() @@ -1117,12 +1129,12 @@ def condense(self, lam=0.1): else: local_index_maps = self.condense_section(lam) new_compartment_indices = self.index_structure() - final_index_list = np.empty([2, 0]) + global_index_map = np.empty([2, 0]) for s in range(self.total_sections): - section_index_list = np.vstack((original_compartment_indices[s], - new_compartment_indices[s][local_index_maps[s][1, :]])) - final_index_list = np.concatenate((final_index_list, section_index_list),axis=1) - return final_index_list.astype(int) + 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): @@ -1880,6 +1892,23 @@ def copy_section(self): 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 @@ -1904,15 +1933,14 @@ def condense_section(self, lam): def condensation_update(self, i): ''' - condenses compartment i+1 into compartment i, it's not important which one is the origin and which the target - :param i: - :return: + 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 ''' - # - # attention: need to catch cases where i-1,i+2 do not exist! - # calculate new r_length_1: distribute internal resistance between - # proportionally to r_m_2/(r_m_1 + r_m_2) - # check formulas! 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])