From 2aa7a632f0dd239af44d8c14ec34260ee239be6f Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 18 Mar 2026 08:53:46 +0000
Subject: [PATCH 1/3] Initial plan


From 32a4d99b0662e9b11653fad5876a98a73256828a Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 18 Mar 2026 09:32:14 +0000
Subject: [PATCH 2/3] Fix multiple code quality and correctness issues from
 reviewer feedback

Co-authored-by: Pablo1990 <1974224+Pablo1990@users.noreply.github.com>
---
 Tests/test_geo.py                             |  14 +-
 Tests/test_vertexModel.py                     |  15 +-
 Tests/tests.py                                |  20 +--
 src/pyVertexModel/__init__.py                 |   4 +-
 src/pyVertexModel/algorithm/integrators.py    |  38 ++---
 src/pyVertexModel/algorithm/vertexModel.py    |   9 +-
 .../analysis/analyse_simulation.py            |   1 -
 .../analysis/find_required_purse_string.py    |   9 +-
 src/pyVertexModel/geometry/geo.py             | 144 ++++++++++--------
 src/pyVertexModel/parameters/set.py           |  23 +++
 src/pyVertexModel/util/utils.py               |   2 +-
 11 files changed, 156 insertions(+), 123 deletions(-)

diff --git a/Tests/test_geo.py b/Tests/test_geo.py
index 899f286..46b065c 100644
--- a/Tests/test_geo.py
+++ b/Tests/test_geo.py
@@ -382,19 +382,19 @@ def test_geometry_is_correct(self):
         :return:
         """
         # Load data
-        vModel_test = load_data('geometry_correct_1.pkl')
+        _, _, vModel_test = load_data('geometry_correct_1.pkl')
 
         # Check if geometry is correct
         self.assertTrue(vModel_test.geo.geometry_is_correct())
 
         # Load data
-        vModel_test = load_data('geometry_correct_2.pkl')
+        _, _, vModel_test = load_data('geometry_correct_2.pkl')
 
         # Check if geometry is correct
         self.assertTrue(vModel_test.geo.geometry_is_correct())
 
         # Load data
-        vModel_test = load_data('geometry_correct_3.pkl')
+        _, _, vModel_test = load_data('geometry_correct_3.pkl')
 
         # Check if geometry is correct
         self.assertTrue(vModel_test.geo.geometry_is_correct())
@@ -405,25 +405,25 @@ def test_geometry_is_incorrect(self):
         :return:
         """
         # Load data
-        vModel_test = load_data('vertices_going_wild_1.pkl')
+        _, _, vModel_test = load_data('vertices_going_wild_1.pkl')
 
         # Check if geometry is correct
         self.assertFalse(vModel_test.geo.geometry_is_correct())
 
         # Another test with a different geometry
-        vModel_test = load_data('vertices_going_wild_2.pkl')
+        _, _, vModel_test = load_data('vertices_going_wild_2.pkl')
 
         # Check if geometry is correct
         self.assertFalse(vModel_test.geo.geometry_is_correct())
 
         # Another test with a different geometry
-        vModel_test = load_data('vertices_going_wild_3.pkl')
+        _, _, vModel_test = load_data('vertices_going_wild_3.pkl')
 
         # Check if geometry is correct
         self.assertFalse(vModel_test.geo.geometry_is_correct())
 
         # Another test with a different geometry
-        vModel_test = load_data('vertices_going_wild_4.pkl')
+        _, _, vModel_test = load_data('vertices_going_wild_4.pkl')
 
         # Check if geometry is correct
         self.assertFalse(vModel_test.geo.geometry_is_correct())
diff --git a/Tests/test_vertexModel.py b/Tests/test_vertexModel.py
index d72d438..a002bc9 100644
--- a/Tests/test_vertexModel.py
+++ b/Tests/test_vertexModel.py
@@ -1,3 +1,4 @@
+import os
 from os.path import exists
 
 import numpy as np
