Fix test suite failures: NumPy 2.0 compatibility, API regressions, and FIRE algorithm initialization

.github/workflows/tests.yaml

@@ -4,6 +4,8 @@ on:
   workflow_dispatch:
   push:
     branches: [ main ]
+  pull_request:
+    branches: [ main ]
 
 jobs:
   build:

Tests/test_vertexModel.py

@@ -1,3 +1,4 @@
+import os
 from os.path import exists
 
 import numpy as np
@@ -189,21 +190,21 @@ def test_iteration_did_not_converged(self):
         check_if_cells_are_the_same(geo_original, v_model_test.geo)
 
         # Save backup variables
-        geo_test.Cells[0].Y[0, 0] = np.Inf
-        geo_test.Cells[0].Faces[0].Centre[0] = np.Inf
+        geo_test.Cells[0].Y[0, 0] = np.inf
+        geo_test.Cells[0].Faces[0].Centre[0] = np.inf
         v_model_test.backupVars = save_backup_vars(geo_test, geo_test, geo_test, 0, DegreesOfFreedom(mat_info['Dofs']))
         v_model_test.set.iter = 1000000
         v_model_test.set.MaxIter0 = v_model_test.set.iter
         v_model_test.set.last_t_converged = 0.5
 
         v_model_test.iteration_did_not_converged()
 
-        geo_test.Cells[0].Y[0, 0] = -np.Inf
-        geo_test.Cells[0].Faces[0].Centre[0] = -np.Inf
+        geo_test.Cells[0].Y[0, 0] = -np.inf
+        geo_test.Cells[0].Faces[0].Centre[0] = -np.inf
 
         # Check if the cells are initialized correctly
-        np.testing.assert_equal(v_model_test.geo.Cells[0].Y[0, 0], np.Inf)
-        np.testing.assert_equal(v_model_test.geo.Cells[0].Faces[0].Centre[0], np.Inf)
+        np.testing.assert_equal(v_model_test.geo.Cells[0].Y[0, 0], np.inf)
+        np.testing.assert_equal(v_model_test.geo.Cells[0].Faces[0].Centre[0], np.inf)
 
     def test_newton_raphson_cyst(self):
         """
@@ -251,11 +252,8 @@ def test_obtain_initial_x_and_tetrahedra(self):
         vModel_test = VertexModelVoronoiFromTimeImage(set_test)
 
         file_name = 'LblImg_imageSequence.mat'
-        test_dir = 'Tests/Tests_data/%s' % file_name
-        if exists(test_dir):
-            Twg_test, X_test = vModel_test.obtain_initial_x_and_tetrahedra(test_dir)
-        else:
-            Twg_test, X_test = vModel_test.obtain_initial_x_and_tetrahedra('Tests_data/%s' % file_name)
+        test_dir = os.path.join(TEST_DIRECTORY, 'Tests_data', file_name)
+        Twg_test, X_test = vModel_test.obtain_initial_x_and_tetrahedra(test_dir)
 
         # Check if the test and expected are the same
         assert_matrix(Twg_test, mat_info_expected['Twg'] - 1)
@@ -361,11 +359,8 @@ def test_process_image(self):
         """
         # Process image
         file_name = 'LblImg_imageSequence.mat'
-        test_dir = 'Tests/Tests_data/%s' % file_name
-        if exists(test_dir):
-            _, imgStackLabelled_test = process_image(test_dir)
-        else:
-            _, imgStackLabelled_test = process_image('Tests_data/%s' % file_name)
+        test_dir = os.path.join(TEST_DIRECTORY, 'Tests_data', file_name)
+        _, imgStackLabelled_test = process_image(test_dir)
 
         # Load expected
         _, _, mat_info_expected = load_data('process_image_wingdisc_expected.mat')
@@ -384,11 +379,8 @@ def test_initialize_voronoi_from_time_image(self):
         # Test if initialize geometry function does not change anything
         vModel_test = VertexModelVoronoiFromTimeImage(set_test)
         file_name = 'voronoi_40cells.pkl'
