added symmetry test and simplified decompose function (#156)

LFPy · Feb 10, 2020 · ae91cbd · ae91cbd
1 parent 3817e0f
commit ae91cbd
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 1 deletion.
diff --git a/LFPy/eegmegcalc.py b/LFPy/eegmegcalc.py
@@ -541,7 +541,8 @@ def _decompose_dipole(self, p):
             Shape (n_timesteps, 3) array, tangential part of p,
             orthogonal to self._rz
         """
-        p_rad = np.dot(p, self._z).reshape(len(p),1)*self._z.reshape(1, len(self._z))
+        z_ = np.expand_dims(self._z, -1) # reshape z-axis vector
+        p_rad = np.dot(np.dot(p, z_), z_.T)
         p_tan = p - p_rad
         return p_rad, p_tan
 

diff --git a/LFPy/test/test_eegmegcalc.py b/LFPy/test/test_eegmegcalc.py
@@ -370,6 +370,42 @@ def test_calc_potential01(self):
         global_error = np.abs(pot_analytical - pot_fem)/(np.max(np.abs(pot_fem)))
         np.testing.assert_array_less(global_error, 0.01)
 
+    def test_calc_potential02(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]
+
+        dips = np.array([[[ 1000., 0., 0.]],
+                         [[-1000., 0., 0.]],
+                         [[0.,  1000., 0.]],
+                         [[0., -1000., 0.]],
+                         [[0., 0.,  1000.]],
+                         [[0., 0., -1000.]]])
+
+        p_locs = np.array([[ 87000., 0., 0.],
+                           [-87000., 0., 0.],
+                           [0.,  87000., 0.],
+                           [0., -87000., 0.],
+                           [0., 0.,  87000],
+                           [0., 0., -87000]])
+
+        el_locs = np.array([[[ 99000., 0., 0.]],
+                            [[-99000., 0., 0.]],
+                            [[0.,  99000., 0.]],
+                            [[0., -99000., 0.]],
+                            [[0., 0.,  99000.]],
+                            [[0., 0., -99000.]]])
+
+
+        for i in range(len(p_locs)):
+            fs = LFPy.FourSphereVolumeConductor(radii, sigmas, el_locs[i])
+            phi = fs.calc_potential(dips[i], p_locs[i])
+            if i == 0:
+                phi0 = phi[0][0]
+            else:
+                np.testing.assert_equal(phi0, phi[0][0])
+
+
     def test_calc_potential_from_multi_dipoles00(self):
         """test comparison between multi-dipoles and single dipole approach"""
         neuron.h('forall delete_section()')