@@ -384,7 +385,7 @@ 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
+        test_dir = os.path.join(TEST_DIRECTORY, 'Tests_data', file_name)
         if exists(test_dir):
             vModel_test.set.initial_filename_state = test_dir
         else:
@@ -397,7 +398,7 @@ def test_initialize_voronoi_from_time_image(self):
 
         vModel_test.set = set_test
 
-        g_test, K_test, energies_test, _ = newtonRaphson.KgGlobal(vModel_test.geo, vModel_test.set, vModel_test.geo)
+        g_test, K_test, energies_test, _ = integrators.KgGlobal(vModel_test.geo, vModel_test.set, vModel_test.geo)
 
         # Check if energies are the same
         assert_array1D(g_test, mat_info['g'])
@@ -491,7 +492,7 @@ def test_weird_bug_should_not_happen(self):
         :return:
         """
         # Load data
-        vModel_test = load_data('vertices_going_wild_1.pkl')
+        _, _, vModel_test = load_data('vertices_going_wild_1.pkl')
 
         # Run for 20 iterations. dt should not decrease to 1e-1
         vModel_test.set.tend = vModel_test.t + 20 * vModel_test.set.dt0
@@ -524,9 +525,9 @@ def test_vertices_shouldnt_be_going_wild(self):
         :return:
         """
         # Load data
-        vModel_test = load_data('vertices_going_wild_2.pkl')
+        _, _, vModel_test = load_data('vertices_going_wild_2.pkl')
 
-        # Run for 10 iterations. dt should not decrease to 1e-1
+        # Run for 20 iterations. dt should not decrease to 1e-1
         vModel_test.set.tend = vModel_test.t + 20 * vModel_test.set.dt0
 
         # Update tolerance