-        test_dir = TEST_DIRECTORY + '/Tests/Tests_data/%s' % file_name
-        if exists(test_dir):
-            vModel_test.set.initial_filename_state = test_dir
-        else:
-            vModel_test.set.initial_filename_state = 'Tests_data/%s' % file_name
+        test_dir = os.path.join(TEST_DIRECTORY, 'Tests_data', file_name)
+        vModel_test.set.initial_filename_state = test_dir
 
         vModel_test.initialize()
 

Tests/tests.py

@@ -1,26 +1,27 @@
+import os
 import unittest
 from os.path import exists, abspath
 
 import numpy as np
 import scipy.io
 
+from Tests import TEST_DIRECTORY
 from pyVertexModel.algorithm.vertexModelVoronoiFromTimeImage import VertexModelVoronoiFromTimeImage
 from pyVertexModel.geometry.geo import Geo
 from pyVertexModel.parameters.set import Set
 from pyVertexModel.util.utils import load_state
 
 
 def load_data(file_name, return_geo=True):
-    test_dir = abspath('Tests/Tests_data/%s' % file_name)
+    test_dir = os.path.join(TEST_DIRECTORY, 'Tests_data', file_name)
     if file_name.endswith('.mat'):
-        if exists(test_dir):
-            mat_info = scipy.io.loadmat(test_dir)
-        else:
-            mat_info = scipy.io.loadmat('Tests_data/%s' % file_name)
+        mat_info = scipy.io.loadmat(test_dir)
 
         if return_geo:
             if 'Geo' in mat_info.keys():
                 geo_test = Geo(mat_info['Geo'])
+                geo_test.update_lmin0()
+                geo_test.update_barrier_tri0()
             else:
                 geo_test = None
 
@@ -37,10 +38,7 @@ def load_data(file_name, return_geo=True):
         return geo_test, set_test, mat_info
     elif file_name.endswith('.pkl'):
         v_model = VertexModelVoronoiFromTimeImage(create_output_folder=False)
-        if exists(test_dir):
-            load_state(v_model, test_dir)
-        else:
-            load_state(v_model, 'Tests_data/%s' % file_name)
+        load_state(v_model, test_dir)
 
         # Set parameters to avoid file output during tests
         v_model.set.OutputFolder = None

src/pyVertexModel/algorithm/integrators.py

@@ -773,12 +773,8 @@ def fire_minimization_loop(geo, c_set, dof, g, t, num_step, selected_cells=None)
     dy = np.zeros(((geo.numY + geo.numF + geo.nCells) * 3, 1), dtype=np.float64)
     geo._fire_velocity = np.zeros((geo.numF + geo.numY + geo.nCells, 3))
 
-    # Reset FIRE iteration counter for this minimization
-    if hasattr(geo, '_fire_velocity'):
-        geo._fire_iteration_count = 0
-    else:
-        # Initialize if not already
-        initialize_fire(geo, c_set)
+    # Always (re-)initialize FIRE state for this minimization step
+    initialize_fire(geo, c_set)
 
     # Store initial gradient for reference
     initial_gradient_norm = np.linalg.norm(g[dof])

src/pyVertexModel/algorithm/vertexModel.py

