diff --git a/CHANGELOG b/CHANGELOG index 179ac4c..54569b1 100644 --- a/CHANGELOG +++ b/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 diff --git a/nanite_community/models/model_hertz_corrected_viscoelasticity_KVM.py b/nanite_community/models/model_hertz_corrected_viscoelasticity_KVM.py index 53012c8..be03ffd 100644 --- a/nanite_community/models/model_hertz_corrected_viscoelasticity_KVM.py +++ b/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() diff --git a/pyproject.toml b/pyproject.toml index d4f2b56..dfdd5ea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,3 +25,8 @@ classifiers=[ "Programming Language :: Python :: 3", "Intended Audience :: Science/Research" ] +[dependency-groups] +test = ["pytest"] +dev = [ + {include-group = "test"} +] diff --git a/tests/test_helper_retract_fraction.py b/tests/test_helper_retract_fraction.py new file mode 100644 index 0000000..c167aca --- /dev/null +++ b/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)