Fix kvm model numpy version compatibility

CHANGELOG

@@ -1,3 +1,5 @@
+0.1.1
+ - fix: KVM hertz model compatible with numpy version 1 and 2 (#1, #5)
 0.1.0
  - feat: add KVM hertz model (#2)
  - setup: initial commit and package creation

nanite_community/models/model_hertz_corrected_viscoelasticity_KVM.py

@@ -113,9 +113,11 @@ def hertz_corrected_viscelasticity_KVM(delta, _mod_constraint, Eu, Eapp,
     root = (contact_point - delta)
     pos = root > 0
     D = (1 - nu ** 2)
-    aa0 = 4 / 3 * np.sqrt(R) * (Eu - (Eu - Eapp) / (1 - np.exp(-0.365 * time_ind / lmd))) / D
+    aa0 = 4 / 3 * np.sqrt(R) * (Eu - (Eu - Eapp) /
+                                (1 - np.exp(-0.365 * time_ind / lmd))) / D
 
-    aa1 = 4 / 3 * np.sqrt(R) * ((Eu - Eapp) / (1 - np.exp(-0.365 * time_ind / lmd))) / D
+    aa1 = 4 / 3 * np.sqrt(R) * ((Eu - Eapp) /
+                                (1 - np.exp(-0.365 * time_ind / lmd))) / D
 
     bb = np.zeros_like(delta)
     bb[pos] = (root[pos]) ** (3 / 2)
@@ -253,7 +255,10 @@ def helper_retract_fraction(tip_position, force, fraction_force):
     fraction = int(0.8 * len(force))
     max_force = np.max(force[:fraction])
     max_position = -1 * tip_position[0]
-    fit_stop = int(np.argwhere(force < (max_force * fraction_force))[0])
+    _mask = force < (max_force * fraction_force)
+    if not np.any(_mask):
+        raise ValueError("No force values below threshold")
+    fit_stop = int(np.argmax(_mask))
     position_seg = tip_position[:fit_stop].copy()
     force_seg = force[:fit_stop].copy()
 

pyproject.toml

@@ -25,3 +25,8 @@ classifiers=[
     "Programming Language :: Python :: 3",
     "Intended Audience :: Science/Research"
 ]
+[dependency-groups]
+test = ["pytest"]
+dev = [
+    {include-group = "test"}
+]

tests/test_helper_retract_fraction.py

@@ -0,0 +1,50 @@
+import numpy as np
+import pytest
+
+
+def old_helper_retract_fraction(tip_position, force, fraction_force):
+    fraction = int(0.8 * len(force))
+    max_force = np.max(force[:fraction])
+    max_position = -1 * tip_position[0]
+    fit_stop = int(np.argwhere(force < (max_force * fraction_force))[0])
+    position_seg = tip_position[:fit_stop].copy()
+    force_seg = force[:fit_stop].copy()
+    return position_seg, force_seg, max_position, max_force
+
+
+def new_helper_retract_fraction(tip_position, force, fraction_force):
+    """Reference implementation using NumPy 1.x and 2.x compatible indexing."""
+    fraction = int(0.8 * len(force))
+    max_force = np.max(force[:fraction])
+    max_position = -1 * tip_position[0]
+
+    # new part
+    _mask = force < (max_force * fraction_force)
+    if not np.any(_mask):
+        raise ValueError("No force values below threshold")
+    fit_stop = int(np.argmax(_mask))
+
+    position_seg = tip_position[:fit_stop].copy()
+    force_seg = force[:fit_stop].copy()
+    return position_seg, force_seg, max_position, max_force
+
+
+def test_helper_behavior_differs_between_numpy_1_and_2():
+    # Deterministic example with a clear cutoff.
+    tip_position = np.linspace(0.0, -9.0, 10)
+    force = np.array([10.0, 9.0, 8.0, 6.0, 4.0, 3.0, 2.0, 1.0, 0.0, 0.0])
+    fraction_force = 0.5  # threshold = 5.0 -> first below threshold at index 4
+
+    major = int(np.__version__.split(".", 1)[0])
+    if major < 2:
+        old = old_helper_retract_fraction(tip_position, force, fraction_force)
+        new = new_helper_retract_fraction(tip_position, force, fraction_force)
+        np.testing.assert_allclose(old[0], new[0])
+        np.testing.assert_allclose(old[1], new[1])
+        assert old[2] == new[2]
+        assert old[3] == new[3]
+    else:
+        # Current implementation creates errors on NumPy 2.x
+        # ("only length-1 arrays..." / "only 0-d arrays...").
+        with pytest.raises(TypeError):
+            new_helper_retract_fraction(tip_position, force, fraction_force)