Merge f423688 into 8ce9964

lfpykit/eegmegcalc.py

@@ -172,10 +172,11 @@ 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 '
-                               'brain center than dipole: r > rz.'
-                               'r = %s, rz = %s', self.r, self._rz)
+        if any(r <= self._rz for r in self.r):
+            mssg = ('Electrode must be farther away from ' +
+                    'brain center than dipole: r > rz.' +
+                    f'r = {self.r}, rz = {self._rz}')
+            raise RuntimeError(mssg)
 
         # compute theta angle between rzloc and rxyz
         self._theta = self._calc_theta()
@@ -197,27 +198,27 @@ def get_dipole_potential(self, p, dipole_location):
         Returns
         -------
         potential: ndarray, dtype=float
-            Shape (n_contacts, n_timesteps) array containing the electric
-            potential at contact point(s) FourSphereVolumeConductor.r in units
-            of (mV) for all timesteps of current dipole moment p.
+            Shape ``(n_contacts, n_timesteps)`` array containing the electric
+            potential at contact point(s) ``FourSphereVolumeConductor.rxyz``
+            in units of (mV) for all timesteps of current dipole moment ``p``.
 
         """
 
         self._rz_params(dipole_location)
         n_contacts = self.r.shape[0]
         n_timesteps = p.shape[1]
 
-        if np.linalg.norm(p) != 0:
+        if np.linalg.norm(p) > 0:
             p_rad, p_tan = self._decompose_dipole(p)
         else:
             p_rad = np.zeros((3, n_timesteps))
             p_tan = np.zeros((3, n_timesteps))
-        if np.linalg.norm(p_rad) != 0.:
+        if np.linalg.norm(p_rad) > 0:
             pot_rad = self._calc_rad_potential(p_rad)
         else:
             pot_rad = np.zeros((n_contacts, n_timesteps))
 
-        if np.linalg.norm(p_tan) != 0.:
+        if np.linalg.norm(p_tan) > 0:
             pot_tan = self._calc_tan_potential(p_tan)
         else:
             pot_tan = np.zeros((n_contacts, n_timesteps))
@@ -228,8 +229,8 @@ def get_dipole_potential(self, p, dipole_location):
     def get_transformation_matrix(self, dipole_location):
         '''
         Get linear response matrix mapping current dipole moment in (nA µm)
-        located in location `rz` to extracellular potential in (mV)
-        at recording sites `FourSphereVolumeConductor.` (µm)
+        located in location ``rz`` to extracellular potential in (mV)
+        at recording sites ``FourSphereVolumeConductor.rxyz`` (µm)
 
         parameters
         ----------
@@ -371,8 +372,11 @@ def _calc_theta(self):
         """
         cos_theta = (self.rxyz @ self._rzloc) / (
             np.linalg.norm(self.rxyz, axis=1) * np.linalg.norm(self._rzloc))
-        theta = np.arccos(cos_theta)
-        return theta
+        # avoid RuntimeWarning: invalid value encountered in arccos
+        # and resulting NaNs
+        with np.errstate(invalid='ignore'):
+            theta = np.arccos(cos_theta)
+        return np.nan_to_num(theta)
 
     def _calc_phi(self, p_tan):
         """
@@ -396,7 +400,8 @@ def _calc_phi(self, p_tan):
         """
 
         # project rxyz onto z-axis (rzloc)
-        proj_rxyz_rz = self.rxyz * self._z
+        proj_rxyz_rz = np.outer(self.rxyz @ self._rzloc / self._rz, self._z)
+
         # find projection of rxyz in xy-plane
         rxy = self.rxyz - proj_rxyz_rz
         # define x-axis
@@ -417,7 +422,9 @@ def _calc_phi(self, p_tan):
                        np.linalg.norm(x, axis=1))[mask]
 
         # compute phi in [0, pi]
-        phi[mask] = np.arccos(cos_phi[mask])
+        with np.errstate(invalid='ignore'):
+            phi[mask] = np.nan_to_num(np.arccos(cos_phi[mask]))
+
 
         # nb: phi in [-pi, pi]. since p_tan defines direction of y-axis,
         # phi < 0 when rxy*p_tan < 0

