Fix edge case in LaminarCurrentSourceDensity

Fixes #196
Fixes #189
Fixes #195

Author: espenhgn

.github/workflows/coveralls.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
         matrix:
             os: [ubuntu-latest]
-            python-version: [3.8]
+            python-version: [3.12]
 
     steps:
     - uses: actions/checkout@v3

.github/workflows/flake8.yml

@@ -12,7 +12,7 @@ jobs:
       - name: Set up Python environment
         uses: actions/setup-python@v3
         with:
-          python-version: "3.8"
+          python-version: "3.12"
       - name: flake8 Lint
         uses: reviewdog/action-flake8@v3
         with:

.github/workflows/python-app.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
         matrix:
             os: [ubuntu-latest]
-            python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
+            python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
 
     steps:
     - uses: actions/checkout@v3

.gitignore

@@ -130,6 +130,7 @@ dmypy.json
 
 # MacOS crap
 .DS_Store
+__MACOSX/
 
 # NEURON
 arm64/
@@ -140,3 +141,4 @@ doc/source/_build
 
 # examples
 examples/src
+L5bPCmodelsEH/

lfpykit/models.py

@@ -2063,8 +2063,34 @@ def get_transformation_matrix(self):
 
             z2 = np.array([z_i[0], z_i[1], z_i[0]])
             z3 = np.array([z_i[0], z_i[1], z_i[1]])
+
+            def _get_fraction_from_z_crossing(k, z_cross, endpoint):
+                z0 = self.cell.z[k, 0]
+                z1 = self.cell.z[k, 1]
+                dz = z1 - z0
+                if np.isclose(dz, 0.):
+                    return None
+
+                t = (z_cross - z0) / dz
+                if not (0. <= t <= 1.):
+                    return None
+
+                if endpoint == 0:
+                    return abs(t)
+                return abs(1. - t)
+
             # iterate over lower, right, upper boundary
             for k in np.where(inds)[0]:
+                # If both endpoints are inside the cylindrical radius, the
+                # in-volume fraction follows directly from z-intersection.
+                if kk1[k]:
+                    z_cross = z_i[1] if self.cell.z[k, 1] >= z_i[1] else z_i[0]
+                    frac = _get_fraction_from_z_crossing(
+                        k=k, z_cross=z_cross, endpoint=0)
+                    if frac is not None:
+                        M[i, k] = frac
+                        continue
+
                 for ll in range(3):
                     Pr, Pz, hit = _PrPz(r0=R[k, 0], z0=self.cell.z[k, 0],
                                         r1=R[k, 1], z1=self.cell.z[k, 1],
@@ -2080,6 +2106,16 @@ def get_transformation_matrix(self):
             inds = (~ll0) & ll1
 
             for k in np.where(inds)[0]:
+                # If both endpoints are inside the cylindrical radius, the
+                # in-volume fraction follows directly from z-intersection.
+                if kk0[k]:
+                    z_cross = z_i[1] if self.cell.z[k, 0] >= z_i[1] else z_i[0]
+                    frac = _get_fraction_from_z_crossing(
+                        k=k, z_cross=z_cross, endpoint=1)
+                    if frac is not None:
+                        M[i, k] = frac
+                        continue
+
                 for ll in range(3):
                     Pr, Pz, hit = _PrPz(r0=R[k, 0], z0=self.cell.z[k, 0],
                                         r1=R[k, 1], z1=self.cell.z[k, 1],
@@ -2106,15 +2142,7 @@ def _PrPz(r0, z0, r1, z1, r2, z2, r3, z3):
                - (r0 - r1) * (r2 * z3 - r3 * z2)) / denom)
         Pz = (((r0 * z1 - z0 * r1) * (z2 - z3)
                - (z0 - z1) * (r2 * z3 - r3 * z2)) / denom)
-    # check if intersection point lies on lines
-    if (Pr >= r0) & (Pr <= r1) & (Pz >= z0) & (Pz <= z1):
-        hit = True
-    elif (Pr <= r0) & (Pr >= r1) & (Pz >= z0) & (Pz <= z1):
-        hit = True
-    elif (Pr >= r0) & (Pr <= r1) & (Pz <= z0) & (Pz >= z1):
-        hit = True
-    elif (Pr <= r0) & (Pr >= r1) & (Pz <= z0) & (Pz >= z1):
-        hit = True
-    else:
-        hit = False
+    # check if intersection point lies on line segment (r0, z0) -> (r1, z1)
+    hit = (min(r0, r1) <= Pr <= max(r0, r1)
+           and min(z0, z1) <= Pz <= max(z0, z1))
     return Pr, Pz, hit