@@ -544,9 +545,9 @@ def test_vertices_shouldnt_be_going_wild_3(self):
         :return:
         """
         # Load data
-        vModel_test = load_data('vertices_going_wild_3.pkl')
+        _, _, vModel_test = load_data('vertices_going_wild_3.pkl')
 
-        # Run for 10 iterations. dt should not decrease to 1e-1
+        # Run for 20 iterations. dt should not decrease to 1e-1
         vModel_test.set.tend = vModel_test.t + 20 * vModel_test.set.dt0
 
         # Update tolerance
diff --git a/Tests/tests.py b/Tests/tests.py
index f49b910..2604921 100644
--- a/Tests/tests.py
+++ b/Tests/tests.py
@@ -1,9 +1,11 @@
+import os
 import unittest
-from os.path import exists, abspath
+from os.path import exists
 
 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
@@ -11,12 +13,9 @@
 
 
 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():
@@ -26,7 +25,7 @@ def load_data(file_name, return_geo=True):
 
             if 'Set' in mat_info.keys():
                 set_test = Set(mat_info['Set'])
-                if set_test.OutputFolder.__eq__(b'') or set_test.OutputFolder is None:
+                if set_test.OutputFolder is None or set_test.OutputFolder == b'' or set_test.OutputFolder == '':
                     set_test.OutputFolder = '../Result/Test'
             else:
                 set_test = None
@@ -37,17 +36,14 @@ 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
         v_model.set.export_images = False
         v_model.set.VTK = False
 
-        return v_model
+        return None, None, v_model
     else:
         raise FileNotFoundError('File %s not found' % file_name)
 
diff --git a/src/pyVertexModel/__init__.py b/src/pyVertexModel/__init__.py
index b8cc1b0..3dd5d91 100644
--- a/src/pyVertexModel/__init__.py
+++ b/src/pyVertexModel/__init__.py
@@ -1,5 +1,6 @@
 import logging
 import os
+import sys
 import warnings
 from pathlib import Path
 
@@ -18,7 +19,8 @@
     datefmt="%Y/%m/%d %I:%M:%S %p",
 )
 # Check if logger already has a console handler to avoid adding multiple handlers
-if not any(isinstance(handler, logging.StreamHandler) for handler in logger.handlers):
+if not any(isinstance(h, logging.StreamHandler) and getattr(h, 'stream', None) in (sys.stdout, sys.stderr)
+           for h in logger.handlers):
     console_handler = logging.StreamHandler()
     console_handler.setFormatter(formatter)
     logger.addHandler(console_handler)
diff --git a/src/pyVertexModel/algorithm/integrators.py b/src/pyVertexModel/algorithm/integrators.py
index f38aafd..2d1b77e 100644
--- a/src/pyVertexModel/algorithm/integrators.py
+++ b/src/pyVertexModel/algorithm/integrators.py
@@ -73,7 +73,7 @@ def solve_remodeling_step(geo_0, geo_n, geo, dofs, c_set):
             f'Local Problem ->Iter: 0, ||gr||= {gr:.3e} ||dyr||= {dyr:.3e}  nu/nu0={c_set.nu / c_set.nu0:.3e}  '
             f'dt/dt0={c_set.dt / c_set.dt0:.3g}')
 
-        geo, g, k, energy, c_set, gr, dyr, dy, _ = newton_raphson(geo_0, geo_n, geo, dofs, c_set, k, g, -1, -1)
+        geo, g, k, energy, c_set, gr, dyr, dy = newton_raphson(geo_0, geo_n, geo, dofs, c_set, k, g, -1, -1)
 
         if gr > c_set.tol or dyr > c_set.tol or np.any(np.isnan(g[dofs.Free])) or np.any(np.isnan(dy[dofs.Free])):
             logger.info(f'Local Problem did not converge after {c_set.iter} iterations.')
@@ -591,7 +591,7 @@ def check_if_fire_converged(geo, f_flat, c_set, dy_flat, v_flat, iteration_count
 
     # Check maximum iterations first
     if iteration_count >= max_iter:
-        converged = True
+        converged = False
         reason = f"Reached max iterations ({max_iter})"
         logger.warning(f"FIRE: {reason} - maxF={max_force:.3e}")
         return converged, reason
@@ -638,12 +638,9 @@ def single_iteration_fire(geo, c_set, dof, dy, g, selected_cells=None):
 
     # Initialize FIRE state variables if not present
     if not hasattr(geo, '_fire_velocity'):
-        # Initialize velocity to zero
         initialize_fire(geo, c_set)
         logger.info("FIRE algorithm initialized")
 
-    geo._fire_velocity = np.zeros((geo.numF + geo.numY + geo.nCells, 3))
-
     # Increment iteration counter
     geo._fire_iteration_count += 1
 
@@ -654,8 +651,9 @@ def single_iteration_fire(geo, c_set, dof, dy, g, selected_cells=None):
     # Forces are negative gradient (F = -∇E = -g)
     F = -g.copy()
 
-    # Extract free DOF velocities and forces
-    v_flat = geo._fire_velocity.flatten()[dof]
+    # Extract free DOF velocities and forces using reshape for a writable view
+    velocity_view = geo._fire_velocity.reshape(-1)
+    v_flat = velocity_view[dof].copy()
     F_flat = F[dof]
 
     # ============================================
@@ -710,14 +708,14 @@ def single_iteration_fire(geo, c_set, dof, dy, g, selected_cells=None):
     # STEP 3: UPDATE GEOMETRY
     # ============================================
 
-    # Store updated velocity
-    geo._fire_velocity.flatten()[dof] = v_flat
+    # Store updated velocity back into geo._fire_velocity through the view
+    velocity_view[dof] = v_flat
 
     # Also update constrained DOFs if any
     if len(g_constrained) > 0:
         F_constrained = -g[g_constrained]
         dy_constrained = geo._fire_dt * F_constrained * (c_set.nu / c_set.nu_bottom)  # More conservative for constrained
-        dy.flatten()[g_constrained] = dy_constrained
+        dy.reshape(-1)[g_constrained] = dy_constrained
 
     # Build full displacement vector
     dy[dof, 0] = dy_flat
@@ -771,14 +769,11 @@ def fire_minimization_loop(geo, c_set, dof, g, t, num_step, selected_cells=None)
     """
     logger.info(f"Starting FIRE minimization for timestep t={t}")
     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)