lfpykit/models.py

@@ -342,7 +342,7 @@ class LineSourcePotential(LinearModel):
 
         - the extracellular conductivity :math:`\\sigma` is infinite,
           homogeneous, frequency independent (linear) and isotropic
-        - each segment is treated as a straigh line source with homogeneous
+        - each segment is treated as a straight line source with homogeneous
           current density between its start and end point coordinate
         - each measurement site :math:`\\mathbf{r}_j = (x_j, y_j, z_j)` is
           treated as a point

lfpykit/tests/test_eegmegcalc.py

@@ -677,6 +677,141 @@ def test_get_dipole_potential_02(self):
             else:
                 np.testing.assert_equal(phi0, phi[0][0])
 
+    def test_get_dipole_potential_03(self):
+        # check that predictions are rotation invariant
+        radii = [79000., 80000., 85000., 90000.]  # (µm)
+        sigmas = [0.3, 1.5, 0.015, 0.3]  # (S/m)
+
+        # locations in xz-plane along outer layer surface
+        r_e = radii[-1] - 1  # radius for prediction sites (µm)
+        theta = np.linspace(0, 2 * np.pi, 72, endpoint=False)  # polar (rad)
+        phi = 0  # azimuth angle (rad)
+        r_el = r_e * np.c_[np.sin(theta) * np.cos(phi),
+                           np.sin(theta) * np.sin(phi),
+                           np.cos(theta)]
+        sphere_model = eegmegcalc.FourSphereVolumeConductor(
+            r_electrodes=r_el,
+            radii=radii,
+            sigmas=sigmas)
+
+        # radial dipole locations in xz-plane
+        r = radii[0] - 1000  # dipole location(µm)
+        theta_p = np.linspace(0, 2 * np.pi, 8, endpoint=False)  # polar(rad)
+        phi_p = 0  # azimuth angle (rad)
+        r_p = r * np.c_[np.sin(theta_p) * np.cos(phi_p),
+                        np.sin(theta_p) * np.sin(phi_p),
+                        np.cos(theta_p)]
+
+        # unit radial current dipoles at each location:
+        p = (r_p.T / np.linalg.norm(r_p, axis=-1)).T  # (nAµm)
+
+        def R_y(theta=0):
+            '''rotation matrix around y-axis by some angle theta (rad)'''
+            return np.c_[[np.cos(theta), 0, np.sin(theta)],
+                         [0, 1, 0],
+                         [-np.sin(theta), 0, np.cos(theta)]].T
+
+        R_y_45 = R_y(theta=np.pi / 4)  # rotate by 45 deg
+
+        V_e = np.zeros((theta_p.size, theta.size))
+        for i, (r_p_, p_, theta_p_) in enumerate(
+                zip(r_p, p @ R_y_45, theta_p)):
+            V_e[i] = np.roll(
+                sphere_model.get_dipole_potential(np.expand_dims(p_, -1),
+                                                  r_p_).ravel(),
+                -i * theta.size // theta_p.size)
+        assert np.allclose(V_e, V_e[0])
+
+    def test_get_dipole_potential_04(self):
+        # check that predictions are rotation invariant
+        radii = [79000., 80000., 85000., 90000.]  # (µm)
+        sigmas = [0.3, 1.5, 0.015, 0.3]  # (S/m)
+
+        # locations in xy-plane along outer layer surface
+        r_e = radii[-1] - 1  # radius for prediction sites (µm)
+        theta = np.linspace(0, 2 * np.pi, 72, endpoint=False)  # polar (rad)
+        phi = 0  # azimuth angle (rad)
+        r_el = r_e * np.c_[np.sin(theta) * np.cos(phi),
+                           np.cos(theta),
+                           np.sin(theta) * np.sin(phi)]
+        sphere_model = eegmegcalc.FourSphereVolumeConductor(
+            r_electrodes=r_el,
+            radii=radii,
+            sigmas=sigmas)
+
+        # radial dipole locations in xz-plane
+        r = radii[0] - 1000  # dipole location(µm)
+        theta_p = np.linspace(0, 2 * np.pi, 8, endpoint=False)  # polar(rad)
+        phi_p = 0  # azimuth angle (rad)
+        r_p = r * np.c_[np.sin(theta_p) * np.cos(phi_p),
+                        np.cos(theta_p),
+                        np.sin(theta_p) * np.sin(phi_p)]
+
+        # unit radial current dipoles at each location:
+        p = (r_p.T / np.linalg.norm(r_p, axis=-1)).T  # (nAµm)
+
+        def R_z(theta=0):
+            '''rotation matrix around z-axis by some angle theta (rad)'''
+            return np.c_[[np.cos(theta), -np.sin(theta), 0],
+                         [np.sin(theta), np.cos(theta), 0],
+                         [0, 0, 1]].T
+
+        R_z_45 = R_z(theta=np.pi / 4)  # rotate by 45 deg
+
+        V_e = np.zeros((theta_p.size, theta.size))
+        for i, (r_p_, p_, theta_p_) in enumerate(
+                zip(r_p, p @ R_z_45, theta_p)):
+            V_e[i] = np.roll(
+                sphere_model.get_dipole_potential(np.expand_dims(p_, -1),
+                                                  r_p_).ravel(),
+                -i * theta.size // theta_p.size)
+        assert np.allclose(V_e, V_e[0])
+
+    def test_get_dipole_potential_05(self):
+        # check that predictions are rotation invariant
+        radii = [79000., 80000., 85000., 90000.]  # (µm)
+        sigmas = [0.3, 1.5, 0.015, 0.3]  # (S/m)
+
+        # locations in yz-plane along outer layer surface
+        r_e = radii[-1] - 1  # radius for prediction sites (µm)
+        theta = np.linspace(0, 2 * np.pi, 72, endpoint=False)  # polar (rad)
+        phi = 0  # azimuth angle (rad)
+        r_el = r_e * np.c_[np.sin(theta) * np.sin(phi),
+                           np.sin(theta) * np.cos(phi),
+                           np.cos(theta)]
+        sphere_model = eegmegcalc.FourSphereVolumeConductor(
+            r_electrodes=r_el,
+            radii=radii,
+            sigmas=sigmas)
+
+        # radial dipole locations in yz-plane
+        r = radii[0] - 1000  # dipole location(µm)
+        theta_p = np.linspace(0, 2 * np.pi, 8, endpoint=False)  # polar(rad)
+        phi_p = 0  # azimuth angle (rad)
+        r_p = r * np.c_[np.sin(theta_p) * np.sin(phi_p),
+                        np.sin(theta_p) * np.cos(phi_p),
+                        np.cos(theta_p)]
+
+        # unit radial current dipoles at each location:
+        p = (r_p.T / np.linalg.norm(r_p, axis=-1)).T  # (nAµm)
+
+        def R_x(theta=0):
+            '''rotation matrix around x-axis by some angle theta (rad)'''
+            return np.c_[[1, 0, 0],
+                         [0, np.cos(theta), -np.sin(theta)],
+                         [0, np.sin(theta), np.cos(theta)]].T
+
+        R_x_45 = R_x(theta=np.pi / 4)  # rotate by 45 deg
+
+        V_e = np.zeros((theta_p.size, theta.size))
+        for i, (r_p_, p_, theta_p_) in enumerate(
+                zip(r_p, p @ R_x_45, theta_p)):
+            V_e[i] = np.roll(
+                sphere_model.get_dipole_potential(np.expand_dims(p_, -1),
+                                                  r_p_).ravel(),
+                -i * theta.size // theta_p.size)
+        assert np.allclose(V_e, V_e[0])
+
     def test_get_transformation_matrix_00(self):
         '''Test radial and tangential parts of dipole sums to dipole'''
         radii = [88000, 90000, 95000, 100000]