@@ -248,6 +248,11 @@ def __init__(self, set_option='wing_disc', c_set=None, create_output_folder=True
         self.geo = None
 
         # Set definition
+        # Allow a Set object to be passed as the first positional argument for backward compatibility
+        if isinstance(set_option, Set):
+            c_set = set_option
+            set_option = 'wing_disc'
+
         if c_set is not None:
             self.set = c_set
         else:
@@ -288,6 +293,14 @@ def initialize(self, img_input=None):
         Raises:
             FileNotFoundError: If the configured initial filename does not exist and img_input is None.
         """
+        if self.set.initial_filename_state is None:
+            # No saved state file configured: go straight to initialize_cells
+            if img_input is None:
+                self.initialize_cells(None)
+            else:
+                self.initialize_cells(img_input)
+            return
+
         filename = os.path.join(PROJECT_DIRECTORY, self.set.initial_filename_state)
 
         if not os.path.exists(filename) and img_input is None:

src/pyVertexModel/algorithm/vertexModelBubbles.py

@@ -7,6 +7,7 @@
 from scipy.spatial import Delaunay
 
 from pyVertexModel.algorithm.vertexModel import VertexModel
+from pyVertexModel.geometry.geo import Geo
 from pyVertexModel.util.utils import save_state
 
 logger = logging.getLogger("pyVertexModel")
@@ -387,8 +388,15 @@ def initialize_cells(self, filename):
         Initialize the geometry and the topology of the model.
         :return:
         """
-        # Build nodal mesh
-        self.generate_Xs(self.geo.nx, self.geo.ny, self.geo.nz)
+        # Ensure geo is initialized before generating X positions
+        if self.geo is None:
+            self.geo = Geo()
+
+        # Build nodal mesh (geo.nx/ny/nz may be None for non-Bubbles inputs)
+        nx = getattr(self.geo, 'nx', None)
+        ny = getattr(self.geo, 'ny', None)
+        nz = getattr(self.geo, 'nz', None)
+        self.generate_Xs(nx, ny, nz)
 
         # This code is to match matlab's output and python's
         # N = 3  # The dimensions of our points
@@ -414,23 +422,26 @@ def initialize_cells(self, filename):
             self.geo.XgBottom = self.geo.XgID[Xg[:, 2] < np.mean(self.X[:, 2])]
             self.geo.XgTop = self.geo.XgID[Xg[:, 2] > np.mean(self.X[:, 2])]
 
-        self.geo.Main_cells = range(len(self.geo.nCells))
+        self.geo.Main_cells = range(self.geo.nCells)
         self.geo.build_cells(self.set, self.X, Twg)
 
         if self.set.InputGeo == 'Bubbles_Cyst':
             # Extrapolate Face centres and Ys to the ellipsoid
             self.geo = extrapolate_ys_faces_ellipsoid(self.geo, self.set)
 
         # Save state with filename using the number of cells
-        filename = filename.replace('.tif', f'_{self.set.TotalCells}cells.pkl')
-        save_state(self.geo, filename)
+        if filename is not None:
+            filename = filename.replace('.tif', f'_{self.set.TotalCells}cells.pkl')
+            save_state(self.geo, filename)
 
     def generate_Xs(self, nx=None, ny=None, nz=None):
         """
         Generate the nodal positions of the mesh based on the input geometry
         :return:
         """
         self.X, X_IDs = build_topo(self.set, nx, ny, nz)
+        if self.geo is None:
+            self.geo = Geo()
         self.geo.nCells = self.X.shape[0]
         # Centre Nodal position at (0,0)
         self.X[:, 0] = self.X[:, 0] - np.mean(self.X[:, 0])

src/pyVertexModel/geometry/face.py

@@ -111,7 +111,7 @@ def build_face(self, ci, cj, face_ids, nCells, Cell, XgID, Set, XgTop, XgBottom,
         if oldFace is not None and getattr(oldFace, 'ij', None) is not None:
             self.Area0 = oldFace.Area0
         else:
-            self.Area0 = self.Area * Set.ref_A0
+            self.Area0 = self.Area * (Set.ref_A0 if Set.ref_A0 is not None else 1.0)
 
 
     def build_interface_type(self, ij, XgID, XgTop, XgBottom):

src/pyVertexModel/parameters/set.py

@@ -12,14 +12,15 @@
 logger = logging.getLogger("pyVertexModel")
 
 class Set:
-    def __init__(self, mat_file=None):
+    def __init__(self, mat_file=None, set_option=None):
         """
         Initialize simulation configuration attributes with sensible defaults.
 
         When mat_file is None, populate a comprehensive set of simulation parameters (topology, geometry, mechanics, time stepping, substrate, viscosity, remodeling, solver, boundary/loading, postprocessing, and ablation defaults). When mat_file is provided, load parameter values from the given MATLAB-like structure via read_mat_file(mat_file).
 
         Parameters:
             mat_file (optional): A MATLAB-style struct or object containing saved Set parameters; when provided, values are read and assigned to this instance.
+            set_option (str, optional): Name of a preset configuration method to call after initialization (e.g. 'wing_disc', 'cyst'). Ignored if mat_file is provided.
         """
         self.dt_tolerance = 1e-6
         self.min_3d_neighbours = None
@@ -184,6 +185,9 @@ def __init__(self, mat_file=None):
             self.VTK_iter = False
             self.SaveWorkspace = False
             self.SaveSetting = False
+
+            if set_option is not None and hasattr(self, set_option):
+                getattr(self, set_option)()
         else:
             self.read_mat_file(mat_file)
 
@@ -259,7 +263,12 @@ def read_mat_file(self, mat_file):
             if len(mat_file[param][0][0]) == 0:
                 setattr(self, param, None)
             else:
-                setattr(self, param, mat_file[param][0][0][0][0])
+                value = mat_file[param][0][0][0][0]
+                # Convert numpy integer/float scalars to Python native types to avoid
+                # integer overflow when used in arithmetic (e.g., uint8 * 600).
+                if hasattr(value, 'item'):
+                    value = value.item()
+                setattr(self, param, value)
 
     def define_if_not_defined(self, param, value):
         """
@@ -294,6 +303,7 @@ def update_derived_parameters(self):
         self.define_if_not_defined("MaxIter0", self.MaxIter)
         self.define_if_not_defined("contributionOldFaceCentre", self.contributionOldYs)
         self.define_if_not_defined("nu_bottom", self.nu * 600)
+        self.define_if_not_defined("ref_A0", 0.99)
 
         current_datetime = datetime.now()
         new_outputFolder = ''.join([PROJECT_DIRECTORY, '/Result/', str(current_datetime.strftime("%m-%d_%H%M%S_")),
@@ -507,4 +517,34 @@ def copy(self):
 
         copy_non_mutable_attributes(self, '', set_copy)
 
-        return set_copy
+        return set_copy
+
+    # Default values for attributes added after initial release, to maintain
+    # backward compatibility when loading old .pkl files.
+    _DEFAULTS = {
+        'integrator': 'euler',
+        'dt_tolerance': 1e-6,
+        'fire_dt_max': 0.5,
+        'fire_dt_min': 1e-8,
+        'fire_N_min': 5,
+        'fire_f_inc': 1.1,
+        'fire_f_dec': 0.5,
+        'fire_alpha_start': 0.1,
+        'fire_f_alpha': 0.99,
+        'fire_force_tol': 1e-6,
+        'fire_disp_tol': 1e-8,
+        'fire_vel_tol': 1e-10,
+        'fire_max_iterations': 100,
+        # Defaults for attributes not stored in legacy .mat files:
+        'noise_random': 0,
+        'model_name': '',
+        'percentage_scutoids': 0,
+        'purseStringStrength': 0,
+        'lateralCablesStrength': 0,
+        'export_images': True,
+    }
+
+    def __getattr__(self, name):
+        if name in Set._DEFAULTS:
+            return Set._DEFAULTS[name]
+        raise AttributeError(f"'Set' object has no attribute '{name}'")

src/pyVertexModel/util/utils.py

@@ -338,8 +338,8 @@ def ismember_rows(a, b):
         b = np.sort(b, axis=1)
         void_b = np.ascontiguousarray(b).view(np.dtype((np.bytes_, b.dtype.itemsize * b.shape[1])))
 
-    # Using numpy's in1d method for finding the presence of 'a' rows in 'b'
-    bool_array = np.in1d(void_a, void_b)
+    # Using numpy's isin method for finding the presence of 'a' rows in 'b'
+    bool_array = np.isin(void_a, void_b)
 
     # Finding the indices where the rows of 'a' are found in 'b'
     index_array = np.array([np.where(void_b == row)[0][0] if row in void_b else -1 for row in void_a])