+    # Initialize FIRE state variables (always re-initialize for a new minimization)
+    initialize_fire(geo, c_set)
+    geo._fire_velocity = np.zeros((geo.numF + geo.numY + geo.nCells, 3))
+    geo._fire_iteration_count = 0
 
     # Store initial gradient for reference
     initial_gradient_norm = np.linalg.norm(g[dof])
@@ -821,12 +816,9 @@ def fire_minimization_loop(geo, c_set, dof, g, t, num_step, selected_cells=None)
         logger.warning(f"FIRE minimization incomplete: {fire_iterations} iterations, "
                        f"gradient {final_gradient_norm:.3e}, force tolerance {c_set.fire_force_tol:.1e}")
 
-    # Reset FIRE state for next timestep (optional - comment out to maintain momentum)
-    if hasattr(geo, '_fire_velocity'):
-        geo._fire_velocity[:] = 0  # Reset velocities
-        geo._fire_dt = c_set.dt      # Reset timestep
-        geo._fire_alpha = c_set.fire_alpha_start
-        geo._fire_n_positive = 0
+    # Reset FIRE state for next timestep
+    initialize_fire(geo, c_set)
+    geo._fire_velocity[:] = 0
 
     return geo, converged, final_gradient_norm
 
diff --git a/src/pyVertexModel/algorithm/vertexModel.py b/src/pyVertexModel/algorithm/vertexModel.py
index b0c7790..4a1686f 100644
--- a/src/pyVertexModel/algorithm/vertexModel.py
+++ b/src/pyVertexModel/algorithm/vertexModel.py
@@ -541,6 +541,8 @@ def single_iteration(self, post_operations=True):
                                                                          self.numStep)
             if converged:
                 self.iteration_converged()
+            else:
+                self.iteration_did_not_converged()
         else:
             self.geo, g, __, __, self.set, gr, dyr, dy = integrators.newton_raphson(self.geo_0, self.geo_n, self.geo,
                                                                                       self.Dofs, self.set, K, g,
@@ -591,7 +593,7 @@ def iteration_did_not_converged(self):
                 self.set.nu = 10 * self.set.nu0
             else:
                 if (self.set.iter >= self.set.MaxIter and
-                        (self.set.dt / self.set.dt0) > self.set.dt_tolerance):
+                        (self.set.dt / self.set.dt0) > getattr(self.set, 'dt_tolerance', 1e-6)):
                     self.set.MaxIter = self.set.MaxIter0
                     self.set.nu = self.set.nu0
                     self.set.dt = self.set.dt / 2
@@ -599,7 +601,10 @@ def iteration_did_not_converged(self):
                 else:
                     self.didNotConverge = True
         else:
