Fix FIRE integrator momentum, convergence, and backward-compatibility issues

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())

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

Tests/tests.py

@@ -1,22 +1,21 @@
+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
 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():
@@ -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)
 

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)

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
 

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,15 +593,18 @@ 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
                     self.t = self.set.last_t_converged + self.set.dt
                 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):
         """

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'):

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 hasattr(vModel.set, 'integrator'):
+            vModel.set.integrator = 'euler'
+        if not hasattr(vModel.set, 'dt_tolerance'):
+            vModel.set.dt_tolerance = 1e-1
 
         # Run the required purse string strength analysis
         current_folder = vModel.set.OutputFolder