From 3817e0fd9bc3610dbad61c1e722324083fc3211a Mon Sep 17 00:00:00 2001 From: solveignaess Date: Thu, 30 Jan 2020 16:35:45 +0100 Subject: [PATCH] fixed bug in decompose_dipole, eegmegcalc.py (#155) * fixed bug in decompose_dipole, added test * removed prints --- LFPy/eegmegcalc.py | 6 ++---- LFPy/test/test_eegmegcalc.py | 29 ++++++++++++++++++++++++++++- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/LFPy/eegmegcalc.py b/LFPy/eegmegcalc.py index 282e5d5a..1dd1e2f6 100644 --- a/LFPy/eegmegcalc.py +++ b/LFPy/eegmegcalc.py @@ -415,7 +415,7 @@ def _rz_params(self, rz): 'Can be avoided by placing dipole further away from brain surface.') if any(r < self._rz for r in self.r): - raise RuntimeError('Electrode must be farther away from' + raise RuntimeError('Electrode must be farther away from ' 'brain center than dipole: r > rz.' 'r = %s, rz = %s', self.r, self._rz) @@ -454,7 +454,6 @@ def calc_potential(self, p, rz): else: p_rad = np.zeros((n_timesteps, 3)) p_tan = np.zeros((n_timesteps, 3)) - if np.linalg.norm(p_rad) != 0.: pot_rad = self._calc_rad_potential(p_rad) else: @@ -542,9 +541,8 @@ def _decompose_dipole(self, p): Shape (n_timesteps, 3) array, tangential part of p, orthogonal to self._rz """ - p_rad = p*self._z + p_rad = np.dot(p, self._z).reshape(len(p),1)*self._z.reshape(1, len(self._z)) p_tan = p - p_rad - return p_rad, p_tan def _calc_rad_potential(self, p_rad): diff --git a/LFPy/test/test_eegmegcalc.py b/LFPy/test/test_eegmegcalc.py index 161a3d08..313c3e3b 100644 --- a/LFPy/test/test_eegmegcalc.py +++ b/LFPy/test/test_eegmegcalc.py @@ -206,12 +206,39 @@ def test_check_params02(self): with np.testing.assert_raises(ValueError): LFPy.FourSphereVolumeConductor(radii, sigmas, r_el2) - def test_decompose_dipole(self): + def test_decompose_dipole01(self): '''Test radial and tangential parts of dipole sums to dipole''' P1 = np.array([[1., 1., 1.]]) p_rad, p_tan = decompose_dipole(P1) np.testing.assert_equal(p_rad + p_tan, P1) + def test_decompose_dipole02(self): + '''Test radial and tangential parts of dipole sums to dipole''' + radii = [88000, 90000, 95000, 100000] + sigmas = [0.3, 1.5, 0.015, 0.3] + ps = np.array([[ 1000., 0., 0.], + [-1000., 0., 0.], + [0., 1000., 0.], + [0., -1000., 0.], + [0., 0., 1000.], + [0., 0., -1000.], + [10., 20., 30.], + [-10., -20., -30.]]) + p_locs = np.array([[ 87000., 0., 0.], + [-87000., 0., 0.], + [0., 87000., 0.], + [0., -87000., 0.], + [0., 0., 87000.], + [0., 0., -87000.], + [80000., 2000., 3000.], + [-2000., -80000., -3000.]]) + el_locs = np.array([[90000., 5000., -5000.]]) + fs = LFPy.FourSphereVolumeConductor(radii, sigmas, el_locs) + for p_loc in p_locs: + fs._rz_params(p_loc) + p_rads, p_tans = fs._decompose_dipole(ps) + np.testing.assert_equal(p_rads + p_tans, ps) + def test_rad_dipole(self): '''Test that radial part of decomposed dipole is correct''' P1 = np.array([[1., 1., 1.]])