-            logger.warning("FIRE did not converge")
+            # FIRE did not converge within its iteration budget: restore geometry and mark failure
+            logger.warning("FIRE minimization did not converge; restoring backup geometry")
+            self.geo, self.geo_n, self.geo_0, self.tr, self.Dofs = load_backup_vars(self.backupVars)
+            self.didNotConverge = True
 
     def iteration_converged(self):
         """
diff --git a/src/pyVertexModel/analysis/analyse_simulation.py b/src/pyVertexModel/analysis/analyse_simulation.py
index 1dd2697..9c54d24 100644
--- a/src/pyVertexModel/analysis/analyse_simulation.py
+++ b/src/pyVertexModel/analysis/analyse_simulation.py
@@ -417,7 +417,6 @@ def analyse_edge_recoil(file_name_v_model, type_of_ablation='recoil_edge_info_ap
                 v_model.set.RemodelingFrequency = 100
             v_model.set.ablation = False
             v_model.set.export_images = False
-            v_model.set.integrator = 'euler'
             v_model.set.purseStringStrength = 0
             v_model.set.lateralCablesStrength = 0
             if v_model.set.export_images and not os.path.exists(v_model.set.OutputFolder + '/images'):
diff --git a/src/pyVertexModel/analysis/find_required_purse_string.py b/src/pyVertexModel/analysis/find_required_purse_string.py
index abcb5a9..b3eb765 100644
--- a/src/pyVertexModel/analysis/find_required_purse_string.py
+++ b/src/pyVertexModel/analysis/find_required_purse_string.py
@@ -135,7 +135,7 @@
 
         # Get directory within directory
         dirs_within = os.listdir(os.path.join(c_folder, ar_dir, directory))
-        dirs_within = [os.path.join(c_folder, ar_dir, directory, d) for d in dirs_within if d.startswith('no_Remodelling_ablating_')]
+        dirs_within = sorted([os.path.join(c_folder, ar_dir, directory, d) for d in dirs_within if d.startswith('no_Remodelling_ablating_')])
         if len(dirs_within) == 0:
             print(f"No directories starting with 'no_Remodelling_ablating_' found in {directory}, skipping.")
             continue
@@ -153,8 +153,11 @@
 
         load_state(vModel, os.path.join(c_folder, ar_dir, directory, 'before_ablation.pkl'))
         t_ablation = vModel.t
-        vModel.set.integrator = 'euler'
-        vModel.set.dt_tolerance = 1e-1
+        # Only set integrator/tolerance if the model doesn't already have them from the saved state
+        if not getattr(vModel.set, 'integrator', None):
+            vModel.set.integrator = 'euler'
+        if not getattr(vModel.set, 'dt_tolerance', None):
+            vModel.set.dt_tolerance = 1e-1
 
         # Run the required purse string strength analysis
         current_folder = vModel.set.OutputFolder
diff --git a/src/pyVertexModel/geometry/geo.py b/src/pyVertexModel/geometry/geo.py
index 815ecd6..213ccae 100644
--- a/src/pyVertexModel/geometry/geo.py
+++ b/src/pyVertexModel/geometry/geo.py
@@ -2300,85 +2300,97 @@ def geometry_is_correct(self):
             return False
         
         # 1. Check cell Y coordinates
-        for cell_idx, cell in enumerate(alive_cells):
-            if cell.Y is None:
-                logger.debug(f"Cell {cell_idx} has Y=None")
+        for c_cell in alive_cells:
+            if c_cell.Y is None:
+                logger.debug(f"Cell {c_cell.ID} has Y=None")
                 return False
-            
+
             # Check for NaN or Inf
-            if np.any(np.isnan(cell.Y)) or np.any(np.isinf(cell.Y)):
-                logger.debug(f"Cell {cell_idx} has NaN or Inf in Y coordinates")
+            if np.any(np.isnan(c_cell.Y)) or np.any(np.isinf(c_cell.Y)):
+                logger.debug(f"Cell {c_cell.ID} has NaN or Inf in Y coordinates")
                 return False
         
         # 2. Check cell volumes (should be positive and reasonable)
-        for cell_idx, cell in enumerate(alive_cells):
-            if not hasattr(cell, 'Vol') or cell.Vol is None:
-                logger.debug(f"Cell {cell_idx} has no volume")
+        for c_cell in alive_cells:
+            if not hasattr(c_cell, 'Vol') or c_cell.Vol is None:
+                logger.debug(f"Cell {c_cell.ID} has no volume")
                 return False
-            
-            if cell.Vol <= 0:
-                logger.debug(f"Cell {cell_idx} has non-positive volume: {cell.Vol}")
+
+            if not np.isfinite(c_cell.Vol):
+                logger.debug(f"Cell {c_cell.ID} has non-finite volume: {c_cell.Vol}")
                 return False
-            
-            if cell.Vol < 1e-10:
-                logger.debug(f"Cell {cell_idx} has tiny volume: {cell.Vol}")
+
+            if c_cell.Vol <= 0:
+                logger.debug(f"Cell {c_cell.ID} has non-positive volume: {c_cell.Vol}")
                 return False
-        
+
+            if c_cell.Vol < 1e-10:
+                logger.debug(f"Cell {c_cell.ID} has tiny volume: {c_cell.Vol}")
+                return False
+
         # 3. Check face areas (should be positive and not too small)
-        for cell_idx, cell in enumerate(alive_cells):
-            if not hasattr(cell, 'Faces'):
+        for c_cell in alive_cells:
+            if not hasattr(c_cell, 'Faces'):
                 continue
-                
-            for face_idx, face in enumerate(cell.Faces):
-                if hasattr(face, 'Area') and face.Area is not None:
-                    if face.Area < 0:
-                        logger.debug(f"Cell {cell_idx}, Face {face_idx} has negative area: {face.Area}")
+
+            for c_face in c_cell.Faces:
+                if hasattr(c_face, 'Area') and c_face.Area is not None:
+                    if not np.isfinite(c_face.Area):
+                        logger.debug(f"Cell {c_cell.ID}, Face has non-finite area: {c_face.Area}")
+                        return False
+
+                    if c_face.Area < 0:
+                        logger.debug(f"Cell {c_cell.ID}, Face has negative area: {c_face.Area}")
                         return False
-                    
-                    if face.Area < 1e-12:
-                        logger.debug(f"Cell {cell_idx}, Face {face_idx} has degenerate area: {face.Area}")
+
+                    if c_face.Area < 1e-12:
+                        logger.debug(f"Cell {c_cell.ID}, Face has degenerate area: {c_face.Area}")
                         return False
-        
+
         # 4. Check triangle areas and degeneracy
-        for cell_idx, cell in enumerate(alive_cells):
-            if not hasattr(cell, 'Faces'):
+        for c_cell in alive_cells:
+            if not hasattr(c_cell, 'Faces'):
                 continue
-                
-            for face_idx, face in enumerate(cell.Faces):
-                if not hasattr(face, 'Tris'):
+
+            for c_face in c_cell.Faces:
+                if not hasattr(c_face, 'Tris'):
                     continue
-                    
-                for tri_idx, tri in enumerate(face.Tris):
+
+                for tri_idx, tri in enumerate(c_face.Tris):
                     # Check for degenerate triangles (same edge vertices)
                     if hasattr(tri, 'Edge') and tri.Edge is not None:
                         if len(tri.Edge) >= 2 and tri.Edge[0] == tri.Edge[1]:
-                            logger.debug(f"Cell {cell_idx}, Face {face_idx}, Tri {tri_idx} is degenerate (same vertices)")
+                            logger.debug(f"Cell {c_cell.ID}, Tri {tri_idx} is degenerate (same vertices)")
                             return False
-                    
+
                     # Check triangle area
                     if hasattr(tri, 'Area') and tri.Area is not None:
+                        if not np.isfinite(tri.Area):
+                            logger.debug(f"Cell {c_cell.ID}, Tri {tri_idx} has non-finite area")
+                            return False
+
                         if tri.Area < 0:
-                            logger.debug(f"Cell {cell_idx}, Face {face_idx}, Tri {tri_idx} has negative area")
+                            logger.debug(f"Cell {c_cell.ID}, Tri {tri_idx} has negative area")
                             return False
-        
+
         # 4b. Check for excessive degenerate/tiny triangles
         # Some tiny triangles are normal, but too many indicate geometric issues
         tiny_triangle_count = 0
         total_triangles = 0
-        
-        for cell_idx, cell in enumerate(alive_cells):
-            if not hasattr(cell, 'Faces'):
+
+        for c_cell in alive_cells:
+            if not hasattr(c_cell, 'Faces'):
                 continue
-                
-            for face_idx, face in enumerate(cell.Faces):
-                if not hasattr(face, 'Tris'):
+
+            for c_face in c_cell.Faces:
+                if not hasattr(c_face, 'Tris'):
                     continue
-                    
-                for tri in face.Tris:
+
+                for tri in c_face.Tris:
                     total_triangles += 1
-                    
-                    # Count tiny triangles
-                    if hasattr(tri, 'Area') and tri.Area is not None and tri.Area < 1e-10:
+
+                    # Count tiny triangles (skip non-finite values)
+                    if hasattr(tri, 'Area') and tri.Area is not None and np.isfinite(tri.Area) and tri.Area < 1e-10:
                         tiny_triangle_count += 1
         
         # If more than 0.2% of triangles are tiny, geometry is problematic
@@ -2394,47 +2406,47 @@ def geometry_is_correct(self):
         # Non-planar faces indicate vertices sticking out, creating spikes.
         planarity_threshold = 0.08  # Eigenvalue ratio threshold
         non_planar_count = 0
-        
-        for cell_idx, cell in enumerate(alive_cells):
-            if not hasattr(cell, 'Faces'):
+
+        for c_cell in alive_cells:
+            if not hasattr(c_cell, 'Faces'):
                 continue
-                
-            for face_idx, face in enumerate(cell.Faces):
-                if not hasattr(face, 'Tris'):
+
+            for c_face in c_cell.Faces:
+                if not hasattr(c_face, 'Tris'):
                     continue
-                    
+
                 # Get all unique vertices in the face
                 vertices = set()
-                for tri in face.Tris:
+                for tri in c_face.Tris:
                     if hasattr(tri, 'Edge') and tri.Edge is not None:
                         vertices.add(tri.Edge[0])
                         vertices.add(tri.Edge[1])
-                
+
                 if len(vertices) < 3:
                     continue
-                
+
                 # Get vertex coordinates
                 try:
-                    vertex_coords = np.array([cell.Y[v] for v in vertices])
-                    
+                    vertex_coords = np.array([c_cell.Y[v] for v in vertices])
+
                     # Check planarity using PCA
                     # For planar faces, the smallest eigenvalue should be very small
                     if len(vertex_coords) >= 3:
                         centered = vertex_coords - np.mean(vertex_coords, axis=0)
                         cov = np.cov(centered.T)
                         eigenvalues = np.linalg.eigvalsh(cov)
-                        
+
                         # Smallest eigenvalue relative to largest
                         # This detects non-planar faces that create "spiky cells"
                         if eigenvalues[2] > 1e-15:  # Avoid division by zero
                             planarity_ratio = eigenvalues[0] / eigenvalues[2]
-                            
+
                             if planarity_ratio > planarity_threshold:
                                 non_planar_count += 1
-                                logger.debug(f"Cell {cell_idx}, Face {face_idx} is non-planar (ratio={planarity_ratio:.4f} > {planarity_threshold})")
-                
+                                logger.debug(f"Cell {c_cell.ID}, Face is non-planar (ratio={planarity_ratio:.4f} > {planarity_threshold})")
+
                 except (IndexError, ValueError) as e:
-                    logger.debug(f"Cell {cell_idx}, Face {face_idx}: Error checking planarity: {e}")
+                    logger.debug(f"Cell {c_cell.ID}: Error checking planarity: {e}")
                     return False
         
         # Any non-planar faces indicate geometry is going wild
diff --git a/src/pyVertexModel/parameters/set.py b/src/pyVertexModel/parameters/set.py
index bda8ad1..8dd871a 100644
--- a/src/pyVertexModel/parameters/set.py
+++ b/src/pyVertexModel/parameters/set.py
@@ -187,6 +187,29 @@ def __init__(self, mat_file=None):
         else:
             self.read_mat_file(mat_file)
 
+    def __setstate__(self, state):
+        """Ensure backward compatibility when loading older pickle files."""
+        self.__dict__.update(state)
+        # New attributes not present in older pickled instances get their defaults
+        defaults = {
+            'integrator': 'euler',
+            'dt_tolerance': 1e-6,
+            'fire_max_iterations': 100,
+            'fire_alpha_start': 0.1,
+            '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_f_alpha': 0.99,
+            'fire_force_tol': 1e-6,
+            'fire_disp_tol': 1e-8,
+            'fire_vel_tol': 1e-10,
+        }
+        for key, value in defaults.items():
+            if not hasattr(self, key):
+                setattr(self, key, value)
+
     def check_for_non_used_parameters(self):
         """
         Adjust configuration attributes based on feature toggles, disabling unused options and setting related defaults.