lfpykit/tests/test_module.py

@@ -732,7 +732,7 @@ def test_LaminarCurrentSourceDensity_02(self):
 
         np.testing.assert_allclose(M_gt, M)
 
-    def test_LaminarCurrentSourceDensity_3(self):
+    def test_LaminarCurrentSourceDensity_03(self):
         '''test LaminarCurrentSourceDensity implementation'''
         cell = get_cell(n_seg=3)
         cell.z = cell.z * 10
@@ -754,7 +754,7 @@ def test_LaminarCurrentSourceDensity_3(self):
 
         np.testing.assert_allclose(M_gt, M)
 
-    def test_LaminarCurrentSourceDensity_4(self):
+    def test_LaminarCurrentSourceDensity_04(self):
         '''test LaminarCurrentSourceDensity implementation'''
         cell = get_cell(n_seg=4)
         cell.z = cell.z * 10
@@ -826,6 +826,55 @@ def test_LaminarCurrentSourceDensity_06(self):
 
         np.testing.assert_allclose(M_gt, M)
 
+
+    def test_LaminarCurrentSourceDensity_07(self):
+        '''test LaminarCurrentSourceDensity implementation
+
+        Issue #196 (https://github.com/LFPy/LFPykit/issues/196)
+        '''
+        cell = lfp.CellGeometry(x=np.array([[50, 50], [50, -50]]),
+                                y=np.array([[0, 0], [0, 0]]),
+                                z=np.array([[-100, -50], [-50, 50]]),
+                                d=np.array([1, 1]))
+
+        h = 150.
+        r = 200.
+        M_gt = np.array([[1., 0.5],
+                         [0., 0.5]]) / (np.pi * r**2 * h)
+        
+        z = np.array([[-h, 0.], [0., h]])
+        r = np.array([r, r])
+
+        # instantiate electrode, get linear response matrix
+        csd_lam = lfp.LaminarCurrentSourceDensity(cell=cell, z=z, r=r)
+        M = csd_lam.get_transformation_matrix()
+
+        np.testing.assert_allclose(M_gt, M)
+
+    def test_LaminarCurrentSourceDensity_08(self):
+        '''test LaminarCurrentSourceDensity implementation
+
+        Issue #196 (https://github.com/LFPy/LFPykit/issues/196)
+        '''
+        cell = lfp.CellGeometry(x=np.array([[0, 0], [0, 100]]),
+                                y=np.array([[0, 0], [0, 0]]),
+                                z=np.array([[-100, -50], [-50, 50]]),
+                                d=np.array([1, 1]))
+
+        h = 150.
+        r = 200.
+        M_gt = np.array([[1., 0.5],
+                         [0., 0.5]]) / (np.pi * r**2 * h)
+        
+        z = np.array([[-h, 0.], [0., h]])
+        r = np.array([r, r])
+
+        # instantiate electrode, get linear response matrix
+        csd_lam = lfp.LaminarCurrentSourceDensity(cell=cell, z=z, r=r)
+        M = csd_lam.get_transformation_matrix()
+
+        np.testing.assert_allclose(M_gt, M)
+
     def test_VolumetricCurrentSourceDensity_00(self):
         cell = get_cell(n_seg=1)
         cell.z = cell.z * 10.

lfpykit/version.py

@@ -1 +1 @@
-version = "0.5.1"
+version = "0.6.0"

pyproject.toml

@@ -0,0 +1,21 @@
+[build-system]
+requires = ["poetry-core>=1.0.0"]
+build-backend = "poetry.core.masonry.api"
+
+[tool.poetry]
+name = "LFPykit"
+version = "0.6.0"
+description = "Electrostatic models for multicompartment neuron models"
+authors = ["LFPy-team <lfpy@users.noreply.github.com>"]
+
+[tool.poetry.dependencies]
+python = ">=3.10"
+numpy = ">=1.15.2"
+scipy = ">=1.5.2"
+meautility = ">=1.5.1"
+
+[tool.poetry.extras]
+tests = ["pytest", "sympy"]
+docs = ["sphinx", "numpydoc", "sphinx_rtd_theme", "recommonmark"]
+
+[tool.poetry.dev-dependencies]