Fix edge case in LaminarCurrentSourceDensity, use Poetry for builds, Python >= 3.10

.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

.github/workflows/python-publish.yml

@@ -13,26 +13,27 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - name: Set up Python
-      uses: actions/setup-python@v3
+      uses: actions/setup-python@v5
       with:
-        python-version: '3.x'
+        python-version: '3.10'
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install setuptools twine wheel
+        pip install poetry twine
+    - name: Build package
+      run: |
+        poetry build
     - name: Build and publish to Test PYPI
       env:
         TWINE_USERNAME: ${{ secrets.TEST_PYPI_USERNAME }}
         TWINE_PASSWORD: ${{ secrets.TEST_PYPI_PASSWORD }}
       run: |
-        python setup.py sdist bdist_wheel
-        twine upload  --repository-url https://test.pypi.org/legacy/ dist/* --verbose
-    - name: Build and publish
+        twine upload --repository-url https://test.pypi.org/legacy/ dist/* --verbose
+    - name: Publish to PYPI
       env:
         TWINE_USERNAME: ${{ secrets.PYPI_USERNAME }}
         TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }}
       run: |
-        python setup.py sdist bdist_wheel
         twine upload dist/* --verbose

.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/

README.md

@@ -1,9 +1,11 @@
 # LFPykit
 
-This Python module contain freestanding implementations of electrostatic
+The online documentation of `LFPykit` can be found here:
+https://lfpykit.readthedocs.io/en/latest
+
+`LFPykit` is a freestanding implementations of electrostatic
 forward models incorporated in LFPy
 (https://github.com/LFPy/LFPy, https://LFPy.readthedocs.io).
-
 The aim of the `LFPykit` module is to provide electrostatic models
 in a manner that facilitates forward-model predictions of extracellular
 potentials and related measures from multicompartment neuron models, but
@@ -266,7 +268,7 @@ https://lfpykit.readthedocs.io/en/latest
 
 ## Dependencies
 
-`LFPykit` is implemented in Python and is written (and continuously tested) for `Python >= 3.7`.
+`LFPykit` is implemented in Python and is written (and continuously tested) for `Python >= 3.10`.
 The main `LFPykit` module depends on `numpy`, `scipy` and `MEAutility` (https://github.com/alejoe91/MEAutility, https://meautility.readthedocs.io/en/latest/).
 
 Running all unit tests and example files may in addition require `py.test`, `matplotlib`,

lfpykit/models.py

@@ -17,6 +17,7 @@
 import sys
 from copy import deepcopy
 import numpy as np
+import warnings
 from . import lfpcalc
 import MEAutility as mu
 
@@ -1994,7 +1995,13 @@ def __init__(self, cell, z, r):
         assert r.ndim == 1, 'r.ndim != 1'
         assert r.shape[0] == z.shape[0], 'r.shape[0] != z.shape[0]'
         assert np.all(r > 0), 'r must be greater than 0'
-
+        if np.max(np.diff(z, axis=-1)) > np.max(cell.length):
+            warn_msg = (
+                "All compartments should be shorter than CSD resolution. "
+                "Consider increasing cell spatial resolution or "
+                "decreasing CSD resolution."
+            )
+            warnings.warn(warn_msg, UserWarning)
         self.z = z
         self.r = r
 
@@ -2063,8 +2070,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 +2113,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 +2149,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,54 @@ 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,27 @@
+[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"
+pytest = {version = "*", optional = true}
+sympy = {version = "*", optional = true}
+sphinx = {version = "*", optional = true}
+numpydoc = {version = "*", optional = true}
+sphinx_rtd_theme = {version = "*", optional = true}
+recommonmark = {version = "*", optional = true}
+
+[tool.poetry.extras]
+tests = ["pytest", "sympy"]
+docs = ["sphinx", "numpydoc", "sphinx_rtd_theme", "recommonmark"]
+
+[tool.poetry.group.dev.dependencies]