diff --git a/src/pyVertexModel/util/utils.py b/src/pyVertexModel/util/utils.py
index 853100b..a1f0471 100644
--- a/src/pyVertexModel/util/utils.py
+++ b/src/pyVertexModel/util/utils.py
@@ -694,7 +694,7 @@ def plot_figure_with_line(best_average_values, scutoids, current_path, x_axis_na
     # From 0 ylim always
     plt.ylim(0, None)
 
-    if y_axis_label == 'Purse string strength (t=' + str(0.1) + ')':
+    if 'Purse string strength' in y_axis_label:
         plt.ylim(0, 1.2e-3)
 
     if y_axis_name == 'top_closure_velocity':

From 73944a5a2b7cb9ac028d68e7e813001e9c07bda0 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 18 Mar 2026 09:35:40 +0000
Subject: [PATCH 3/3] Address code review: use class-level constant for
 backward compat defaults, use hasattr() instead of getattr()

Co-authored-by: Pablo1990 <1974224+Pablo1990@users.noreply.github.com>
---
 .../analysis/find_required_purse_string.py    |  4 +--
 src/pyVertexModel/parameters/set.py           | 34 ++++++++++---------
 2 files changed, 20 insertions(+), 18 deletions(-)

diff --git a/src/pyVertexModel/analysis/find_required_purse_string.py b/src/pyVertexModel/analysis/find_required_purse_string.py
index b3eb765..0652f30 100644
--- a/src/pyVertexModel/analysis/find_required_purse_string.py
+++ b/src/pyVertexModel/analysis/find_required_purse_string.py
@@ -154,9 +154,9 @@
         load_state(vModel, os.path.join(c_folder, ar_dir, directory, 'before_ablation.pkl'))
         t_ablation = vModel.t
         # Only set integrator/tolerance if the model doesn't already have them from the saved state
-        if not getattr(vModel.set, 'integrator', None):
+        if not hasattr(vModel.set, 'integrator'):
             vModel.set.integrator = 'euler'
-        if not getattr(vModel.set, 'dt_tolerance', None):
+        if not hasattr(vModel.set, 'dt_tolerance'):
             vModel.set.dt_tolerance = 1e-1
 
         # Run the required purse string strength analysis
diff --git a/src/pyVertexModel/parameters/set.py b/src/pyVertexModel/parameters/set.py
index 8dd871a..67c5b86 100644
--- a/src/pyVertexModel/parameters/set.py
+++ b/src/pyVertexModel/parameters/set.py
@@ -12,6 +12,23 @@
 logger = logging.getLogger("pyVertexModel")
 
 class Set:
+    # Defaults for attributes added after initial release, used by __setstate__ for backward compatibility
+    _BACKWARD_COMPAT_DEFAULTS = {
+        'integrator': 'euler',
+        'dt_tolerance': 1e-6,
+        'fire_max_iterations': 100,
+        'fire_alpha_start': 0.1,
+        '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_f_alpha': 0.99,
+        'fire_force_tol': 1e-6,
+        'fire_disp_tol': 1e-8,
+        'fire_vel_tol': 1e-10,
+    }
+
     def __init__(self, mat_file=None):
         """
         Initialize simulation configuration attributes with sensible defaults.
@@ -191,22 +208,7 @@ def __setstate__(self, state):
         """Ensure backward compatibility when loading older pickle files."""
         self.__dict__.update(state)
         # New attributes not present in older pickled instances get their defaults
-        defaults = {
-            'integrator': 'euler',
-            'dt_tolerance': 1e-6,
-            'fire_max_iterations': 100,
-            'fire_alpha_start': 0.1,
-            '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_f_alpha': 0.99,
-            'fire_force_tol': 1e-6,
-            'fire_disp_tol': 1e-8,
-            'fire_vel_tol': 1e-10,
-        }
-        for key, value in defaults.items():
+        for key, value in self._BACKWARD_COMPAT_DEFAULTS.items():
             if not hasattr(self, key):
                 setattr(self, key, value)