diff --git a/Tests/__init__.py b/Tests/__init__.py index 9e94e63c..956e0a89 100644 --- a/Tests/__init__.py +++ b/Tests/__init__.py @@ -1,3 +1,29 @@ +import logging import os +import warnings -TEST_DIRECTORY = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) \ No newline at end of file +TEST_DIRECTORY = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + +PROJECT_DIRECTORY = os.getenv('PROJECT_DIR', os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +# get the logger instance +logger = logging.getLogger("pyVertexModel") + +formatter = logging.Formatter( + "%(levelname)s [%(asctime)s] pyVertexModel: %(message)s", + datefmt="%Y/%m/%d %I:%M:%S %p", +) +console_handler = logging.StreamHandler() +console_handler.setFormatter(formatter) +logger.addHandler(console_handler) + +logger.setLevel(logging.DEBUG) +logger.propagate = False + + +# Function to handle warnings +def warning_handler(message, category, filename, lineno, file=None, line=None): + logger.warning(f'{filename}:{lineno}: {category.__name__}: {message}') + +# Set the warnings' showwarning function to the handler +warnings.showwarning = warning_handler \ No newline at end of file diff --git a/Tests/test_geo.py b/Tests/test_geo.py index 3c7e25a9..899f2864 100644 --- a/Tests/test_geo.py +++ b/Tests/test_geo.py @@ -376,3 +376,57 @@ def test_get_node_neighbours_per_domain(self): # Check if node neighbours are the same assert_array1D(node_neighbours_test, np.concatenate(node_neighbours_expected)) + def test_geometry_is_correct(self): + """ + Test the function geometry_is_correct + :return: + """ + # Load data + 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') + + # Check if geometry is correct + self.assertTrue(vModel_test.geo.geometry_is_correct()) + + # Load data + vModel_test = load_data('geometry_correct_3.pkl') + + # Check if geometry is correct + self.assertTrue(vModel_test.geo.geometry_is_correct()) + + def test_geometry_is_incorrect(self): + """ + Test the function geometry_is_correct + :return: + """ + # Load data + 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') + + # 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') + + # 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') + + # Check if geometry is correct + self.assertFalse(vModel_test.geo.geometry_is_correct()) + + + diff --git a/Tests/test_kg.py b/Tests/test_kg.py index 276707e2..a6eb3a2b 100644 --- a/Tests/test_kg.py +++ b/Tests/test_kg.py @@ -9,7 +9,7 @@ from pyVertexModel.Kg.kgTriEnergyBarrier import KgTriEnergyBarrier from pyVertexModel.Kg.kgViscosity import KgViscosity from pyVertexModel.Kg.kgVolume import KgVolume -from pyVertexModel.algorithm.newtonRaphson import KgGlobal +from pyVertexModel.algorithm.integrators import KgGlobal from pyVertexModel.geometry.geo import Geo @@ -18,7 +18,7 @@ def test_kg_global_filename(filename): geo_n_test = Geo(mat_info['Geo_n']) geo_0_test = Geo(mat_info['Geo_0']) # Compute the global K, and g - g, K, E, _ = KgGlobal(geo_0_test, geo_n_test, geo_test, set_test) + g, K, E, _ = KgGlobal(geo_test, set_test, geo_n_test) # Get the expected results from the mat file g_expected = mat_info['g'][:, 0] k_expected = mat_info['K'] @@ -140,7 +140,7 @@ def test_KgGlobal(self): geo_n_test = Geo(mat_info['Geo_n']) geo_0_test = Geo(mat_info['Geo_0']) - g, K, E, _ = KgGlobal(geo_0_test, geo_n_test, geo_test, set_test) + g, K, E, _ = KgGlobal(geo_test, set_test, geo_n_test) g_expected = (mat_expected['gs_full'] + mat_expected['gv_full'] + mat_expected['gf_full'] + mat_expected['gBA_full'] + mat_expected['gBAR_full'] + mat_expected['gC_full'] + diff --git a/Tests/test_newtonRaphson.py b/Tests/test_newtonRaphson.py index c67c8322..9c66fbd9 100644 --- a/Tests/test_newtonRaphson.py +++ b/Tests/test_newtonRaphson.py @@ -2,7 +2,7 @@ from Tests.test_geo import check_if_cells_are_the_same from Tests.tests import load_data, Tests, assert_array1D, assert_matrix -from pyVertexModel.algorithm.newtonRaphson import line_search, newton_raphson, newton_raphson_iteration, ml_divide, \ +from pyVertexModel.algorithm.integrators import line_search, newton_raphson, newton_raphson_iteration, ml_divide, \ solve_remodeling_step from pyVertexModel.geometry.degreesOfFreedom import DegreesOfFreedom from pyVertexModel.geometry.geo import Geo @@ -106,7 +106,7 @@ def test_newton_raphson_wingdisc(self): t_test = mat_info['t'][0][0] Set.iter = 1000000 - geo_test, g_test, k_test, energy_test, set_test, gr_test, dyr_test, dy_test = ( + geo_test, g_test, k_test, energy_test, set_test, gr_test, dyr_test, dy_test, _ = ( newton_raphson(geo_0_test, geo_n_test, geo_test, @@ -143,7 +143,7 @@ def test_newton_raphson(self): t_test = mat_info['t'][0][0] Set.iter = 1000000 - geo_test, g_test, k_test, energy_test, set_test, gr_test, dyr_test, dy_test = ( + geo_test, g_test, k_test, energy_test, set_test, gr_test, dyr_test, dy_test, _ = ( newton_raphson(geo_0_test, geo_n_test, geo_test, @@ -177,7 +177,7 @@ def test_line_search(self): g_test = mat_info['g'][:, 0] dy_test = mat_info['dy'][:, 0] set_test = Set(mat_info['Set']) - alpha = line_search(geo_0_test, geo_n_test, geo_test, dofs_test.Free, set_test, g_test, dy_test) + alpha = line_search(geo_test, dofs_test.Free, set_test, g_test, dy_test, geo_n_test) np.testing.assert_almost_equal(alpha, 1) @@ -202,7 +202,7 @@ def test_line_search_cyst(self): g_test = mat_info['g'][:, 0] dy_test = mat_info['dy'] set_test = Set(mat_info['Set']) - alpha = line_search(geo_0_test, geo_n_test, geo_test, dofs_test.Free, set_test, g_test, dy_test) + alpha = line_search(geo_test, dofs_test.Free, set_test, g_test, dy_test, geo_n_test) np.testing.assert_almost_equal(alpha, mat_info['alpha'][0][0]) diff --git a/Tests/test_vertexModel.py b/Tests/test_vertexModel.py index 5f30fdcf..d72d438c 100644 --- a/Tests/test_vertexModel.py +++ b/Tests/test_vertexModel.py @@ -7,8 +7,8 @@ from Tests import TEST_DIRECTORY from Tests.test_geo import check_if_cells_are_the_same from Tests.tests import Tests, assert_matrix, load_data, assert_array1D -from pyVertexModel.algorithm import newtonRaphson -from pyVertexModel.algorithm.newtonRaphson import newton_raphson +from pyVertexModel.algorithm import integrators +from pyVertexModel.algorithm.integrators import newton_raphson from pyVertexModel.algorithm.vertexModel import create_tetrahedra from pyVertexModel.algorithm.vertexModelBubbles import build_topo, SeedWithBoundingBox, generate_first_ghost_nodes, \ delaunay_compute_entities, VertexModelBubbles @@ -239,21 +239,6 @@ def test_build_triplets_of_neighs(self): # Check if triplets of neighbours are correct assert_matrix(triplets_of_neighs_test, mat_info['neighboursVertices']) - def test_calculate_neighbours(self): - """ - Test the calculate_neighbours function. - :return: - """ - # Load data - _, _, mat_info = load_data('calculate_neighbours_wingdisc.mat') - - neighbours_test = calculate_neighbours(mat_info['labelledImg'], 2) - - neighbours_expected = [np.concatenate(neighbours[0]) for neighbours in mat_info['imgNeighbours']] - - # Check if the cells are initialized correctly - np.testing.assert_equal(neighbours_test[1:], neighbours_expected) - def test_obtain_initial_x_and_tetrahedra(self): """ Test the obtain_initial_x_and_tetrahedra function. @@ -266,11 +251,11 @@ def test_obtain_initial_x_and_tetrahedra(self): vModel_test = VertexModelVoronoiFromTimeImage(set_test) file_name = 'LblImg_imageSequence.mat' - test_dir = 'Tests/data/%s' % file_name + 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('data/%s' % file_name) + Twg_test, X_test = vModel_test.obtain_initial_x_and_tetrahedra('Tests_data/%s' % file_name) # Check if the test and expected are the same assert_matrix(Twg_test, mat_info_expected['Twg'] - 1) @@ -326,83 +311,6 @@ def test_add_tetrahedra_intercalations(self): # Check if the test and expected are the same assert_matrix(Twg, mat_info_expected['Twg']) - def test_build_2d_voronoi_from_image(self): - """ - Test the build_2d_voronoi_from_image function. - :return: - """ - # Load data - _, _, mat_info = load_data('build_2d_voronoi_from_image_wingdisc.mat') - labelled_img = mat_info['labelledImg'] - watershed_img = mat_info['watershedImg'] - main_cells = mat_info['mainCells'][0] - - # Test if initialize geometry function does not change anything - (triangles_connectivity, neighbours_network, cell_edges, vertices_location, border_cells, - border_of_border_cells_and_main_cells) = build_2d_voronoi_from_image(labelled_img, watershed_img, main_cells) - - # Load expected - _, _, mat_info_expected = load_data('build_2d_voronoi_from_image_wingdisc_expected.mat') - - # Assert - np.testing.assert_equal(triangles_connectivity, mat_info_expected['trianglesConnectivity']) - np.testing.assert_equal(neighbours_network, mat_info_expected['neighboursNetwork']) - np.testing.assert_equal([cell_edge+1 for cell_edge in cell_edges if cell_edge is not None], [cell_edge[0] for cell_edge in mat_info_expected['cellEdges']]) - np.testing.assert_equal(border_cells, np.concatenate(mat_info_expected['borderCells'])) - np.testing.assert_equal(border_of_border_cells_and_main_cells, mat_info_expected['borderOfborderCellsAndMainCells'][0]) - - def test_populate_vertices_info(self): - """ - Test the populate_vertices_info function. - :return: - """ - # Load data - _, _, mat_info = load_data('populate_vertices_info_wingdisc.mat') - - # Load data - border_cells_and_main_cells = [border_cell[0] for border_cell in mat_info['borderCellsAndMainCells']] - labelled_img = mat_info['labelledImg'] - img_neighbours_all = [np.concatenate(neighbours[0]) for neighbours in mat_info['imgNeighbours']] - main_cells = mat_info['mainCells'][0] - ratio = 2 - - img_neighbours_all.insert(0, None) - - vertices_info_test = populate_vertices_info(border_cells_and_main_cells, img_neighbours_all, labelled_img, - main_cells, ratio) - - vertices_info_expected_per_cell = [np.concatenate(vertices[0]) for vertices in mat_info['verticesInfo']['PerCell'][0][0] if len(vertices[0][0]) > 0] - vertices_info_expected_edges = [np.concatenate(vertices[0]) for vertices in mat_info['verticesInfo']['edges'][0][0] if len(vertices[0][0]) > 0] - - # Assert - np.testing.assert_equal([vertices + 1 for vertices in vertices_info_test['PerCell'] if vertices is not None], vertices_info_expected_per_cell) - np.testing.assert_equal([np.concatenate(edges) + 1 for edges in vertices_info_test['edges'] if edges is not None], vertices_info_expected_edges) - np.testing.assert_equal(vertices_info_test['connectedCells'], mat_info['verticesInfo']['connectedCells'][0][0]) - - def test_calculate_vertices(self): - """ - Test the calculate_vertices function. - :return: - """ - # Load data - _, _, mat_info = load_data('calculate_vertices_wingdisc.mat') - - # Load data - labelled_img = mat_info['labelledImg'] - img_neighbours_all = [np.concatenate(neighbours[0]) for neighbours in mat_info['neighbours']] - ratio = 2 - - img_neighbours_all.insert(0, None) - - # Test if initialize geometry function does not change anything - vertices_info_test = calculate_vertices(labelled_img, img_neighbours_all, ratio) - - # Load expected - _, _, mat_info_expected = load_data('calculate_vertices_wingdisc_expected.mat') - - # Assert - assert_matrix(vertices_info_test['connectedCells'], mat_info_expected['verticesInfo']['connectedCells'][0][0]) - def test_get_four_fold_vertices(self): """ Test the get_four_fold_vertices function. @@ -453,11 +361,11 @@ def test_process_image(self): """ # Process image file_name = 'LblImg_imageSequence.mat' - test_dir = 'Tests/data/%s' % file_name + test_dir = 'Tests/Tests_data/%s' % file_name if exists(test_dir): _, imgStackLabelled_test = process_image(test_dir) else: - _, imgStackLabelled_test = process_image('data/%s' % file_name) + _, imgStackLabelled_test = process_image('Tests_data/%s' % file_name) # Load expected _, _, mat_info_expected = load_data('process_image_wingdisc_expected.mat') @@ -476,11 +384,11 @@ 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/data/%s' % file_name + 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 = 'data/%s' % file_name + vModel_test.set.initial_filename_state = 'Tests_data/%s' % file_name vModel_test.initialize() @@ -489,8 +397,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_0, vModel_test.geo, vModel_test.geo, - vModel_test.set) + g_test, K_test, energies_test, _ = newtonRaphson.KgGlobal(vModel_test.geo, vModel_test.set, vModel_test.geo) # Check if energies are the same assert_array1D(g_test, mat_info['g']) @@ -577,3 +484,76 @@ def test_initialize_with_numpy_array(self): assert vModel_test.geo.nCells > 0, "Should have cells" + def test_weird_bug_should_not_happen(self): + """ + Test for a weird bug that should not happen. + Uses FIRE algorithm for stable, fast convergence. + :return: + """ + # Load data + 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 + + # Update tolerance + vModel_test.set.dt_tolerance = 0.25 + + # Use FIRE algorithm for better stability and convergence + vModel_test.set.integrator = 'fire' + + # Initialize FIRE parameters tuned for edge case geometries + # More lenient parameters to handle near-zero power cases + vModel_test.set.fire_dt_max = 10 * vModel_test.set.dt0 + vModel_test.set.fire_dt_min = 0.1 * vModel_test.set.dt0 # Higher minimum (was 0.02) + vModel_test.set.fire_N_min = 3 # Accelerate sooner (was 5) + vModel_test.set.fire_f_inc = 1.2 # Faster acceleration (was 1.1) + vModel_test.set.fire_f_dec = 0.7 # Less aggressive decrease (was 0.5) + vModel_test.set.fire_alpha_start = 0.15 # More damping initially (was 0.1) + vModel_test.set.fire_f_alpha = 0.98 # Slower damping reduction (was 0.99) + + # Run the model + vModel_test.iterate_over_time() + + # Check if it did not converge + self.assertFalse(vModel_test.didNotConverge) + + def test_vertices_shouldnt_be_going_wild(self): + """ + Test for another weird bug that should not happen. + :return: + """ + # Load data + vModel_test = load_data('vertices_going_wild_2.pkl') + + # Run for 10 iterations. dt should not decrease to 1e-1 + vModel_test.set.tend = vModel_test.t + 20 * vModel_test.set.dt0 + + # Update tolerance + vModel_test.set.dt_tolerance = 0.25 + + # Run the model + vModel_test.iterate_over_time() + + # Check if it did not converge + self.assertFalse(vModel_test.didNotConverge) + + def test_vertices_shouldnt_be_going_wild_3(self): + """ + Test for another weird bug that should not happen. + :return: + """ + # Load data + vModel_test = load_data('vertices_going_wild_3.pkl') + + # Run for 10 iterations. dt should not decrease to 1e-1 + vModel_test.set.tend = vModel_test.t + 20 * vModel_test.set.dt0 + + # Update tolerance + vModel_test.set.dt_tolerance = 0.25 + + # Run the model + vModel_test.iterate_over_time() + + # Check if it did not converge + self.assertFalse(vModel_test.didNotConverge) \ No newline at end of file diff --git a/Tests/tests.py b/Tests/tests.py index eef2d05d..f49b9107 100644 --- a/Tests/tests.py +++ b/Tests/tests.py @@ -1,36 +1,57 @@ import unittest -import scipy.io -import numpy as np from os.path import exists, abspath +import numpy as np +import scipy.io + +from pyVertexModel.algorithm.vertexModelVoronoiFromTimeImage import VertexModelVoronoiFromTimeImage from pyVertexModel.geometry.geo import Geo from pyVertexModel.parameters.set import Set -from pyVertexModel.Kg import kg_functions +from pyVertexModel.util.utils import load_state + def load_data(file_name, return_geo=True): - test_dir = abspath('Tests/data/%s' % file_name) - if exists(test_dir): - mat_info = scipy.io.loadmat(test_dir) - else: - mat_info = scipy.io.loadmat('data/%s' % file_name) + test_dir = abspath('Tests/Tests_data/%s' % 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) + + if return_geo: + if 'Geo' in mat_info.keys(): + geo_test = Geo(mat_info['Geo']) + else: + geo_test = None - if return_geo: - if 'Geo' in mat_info.keys(): - geo_test = Geo(mat_info['Geo']) + if 'Set' in mat_info.keys(): + set_test = Set(mat_info['Set']) + if set_test.OutputFolder.__eq__(b'') or set_test.OutputFolder is None: + set_test.OutputFolder = '../Result/Test' + else: + set_test = None else: geo_test = None + set_test = None - if 'Set' in mat_info.keys(): - set_test = Set(mat_info['Set']) - if set_test.OutputFolder.__eq__(b'') or set_test.OutputFolder is None: - set_test.OutputFolder = '../Result/Test' + 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: - set_test = None + load_state(v_model, 'Tests_data/%s' % file_name) + + # 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 else: - geo_test = None - set_test = None + raise FileNotFoundError('File %s not found' % file_name) + - return geo_test, set_test, mat_info def assert_matrix(k_expected, k): diff --git a/Tests_data/Geo_3x3_dofs_expected.mat b/Tests_data/Geo_3x3_dofs_expected.mat new file mode 100644 index 00000000..693213c5 Binary files /dev/null and b/Tests_data/Geo_3x3_dofs_expected.mat differ diff --git a/Tests_data/Geo_var_3x3_stretch.mat b/Tests_data/Geo_var_3x3_stretch.mat new file mode 100644 index 00000000..425a7919 Binary files /dev/null and b/Tests_data/Geo_var_3x3_stretch.mat differ diff --git a/Tests_data/Geo_var_3x3_stretch_Iter6_Kgs_expectedResults.mat b/Tests_data/Geo_var_3x3_stretch_Iter6_Kgs_expectedResults.mat new file mode 100644 index 00000000..1806da60 Binary files /dev/null and b/Tests_data/Geo_var_3x3_stretch_Iter6_Kgs_expectedResults.mat differ diff --git a/Tests_data/Geo_var_3x3_stretch_Iter6_expectedResults.mat b/Tests_data/Geo_var_3x3_stretch_Iter6_expectedResults.mat new file mode 100644 index 00000000..9bfb41bd Binary files /dev/null and b/Tests_data/Geo_var_3x3_stretch_Iter6_expectedResults.mat differ diff --git a/Tests_data/Geo_var_3x3_stretch_expectedResults.mat b/Tests_data/Geo_var_3x3_stretch_expectedResults.mat new file mode 100644 index 00000000..f978aa2f Binary files /dev/null and b/Tests_data/Geo_var_3x3_stretch_expectedResults.mat differ diff --git a/Tests_data/Geo_var_cyst.mat b/Tests_data/Geo_var_cyst.mat new file mode 100644 index 00000000..cb003a54 Binary files /dev/null and b/Tests_data/Geo_var_cyst.mat differ diff --git a/Tests_data/Geo_var_cyst_expectedResults_cell1.mat b/Tests_data/Geo_var_cyst_expectedResults_cell1.mat new file mode 100644 index 00000000..2d3d7a7a Binary files /dev/null and b/Tests_data/Geo_var_cyst_expectedResults_cell1.mat differ diff --git a/Tests_data/LblImg_imageSequence.mat b/Tests_data/LblImg_imageSequence.mat new file mode 100644 index 00000000..d266e0b2 Binary files /dev/null and b/Tests_data/LblImg_imageSequence.mat differ diff --git a/Tests_data/Newton_Raphson_3x3_stretch.mat b/Tests_data/Newton_Raphson_3x3_stretch.mat new file mode 100644 index 00000000..9e43d968 Binary files /dev/null and b/Tests_data/Newton_Raphson_3x3_stretch.mat differ diff --git a/Tests_data/Newton_Raphson_3x3_stretch_expected.mat b/Tests_data/Newton_Raphson_3x3_stretch_expected.mat new file mode 100644 index 00000000..2a5702e1 Binary files /dev/null and b/Tests_data/Newton_Raphson_3x3_stretch_expected.mat differ diff --git a/Tests_data/Newton_Raphson_Iteration_3x3_stretch_expected.mat b/Tests_data/Newton_Raphson_Iteration_3x3_stretch_expected.mat new file mode 100644 index 00000000..65f1048a Binary files /dev/null and b/Tests_data/Newton_Raphson_Iteration_3x3_stretch_expected.mat differ diff --git a/Tests_data/Newton_Raphson_Iteration_wingdisc.mat b/Tests_data/Newton_Raphson_Iteration_wingdisc.mat new file mode 100644 index 00000000..e8e10065 Binary files /dev/null and b/Tests_data/Newton_Raphson_Iteration_wingdisc.mat differ diff --git a/Tests_data/Newton_Raphson_Iteration_wingdisc_expected.mat b/Tests_data/Newton_Raphson_Iteration_wingdisc_expected.mat new file mode 100644 index 00000000..fe3efc4d Binary files /dev/null and b/Tests_data/Newton_Raphson_Iteration_wingdisc_expected.mat differ diff --git a/Tests_data/Newton_Raphson_ml_divide_3x3_stretch_expected.mat b/Tests_data/Newton_Raphson_ml_divide_3x3_stretch_expected.mat new file mode 100644 index 00000000..31ebcccf Binary files /dev/null and b/Tests_data/Newton_Raphson_ml_divide_3x3_stretch_expected.mat differ diff --git a/Tests_data/Newton_Raphson_wingdisc_expected.mat b/Tests_data/Newton_Raphson_wingdisc_expected.mat new file mode 100644 index 00000000..843698b1 Binary files /dev/null and b/Tests_data/Newton_Raphson_wingdisc_expected.mat differ diff --git a/Tests_data/add_and_rebuild_cells_wingdisc.mat b/Tests_data/add_and_rebuild_cells_wingdisc.mat new file mode 100644 index 00000000..90769170 Binary files /dev/null and b/Tests_data/add_and_rebuild_cells_wingdisc.mat differ diff --git a/Tests_data/add_and_rebuild_cells_wingdisc_expected.mat b/Tests_data/add_and_rebuild_cells_wingdisc_expected.mat new file mode 100644 index 00000000..aa373238 Binary files /dev/null and b/Tests_data/add_and_rebuild_cells_wingdisc_expected.mat differ diff --git a/Tests_data/add_tetrahedra_from_intercalations_wingdisc.mat b/Tests_data/add_tetrahedra_from_intercalations_wingdisc.mat new file mode 100644 index 00000000..50888996 Binary files /dev/null and b/Tests_data/add_tetrahedra_from_intercalations_wingdisc.mat differ diff --git a/Tests_data/add_tetrahedra_from_intercalations_wingdisc_expected.mat b/Tests_data/add_tetrahedra_from_intercalations_wingdisc_expected.mat new file mode 100644 index 00000000..3130f49d Binary files /dev/null and b/Tests_data/add_tetrahedra_from_intercalations_wingdisc_expected.mat differ diff --git a/Tests_data/add_tetrahedra_wingdisc.mat b/Tests_data/add_tetrahedra_wingdisc.mat new file mode 100644 index 00000000..1d535a1a Binary files /dev/null and b/Tests_data/add_tetrahedra_wingdisc.mat differ diff --git a/Tests_data/add_tetrahedra_wingdisc_expected.mat b/Tests_data/add_tetrahedra_wingdisc_expected.mat new file mode 100644 index 00000000..ef3bfae9 Binary files /dev/null and b/Tests_data/add_tetrahedra_wingdisc_expected.mat differ diff --git a/Tests_data/build_2d_voronoi_from_image_wingdisc.mat b/Tests_data/build_2d_voronoi_from_image_wingdisc.mat new file mode 100644 index 00000000..d13b3114 Binary files /dev/null and b/Tests_data/build_2d_voronoi_from_image_wingdisc.mat differ diff --git a/Tests_data/build_2d_voronoi_from_image_wingdisc_expected.mat b/Tests_data/build_2d_voronoi_from_image_wingdisc_expected.mat new file mode 100644 index 00000000..9ae39346 Binary files /dev/null and b/Tests_data/build_2d_voronoi_from_image_wingdisc_expected.mat differ diff --git a/Tests_data/build_cells_cyst.mat b/Tests_data/build_cells_cyst.mat new file mode 100644 index 00000000..a586c826 Binary files /dev/null and b/Tests_data/build_cells_cyst.mat differ diff --git a/Tests_data/build_cells_cyst_expected.mat b/Tests_data/build_cells_cyst_expected.mat new file mode 100644 index 00000000..537c3927 Binary files /dev/null and b/Tests_data/build_cells_cyst_expected.mat differ diff --git a/Tests_data/build_face_cyst.mat b/Tests_data/build_face_cyst.mat new file mode 100644 index 00000000..cd2d4592 Binary files /dev/null and b/Tests_data/build_face_cyst.mat differ diff --git a/Tests_data/build_face_cyst_expected.mat b/Tests_data/build_face_cyst_expected.mat new file mode 100644 index 00000000..1fe97027 Binary files /dev/null and b/Tests_data/build_face_cyst_expected.mat differ diff --git a/Tests_data/build_globald_ids_wingdisc_expected.mat b/Tests_data/build_globald_ids_wingdisc_expected.mat new file mode 100644 index 00000000..563eea82 Binary files /dev/null and b/Tests_data/build_globald_ids_wingdisc_expected.mat differ diff --git a/Tests_data/build_topo_expected.mat b/Tests_data/build_topo_expected.mat new file mode 100644 index 00000000..8d2d9d38 Binary files /dev/null and b/Tests_data/build_topo_expected.mat differ diff --git a/Tests_data/build_triplets_wingdisc.mat b/Tests_data/build_triplets_wingdisc.mat new file mode 100644 index 00000000..3d57930a Binary files /dev/null and b/Tests_data/build_triplets_wingdisc.mat differ diff --git a/Tests_data/build_x_from_y_wingdisc.mat b/Tests_data/build_x_from_y_wingdisc.mat new file mode 100644 index 00000000..3cb18c36 Binary files /dev/null and b/Tests_data/build_x_from_y_wingdisc.mat differ diff --git a/Tests_data/build_x_from_y_wingdisc_expected.mat b/Tests_data/build_x_from_y_wingdisc_expected.mat new file mode 100644 index 00000000..8da1083b Binary files /dev/null and b/Tests_data/build_x_from_y_wingdisc_expected.mat differ diff --git a/Tests_data/calculate_neighbours_wingdisc.mat b/Tests_data/calculate_neighbours_wingdisc.mat new file mode 100644 index 00000000..3118094c Binary files /dev/null and b/Tests_data/calculate_neighbours_wingdisc.mat differ diff --git a/Tests_data/calculate_vertices_wingdisc.mat b/Tests_data/calculate_vertices_wingdisc.mat new file mode 100644 index 00000000..b3816367 Binary files /dev/null and b/Tests_data/calculate_vertices_wingdisc.mat differ diff --git a/Tests_data/calculate_vertices_wingdisc_expected.mat b/Tests_data/calculate_vertices_wingdisc_expected.mat new file mode 100644 index 00000000..89fec459 Binary files /dev/null and b/Tests_data/calculate_vertices_wingdisc_expected.mat differ diff --git a/Tests_data/check_ys_and_faces_have_not_changed_wingdisc.mat b/Tests_data/check_ys_and_faces_have_not_changed_wingdisc.mat new file mode 100644 index 00000000..a183e1c2 Binary files /dev/null and b/Tests_data/check_ys_and_faces_have_not_changed_wingdisc.mat differ diff --git a/Tests_data/create_tetrahedra_wingdisc.mat b/Tests_data/create_tetrahedra_wingdisc.mat new file mode 100644 index 00000000..e008e42d Binary files /dev/null and b/Tests_data/create_tetrahedra_wingdisc.mat differ diff --git a/Tests_data/delaunay_compute_entities_expected.mat b/Tests_data/delaunay_compute_entities_expected.mat new file mode 100644 index 00000000..adc16d1c Binary files /dev/null and b/Tests_data/delaunay_compute_entities_expected.mat differ diff --git a/Tests_data/delaunay_compute_entities_input.mat b/Tests_data/delaunay_compute_entities_input.mat new file mode 100644 index 00000000..8e89eaa0 Binary files /dev/null and b/Tests_data/delaunay_compute_entities_input.mat differ diff --git a/Tests_data/delaunay_input_expected.mat b/Tests_data/delaunay_input_expected.mat new file mode 100644 index 00000000..13b9e4ca Binary files /dev/null and b/Tests_data/delaunay_input_expected.mat differ diff --git a/Tests_data/delaunay_output_cyst.mat b/Tests_data/delaunay_output_cyst.mat new file mode 100644 index 00000000..821788d5 Binary files /dev/null and b/Tests_data/delaunay_output_cyst.mat differ diff --git a/Tests_data/flipNM_cyst.mat b/Tests_data/flipNM_cyst.mat new file mode 100644 index 00000000..157ff9ff Binary files /dev/null and b/Tests_data/flipNM_cyst.mat differ diff --git a/Tests_data/flipNM_recursive_cyst.mat b/Tests_data/flipNM_recursive_cyst.mat new file mode 100644 index 00000000..a913e187 Binary files /dev/null and b/Tests_data/flipNM_recursive_cyst.mat differ diff --git a/Tests_data/flipNM_recursive_cyst_expected.mat b/Tests_data/flipNM_recursive_cyst_expected.mat new file mode 100644 index 00000000..1b32d31d Binary files /dev/null and b/Tests_data/flipNM_recursive_cyst_expected.mat differ diff --git a/Tests_data/generate_first_ghost_nodes_expected.mat b/Tests_data/generate_first_ghost_nodes_expected.mat new file mode 100644 index 00000000..24439994 Binary files /dev/null and b/Tests_data/generate_first_ghost_nodes_expected.mat differ diff --git a/Tests_data/geo_cyst_expected_extrapolatedYs.mat b/Tests_data/geo_cyst_expected_extrapolatedYs.mat new file mode 100644 index 00000000..c161c19a Binary files /dev/null and b/Tests_data/geo_cyst_expected_extrapolatedYs.mat differ diff --git a/Tests_data/geometry_correct_1.pkl b/Tests_data/geometry_correct_1.pkl new file mode 100644 index 00000000..74b8744e Binary files /dev/null and b/Tests_data/geometry_correct_1.pkl differ diff --git a/Tests_data/geometry_correct_2.pkl b/Tests_data/geometry_correct_2.pkl new file mode 100644 index 00000000..de99fabc Binary files /dev/null and b/Tests_data/geometry_correct_2.pkl differ diff --git a/Tests_data/geometry_correct_3.pkl b/Tests_data/geometry_correct_3.pkl new file mode 100644 index 00000000..4e5c616c Binary files /dev/null and b/Tests_data/geometry_correct_3.pkl differ diff --git a/Tests_data/get_dofs_remodel.mat b/Tests_data/get_dofs_remodel.mat new file mode 100644 index 00000000..e4ad92af Binary files /dev/null and b/Tests_data/get_dofs_remodel.mat differ diff --git a/Tests_data/get_four_fold_vertices_wingdisc.mat b/Tests_data/get_four_fold_vertices_wingdisc.mat new file mode 100644 index 00000000..9e425d53 Binary files /dev/null and b/Tests_data/get_four_fold_vertices_wingdisc.mat differ diff --git a/Tests_data/get_node_neighbours_per_domain_wingdisc.mat b/Tests_data/get_node_neighbours_per_domain_wingdisc.mat new file mode 100644 index 00000000..2c56fb60 Binary files /dev/null and b/Tests_data/get_node_neighbours_per_domain_wingdisc.mat differ diff --git a/Tests_data/initialize_cells_cyst_expected.mat b/Tests_data/initialize_cells_cyst_expected.mat new file mode 100644 index 00000000..0c8bbd32 Binary files /dev/null and b/Tests_data/initialize_cells_cyst_expected.mat differ diff --git a/Tests_data/initialize_voronoi_wingdisc.mat b/Tests_data/initialize_voronoi_wingdisc.mat new file mode 100644 index 00000000..83eee65d Binary files /dev/null and b/Tests_data/initialize_voronoi_wingdisc.mat differ diff --git a/Tests_data/kK_test.mat b/Tests_data/kK_test.mat new file mode 100644 index 00000000..34015607 Binary files /dev/null and b/Tests_data/kK_test.mat differ diff --git a/Tests_data/line_search_cyst.mat b/Tests_data/line_search_cyst.mat new file mode 100644 index 00000000..febeea9c Binary files /dev/null and b/Tests_data/line_search_cyst.mat differ diff --git a/Tests_data/obtain_x_and_twg_wingdisc.mat b/Tests_data/obtain_x_and_twg_wingdisc.mat new file mode 100644 index 00000000..ffc49517 Binary files /dev/null and b/Tests_data/obtain_x_and_twg_wingdisc.mat differ diff --git a/Tests_data/populate_vertices_info_wingdisc.mat b/Tests_data/populate_vertices_info_wingdisc.mat new file mode 100644 index 00000000..48d032fc Binary files /dev/null and b/Tests_data/populate_vertices_info_wingdisc.mat differ diff --git a/Tests_data/post_flip_wingdisc.mat b/Tests_data/post_flip_wingdisc.mat new file mode 100644 index 00000000..b038804f Binary files /dev/null and b/Tests_data/post_flip_wingdisc.mat differ diff --git a/Tests_data/post_flip_wingdisc_2.mat b/Tests_data/post_flip_wingdisc_2.mat new file mode 100644 index 00000000..078b96b4 Binary files /dev/null and b/Tests_data/post_flip_wingdisc_2.mat differ diff --git a/Tests_data/post_flip_wingdisc_2_expected.mat b/Tests_data/post_flip_wingdisc_2_expected.mat new file mode 100644 index 00000000..30cd404f Binary files /dev/null and b/Tests_data/post_flip_wingdisc_2_expected.mat differ diff --git a/Tests_data/post_flip_wingdisc_expected.mat b/Tests_data/post_flip_wingdisc_expected.mat new file mode 100644 index 00000000..18221b41 Binary files /dev/null and b/Tests_data/post_flip_wingdisc_expected.mat differ diff --git a/Tests_data/process_image_wingdisc_expected.mat b/Tests_data/process_image_wingdisc_expected.mat new file mode 100644 index 00000000..f919853d Binary files /dev/null and b/Tests_data/process_image_wingdisc_expected.mat differ diff --git a/Tests_data/rebuild_stretch_3x3_expected.mat b/Tests_data/rebuild_stretch_3x3_expected.mat new file mode 100644 index 00000000..774a663e Binary files /dev/null and b/Tests_data/rebuild_stretch_3x3_expected.mat differ diff --git a/Tests_data/rebuild_wingdisc_expected.mat b/Tests_data/rebuild_wingdisc_expected.mat new file mode 100644 index 00000000..d3aabb18 Binary files /dev/null and b/Tests_data/rebuild_wingdisc_expected.mat differ diff --git a/Tests_data/remodelling_intercalating_cells_wingdisc.mat b/Tests_data/remodelling_intercalating_cells_wingdisc.mat new file mode 100644 index 00000000..62caf47d Binary files /dev/null and b/Tests_data/remodelling_intercalating_cells_wingdisc.mat differ diff --git a/Tests_data/remodelling_intercalating_cells_wingdisc_expected.mat b/Tests_data/remodelling_intercalating_cells_wingdisc_expected.mat new file mode 100644 index 00000000..871cd8e4 Binary files /dev/null and b/Tests_data/remodelling_intercalating_cells_wingdisc_expected.mat differ diff --git a/Tests_data/remove_tetrahedra_wingdisc.mat b/Tests_data/remove_tetrahedra_wingdisc.mat new file mode 100644 index 00000000..74dd62f2 Binary files /dev/null and b/Tests_data/remove_tetrahedra_wingdisc.mat differ diff --git a/Tests_data/seed_with_bbox_expected.mat b/Tests_data/seed_with_bbox_expected.mat new file mode 100644 index 00000000..197709ff Binary files /dev/null and b/Tests_data/seed_with_bbox_expected.mat differ diff --git a/Tests_data/seed_with_bbox_input.mat b/Tests_data/seed_with_bbox_input.mat new file mode 100644 index 00000000..4e7c12c9 Binary files /dev/null and b/Tests_data/seed_with_bbox_input.mat differ diff --git a/Tests_data/solve_remodelling_step_wingdisc.mat b/Tests_data/solve_remodelling_step_wingdisc.mat new file mode 100644 index 00000000..e7bc6366 Binary files /dev/null and b/Tests_data/solve_remodelling_step_wingdisc.mat differ diff --git a/Tests_data/test_gkDet.mat b/Tests_data/test_gkDet.mat new file mode 100644 index 00000000..8367fff9 Binary files /dev/null and b/Tests_data/test_gkDet.mat differ diff --git a/Tests_data/vertices_going_wild_1.pkl b/Tests_data/vertices_going_wild_1.pkl new file mode 100644 index 00000000..bdf9f86e Binary files /dev/null and b/Tests_data/vertices_going_wild_1.pkl differ diff --git a/Tests_data/vertices_going_wild_2.pkl b/Tests_data/vertices_going_wild_2.pkl new file mode 100644 index 00000000..17715ed1 Binary files /dev/null and b/Tests_data/vertices_going_wild_2.pkl differ diff --git a/Tests_data/vertices_going_wild_3.pkl b/Tests_data/vertices_going_wild_3.pkl new file mode 100644 index 00000000..fd0cf758 Binary files /dev/null and b/Tests_data/vertices_going_wild_3.pkl differ diff --git a/Tests_data/vertices_going_wild_4.pkl b/Tests_data/vertices_going_wild_4.pkl new file mode 100644 index 00000000..aa51edc0 Binary files /dev/null and b/Tests_data/vertices_going_wild_4.pkl differ diff --git a/Tests_data/voronoi_40cells.pkl b/Tests_data/voronoi_40cells.pkl new file mode 100644 index 00000000..885360af Binary files /dev/null and b/Tests_data/voronoi_40cells.pkl differ diff --git a/Tests_data/wing_disc_110.mat b/Tests_data/wing_disc_110.mat new file mode 100644 index 00000000..16fd31df Binary files /dev/null and b/Tests_data/wing_disc_110.mat differ diff --git a/Tests_data/wing_disc_150.mat b/Tests_data/wing_disc_150.mat new file mode 100644 index 00000000..a91d3036 Binary files /dev/null and b/Tests_data/wing_disc_150.mat differ diff --git a/pyproject.toml b/pyproject.toml index 101e41b9..3acf313f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,8 @@ dependencies = [ "scikit-learn", "matplotlib", "seaborn", + "opencv-python>=4.13.0.92", + "openpyxl>=3.1.5", ] [tool.setuptools_scm] @@ -135,4 +137,4 @@ legacy_tox_ini = """ [tox] envlist = py3{9,10}-{linux,macos,windows} isolated_build=true -""" \ No newline at end of file +""" diff --git a/report.md b/report.md deleted file mode 100644 index 8b137891..00000000 --- a/report.md +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/pyVertexModel/__init__.py b/src/pyVertexModel/__init__.py index b7e540b5..b8cc1b0e 100644 --- a/src/pyVertexModel/__init__.py +++ b/src/pyVertexModel/__init__.py @@ -10,8 +10,6 @@ # We want the project root directory (the parent of 'src') PROJECT_DIRECTORY = os.getenv('PROJECT_DIR', str(Path(__file__).parent.parent.parent)) -__all__ = ["PROJECT_DIRECTORY"] - # get the logger instance logger = logging.getLogger("pyVertexModel") @@ -19,9 +17,11 @@ "%(levelname)s [%(asctime)s] pyVertexModel: %(message)s", datefmt="%Y/%m/%d %I:%M:%S %p", ) -console_handler = logging.StreamHandler() -console_handler.setFormatter(formatter) -logger.addHandler(console_handler) +# 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): + console_handler = logging.StreamHandler() + console_handler.setFormatter(formatter) + logger.addHandler(console_handler) logger.setLevel(logging.DEBUG) logger.propagate = False @@ -35,4 +35,4 @@ def warning_handler(message, category, filename, lineno, file=None, line=None): # Set the warnings' showwarning function to the handler warnings.showwarning = warning_handler -__all__ = ['__version__', 'PROJECT_DIRECTORY', 'logger'] \ No newline at end of file +__all__ = ['__version__', 'PROJECT_DIRECTORY'] \ No newline at end of file diff --git a/src/pyVertexModel/algorithm/newtonRaphson.py b/src/pyVertexModel/algorithm/integrators.py similarity index 58% rename from src/pyVertexModel/algorithm/newtonRaphson.py rename to src/pyVertexModel/algorithm/integrators.py index e4802eed..f38aafdf 100644 --- a/src/pyVertexModel/algorithm/newtonRaphson.py +++ b/src/pyVertexModel/algorithm/integrators.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import logging import numpy as np @@ -62,7 +64,7 @@ def solve_remodeling_step(geo_0, geo_n, geo, dofs, c_set): for n_id, nu_factor in enumerate(np.linspace(10, 1, 3)): c_set.nu0 = original_nu0 * nu_factor c_set.nu = original_nu * nu_factor - g, k, _, _ = KgGlobal(geo_0, geo_n, geo, c_set) + g, k, _, _ = KgGlobal(geo, c_set, geo_n) dy = np.zeros(((geo.numY + geo.numF + geo.nCells) * 3, 1), dtype=np.float64) dyr = np.linalg.norm(dy[dofs.remodel, 0]) @@ -71,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.') @@ -91,7 +93,6 @@ def solve_remodeling_step(geo_0, geo_n, geo, dofs, c_set): return geo, c_set, has_converged - def newton_raphson(Geo_0, Geo_n, Geo, Dofs, Set, K, g, numStep, t, implicit_method=True): """ Newton-Raphson method @@ -166,11 +167,11 @@ def newton_raphson_iteration(Dofs, Geo, Geo_0, Geo_n, K, Set, aux_gr, dof, dy, g """ dy[dof, 0] = ml_divide(K, dof, g) - alpha = line_search(Geo_0, Geo_n, Geo, dof, Set, g, dy) + alpha = line_search(Geo, dof, Set, g, dy, Geo_n) dy_reshaped = np.reshape(dy * alpha, (Geo.numF + Geo.numY + Geo.nCells, 3)) Geo.update_vertices(dy_reshaped) Geo.update_measures() - g, K, energy_total, _ = KgGlobal(Geo_0, Geo_n, Geo, Set) + g, K, energy_total, _ = KgGlobal(Geo, Set, Geo_n) dyr = np.linalg.norm(dy[dof, 0]) gr = np.linalg.norm(g[dof]) @@ -196,7 +197,7 @@ def newton_raphson_iteration(Dofs, Geo, Geo_0, Geo_n, K, Set, aux_gr, dof, dy, g def newton_raphson_iteration_explicit(Geo, Set, dof, dy, g, selected_cells=None): """ - Explicit update method + Explicit update method with adaptive step scaling for stability :param selected_cells: :param Geo_0: :param Geo_n: @@ -208,7 +209,11 @@ def newton_raphson_iteration_explicit(Geo, Set, dof, dy, g, selected_cells=None) # Bottom nodes g_constrained = constrain_bottom_vertices_x_y(Geo) - # Update the bottom nodes with the same displacement as the corresponding real nodes + # Adaptive step size scaling to prevent gradient explosion + # ALWAYS use conservative scaling to prevent any gradient growth + gr = np.linalg.norm(g[dof]) + + # Update the bottom nodes with scaled displacement dy[dof, 0] = -Set.dt / Set.nu * g[dof] dy[g_constrained, 0] = -Set.dt / Set.nu_bottom * g[g_constrained] @@ -229,6 +234,8 @@ def newton_raphson_iteration_explicit(Geo, Set, dof, dy, g, selected_cells=None) return Geo, dy, dyr + + def map_vertices_periodic_boundaries(Geo, dy): """ Update the border ghost nodes with the same displacement as the corresponding real nodes @@ -341,11 +348,10 @@ def ml_divide(K, dof, g): return -np.linalg.solve(K[np.ix_(dof, dof)], g[dof]) -def line_search(Geo_0, Geo_n, geo, dof, Set, gc, dy): +def line_search(geo, dof, Set, gc, dy, geo_n): """ Line search to find the best alpha to minimize the energy - :param Geo_0: Initial geometry. - :param Geo_n: Geometry at the previous time step + :param geo_n: Geometry at the previous time step :param geo: Current geometry :param Dofs: Dofs object :param Set: Set object @@ -361,7 +367,10 @@ def line_search(Geo_0, Geo_n, geo, dof, Set, gc, dy): Geo_copy.update_vertices(dy_reshaped) Geo_copy.update_measures() - g, _ = gGlobal(Geo_0, Geo_n, Geo_copy, Set) + if Set.implicit_method: + g, _ = g_global(Geo_copy, Set, geo_n) + else: + g, _ = g_global(Geo_copy, Set, None, False) gr0 = np.linalg.norm(gc[dof]) gr = np.linalg.norm(g[dof]) @@ -386,10 +395,9 @@ def line_search(Geo_0, Geo_n, geo, dof, Set, gc, dy): return alpha -def KgGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True): +def KgGlobal(Geo, Set, Geo_n=None, implicit_method=True): """ Compute the global Jacobian matrix and the global gradient - :param Geo_0: :param Geo_n: :param Geo: :param Set: @@ -413,20 +421,6 @@ def KgGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True): energy_total += kg_Vol.energy energies["Volume"] = kg_Vol.energy - # # TODO: Plane Elasticity - # if Set.InPlaneElasticity: - # gt, Kt, EBulk = KgBulk(geo_0, Geo, Set) - # K += Kt - # g += gt - # E += EBulk - # Energies["Bulk"] = EBulk - - # Bending Energy - # TODO - - # Propulsion Forces - # TODO - # Contractility if Set.Contractility: kg_lt = KgContractility(Geo) @@ -471,7 +465,7 @@ def KgGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True): energy_total += kg_TriAR.energy energies["TriEnergyBarrierAR"] = kg_TriAR.energy - if implicit_method is True: + if implicit_method: # Viscous Energy kg_Viscosity = KgViscosity(Geo) kg_Viscosity.compute_work(Geo, Set, Geo_n) @@ -483,95 +477,368 @@ def KgGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True): return g, K, energy_total, energies -def gGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True, num_step=None): +def g_global(geo, c_set, geo_n=None, implicit_method=True, num_step=None): """ Compute the global gradient - :param Geo_0: - :param Geo_n: - :param Geo: - :param Set: + :param num_step: + :param geo_n: + :param geo: + :param c_set: :param implicit_method: :return: """ - g = np.zeros((Geo.numY + Geo.numF + Geo.nCells) * 3, dtype=np.float64) + g = np.zeros((geo.numY + geo.numF + geo.nCells) * 3, dtype=np.float64) energies = {} # Surface Energy - if Set.lambdaS1 > 0: - kg_SA = KgSurfaceCellBasedAdhesion(Geo) - kg_SA.compute_work(Geo, Set, None, False) + if c_set.lambdaS1 > 0: + kg_SA = KgSurfaceCellBasedAdhesion(geo) + kg_SA.compute_work(geo, c_set, None, False) g += kg_SA.g - Geo.create_vtk_cell(Set, num_step, 'Arrows_surface', -kg_SA.g) + geo.create_vtk_cell(c_set, num_step, 'Arrows_surface', -kg_SA.g) energies["Surface"] = kg_SA.energy # Volume Energy - if Set.lambdaV > 0: - kg_Vol = KgVolume(Geo) - kg_Vol.compute_work(Geo, Set, None, False) + if c_set.lambdaV > 0: + kg_Vol = KgVolume(geo) + kg_Vol.compute_work(geo, c_set, None, False) g += kg_Vol.g[:] - Geo.create_vtk_cell(Set, num_step, 'Arrows_volume', -kg_Vol.g) + geo.create_vtk_cell(c_set, num_step, 'Arrows_volume', -kg_Vol.g) energies["Volume"] = kg_Vol.energy if implicit_method: # Viscous Energy - kg_Viscosity = KgViscosity(Geo) - kg_Viscosity.compute_work(Geo, Set, Geo_n, False) + kg_Viscosity = KgViscosity(geo) + kg_Viscosity.compute_work(geo, c_set, geo_n, False) g += kg_Viscosity.g - Geo.create_vtk_cell(Set, num_step, 'Arrows_viscosity', -kg_Viscosity.g) + geo.create_vtk_cell(c_set, num_step, 'Arrows_viscosity', -kg_Viscosity.g) energies["Viscosity"] = kg_Viscosity.energy - # # TODO: Plane Elasticity - # if Set.InPlaneElasticity: - # gt, Kt, EBulk = KgBulk(geo_0, Geo, Set) - # K += Kt - # g += gt - # E += EBulk - # Energies["Bulk"] = EBulk - - # Bending Energy - # TODO - # Triangle Energy Barrier - if Set.EnergyBarrierA: - kg_Tri = KgTriEnergyBarrier(Geo) - kg_Tri.compute_work(Geo, Set, None, False) + if c_set.EnergyBarrierA: + kg_Tri = KgTriEnergyBarrier(geo) + kg_Tri.compute_work(geo, c_set, None, False) g += kg_Tri.g - Geo.create_vtk_cell(Set, num_step, 'Arrows_tri', -kg_Tri.g) + geo.create_vtk_cell(c_set, num_step, 'Arrows_tri', -kg_Tri.g) energies["TriEnergyBarrier"] = kg_Tri.energy # Triangle Energy Barrier Aspect Ratio - if Set.EnergyBarrierAR: - kg_TriAR = KgTriAREnergyBarrier(Geo) - kg_TriAR.compute_work(Geo, Set, None, False) + if c_set.EnergyBarrierAR: + kg_TriAR = KgTriAREnergyBarrier(geo) + kg_TriAR.compute_work(geo, c_set, None, False) g += kg_TriAR.g - Geo.create_vtk_cell(Set, num_step, 'Arrows_tri_ar', -kg_TriAR.g) + geo.create_vtk_cell(c_set, num_step, 'Arrows_tri_ar', -kg_TriAR.g) energies["TriEnergyBarrierAR"] = kg_TriAR.energy - # Propulsion Forces - # TODO - # Contractility - if Set.Contractility: - kg_lt = KgContractility(Geo) - kg_lt.compute_work(Geo, Set, None, False) + if c_set.Contractility: + kg_lt = KgContractility(geo) + kg_lt.compute_work(geo, c_set, None, False) g += kg_lt.g - Geo.create_vtk_cell(Set, num_step, 'Arrows_contractility', -kg_lt.g) + geo.create_vtk_cell(c_set, num_step, 'Arrows_contractility', -kg_lt.g) energies["Contractility"] = kg_lt.energy # Contractility as external force - if Set.Contractility_external: - kg_ext = KgContractilityExternal(Geo) - kg_ext.compute_work(Geo, Set, None, False) + if c_set.Contractility_external: + kg_ext = KgContractilityExternal(geo) + kg_ext.compute_work(geo, c_set, None, False) g += kg_ext.g - Geo.create_vtk_cell(Set, num_step, 'Arrows_contractility_external', -kg_ext.g) + geo.create_vtk_cell(c_set, num_step, 'Arrows_contractility_external', -kg_ext.g) energies["Contractility_external"] = kg_ext.energy # Substrate - if Set.Substrate == 2: - kg_subs = KgSubstrate(Geo) - kg_subs.compute_work(Geo, Set, None, False) + if c_set.Substrate == 2: + kg_subs = KgSubstrate(geo) + kg_subs.compute_work(geo, c_set, None, False) g += kg_subs.g - Geo.create_vtk_cell(Set, num_step, 'Arrows_substrate', -kg_subs.g) + geo.create_vtk_cell(c_set, num_step, 'Arrows_substrate', -kg_subs.g) energies["Substrate"] = kg_subs.energy - return g, energies \ No newline at end of file + return g, energies + +def check_if_fire_converged(geo, f_flat, c_set, dy_flat, v_flat, iteration_count): + """ + Realistic FIRE convergence for vertex models. + + Args: + f_flat: Force vector (negative gradient) + c_set: Simulation settings object + dy_flat: Displacement vector + v_flat: Velocity vector + iteration_count: Current FIRE iteration number + + Returns: + converged: Boolean indicating if minimization converged + reason: String explaining convergence reason + """ + # 1. Force-based convergence (most important) + max_force = np.max(np.abs(f_flat)) + force_tol = c_set.fire_force_tol + + # 2. Displacement-based (secondary) + max_disp = np.max(np.abs(dy_flat)) + disp_tol = c_set.fire_disp_tol + + # 3. Velocity-based (system at rest) + v_mag = np.linalg.norm(v_flat) + vel_tol = c_set.fire_vel_tol + + # 4. Maximum iterations + max_iter = c_set.fire_max_iterations + + converged = False + reason = "" + + # Check maximum iterations first + if iteration_count >= max_iter: + converged = True + reason = f"Reached max iterations ({max_iter})" + logger.warning(f"FIRE: {reason} - maxF={max_force:.3e}") + return converged, reason + + # Primary: Force convergence + if max_force < force_tol: + converged = True + reason = f"Max force {max_force:.2e} < {force_tol:.0e}" + + # Secondary: Small velocities AND small displacements + elif v_mag < vel_tol and max_disp < disp_tol: + converged = True + reason = f"System at rest: |v|={v_mag:.2e}, max disp={max_disp:.2e}" + + # Log progress every 10 iterations + if iteration_count % 10 == 0: + logger.info(f"FIRE iter {iteration_count}: maxF={max_force:.3e}, |v|={v_mag:.3e}, dt={geo._fire_dt:.4f}") + + return converged, reason + + +def single_iteration_fire(geo, c_set, dof, dy, g, selected_cells=None): + """ + CORRECTED FIRE algorithm with inner loop support. + + This function handles ONE FIRE iteration. For proper minimization, + it should be called repeatedly until convergence within a timestep. + + Algorithm (per iteration): + 1. Compute forces: F = -∇E + 2. MD integration: v = v + dt*F, x = x + dt*v + 3. Compute power: P = F·v + 4. If P > 0: apply velocity mixing, potentially accelerate + 5. If P <= 0: reset velocities, decrease dt + + This is the standard FIRE from Bitzek et al. (2006). + """ + # ============================================ + # INITIALIZATION & STATE MANAGEMENT + # ============================================ + + # Bottom nodes constraints + g_constrained = constrain_bottom_vertices_x_y(geo) + + # 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 + + # ============================================ + # FORCE COMPUTATION + # ============================================ + + # Forces are negative gradient (F = -∇E = -g) + F = -g.copy() + + # Extract free DOF velocities and forces + v_flat = geo._fire_velocity.flatten()[dof] + F_flat = F[dof] + + # ============================================ + # STEP 1: MD INTEGRATION + # ============================================ + + # Velocity update (Euler): v_new = v_old + dt*F + v_flat = v_flat + geo._fire_dt * F_flat + + # Position update (Euler): x_new = x_old + dt*v_new + dy_flat = geo._fire_dt * v_flat + + # ============================================ + # STEP 2: FIRE ADAPTATION + # ============================================ + + # Compute power P = F·v (using velocities AFTER integration) + P = np.dot(F_flat, v_flat) + + if P > 1e-16: # GOOD: moving downhill (add small tolerance) + geo._fire_n_positive += 1 + + # Apply FIRE velocity mixing: v = (1-α)v + α*|v|*(F/|F|) + v_mag = np.linalg.norm(v_flat) + F_mag = np.linalg.norm(F_flat) + + if v_mag > 1e-10 and F_mag > 1e-10: + F_hat = F_flat / F_mag + v_flat = (1.0 - geo._fire_alpha) * v_flat + geo._fire_alpha * v_mag * F_hat + + # Accelerate if we've had N_min consecutive good steps + if geo._fire_n_positive > c_set.fire_N_min: + # Increase timestep (up to maximum) + if geo._fire_dt * c_set.fire_f_inc < c_set.fire_dt_max: + geo._fire_dt = geo._fire_dt * c_set.fire_f_inc + else: + geo._fire_dt = c_set.fire_dt_max + # Decrease mixing parameter + geo._fire_alpha = geo._fire_alpha * c_set.fire_f_alpha + + else: # BAD: overshoot or oscillation + logger.debug(f"FIRE reset: P={P:.3e} <= 0") + geo._fire_n_positive = 0 + # RESET velocities to zero (critical!) + v_flat = np.zeros_like(v_flat) + # Decrease timestep + geo._fire_dt = max(geo._fire_dt * c_set.fire_f_dec, c_set.fire_dt_min) + # Reset mixing parameter + geo._fire_alpha = c_set.fire_alpha_start + + # ============================================ + # STEP 3: UPDATE GEOMETRY + # ============================================ + + # Store updated velocity + geo._fire_velocity.flatten()[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 + + # Build full displacement vector + dy[dof, 0] = dy_flat + dy = map_vertices_periodic_boundaries(geo, dy) + dyr = np.linalg.norm(dy[dof, 0]) + dy_reshaped = np.reshape(dy, (geo.numF + geo.numY + geo.nCells, 3)) + + # Update geometry + geo.update_vertices(dy_reshaped, selected_cells) + if c_set.frozen_face_centres or c_set.frozen_face_centres_border_cells: + for cell in geo.Cells: + if cell.AliveStatus is not None and ( + (cell.ID in geo.BorderCells and c_set.frozen_face_centres_border_cells) + or c_set.frozen_face_centres + ): + face_centres_to_middle_of_neighbours_vertices(geo, cell.ID) + + geo.update_measures() + c_set.iter = c_set.MaxIter + + # ============================================ + # CONVERGENCE CHECK + # ============================================ + + converged, reason = check_if_fire_converged( + geo, F_flat, c_set, dy_flat, v_flat, geo._fire_iteration_count + ) + + logger.debug(f"FIRE: dt={geo._fire_dt:.4f}, α={geo._fire_alpha:.4f}, " + f"P={P:.3e}, |v|={np.linalg.norm(v_flat):.3e}, " + f"iter={geo._fire_iteration_count}") + + if converged: + logger.info(f"FIRE converged: {reason}") + + return geo, dy, dyr, converged + + +def fire_minimization_loop(geo, c_set, dof, g, t, num_step, selected_cells=None): + """ + Complete FIRE minimization loop for one timestep. + + This function runs FIRE repeatedly until convergence or max iterations. + It should be called ONCE per simulation timestep. + + Returns: + Geo: Minimized geometry + converged: Whether minimization converged + iterations: Number of FIRE iterations performed + final_gradient_norm: Norm of gradient at convergence + """ + 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) + + # Store initial gradient for reference + initial_gradient_norm = np.linalg.norm(g[dof]) + logger.info(f"Initial gradient norm: {initial_gradient_norm:.3e}") + + # Main minimization loop + converged = False + fire_iterations = 0 + max_iterations = c_set.fire_max_iterations + + while not converged and fire_iterations < max_iterations: + # One FIRE iteration + geo, dy, dyr, converged = single_iteration_fire( + geo, c_set, dof, dy, g, selected_cells + ) + + fire_iterations += 1 + + # Recompute gradient for next iteration + # (This depends on your gradient computation function) + g, energies = g_global(geo, c_set, None, False, num_step) + + # Optional: Break if stuck + if fire_iterations > 50 and geo._fire_n_positive == 0: + logger.warning("FIRE: No positive P steps for 50 iterations, likely stuck") + break + + # Check for NaNs as a safety + if np.any(np.isnan(g[dof])) or np.any(np.isnan(dy[dof])): + logger.error("FIRE: NaN detected in gradient or displacement!") + converged = False + + # Final statistics + final_gradient_norm = np.linalg.norm(g[dof]) if hasattr(geo, '_fire_velocity') else float('inf') + + if converged: + logger.info(f"FIRE minimization successful: {fire_iterations} iterations, " + f"gradient reduced from {initial_gradient_norm:.3e} to {final_gradient_norm:.3e}") + else: + 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 + + return geo, converged, final_gradient_norm + + +def initialize_fire(geo, c_set): + """ + Initialize FIRE algorithm state variables on the Geo object. + :param geo: + :param c_set: + :return: + """ + geo._fire_dt = c_set.dt + geo._fire_alpha = c_set.fire_alpha_start + geo._fire_n_positive = 0 + geo._fire_iteration_count = 0 \ No newline at end of file diff --git a/src/pyVertexModel/algorithm/vertexModel.py b/src/pyVertexModel/algorithm/vertexModel.py index 79d8e1c3..b0c77901 100644 --- a/src/pyVertexModel/algorithm/vertexModel.py +++ b/src/pyVertexModel/algorithm/vertexModel.py @@ -15,7 +15,7 @@ from scipy.optimize import minimize from skimage.measure import regionprops -from pyVertexModel.algorithm import newtonRaphson +from pyVertexModel.algorithm import integrators from pyVertexModel.geometry import degreesOfFreedom from pyVertexModel.geometry.geo import ( Geo, @@ -495,7 +495,7 @@ def iterate_over_time(self): gr = self.single_iteration() if np.isnan(gr): - break + self.didNotConverge = True self.save_v_model_state() @@ -527,22 +527,28 @@ def single_iteration(self, post_operations=True): # up-to-date self.geo.update_measures() - if self.set.implicit_method is True: - g, K, _, energies = newtonRaphson.KgGlobal(self.geo_0, self.geo_n, self.geo, self.set, - self.set.implicit_method) + if self.set.implicit_method: + g, K, _, energies = integrators.KgGlobal(self.geo, self.set, self.geo_n, self.set.implicit_method) else: K = 0 - g, energies = newtonRaphson.gGlobal(self.geo_0, self.geo_n, self.geo, self.set, - self.set.implicit_method, self.numStep) + g, energies = integrators.g_global(self.geo, self.set, self.geo_n, self.set.implicit_method, self.numStep) for key, energy in energies.items(): logger.info(f"{key}: {energy}") - self.geo, g, __, __, self.set, gr, dyr, dy = newtonRaphson.newton_raphson(self.geo_0, self.geo_n, self.geo, - self.Dofs, self.set, K, g, - self.numStep, self.t, - self.set.implicit_method) - if not np.isnan(gr) and post_operations: - self.post_newton_raphson(dy, g, gr) + + if self.set.integrator == 'fire': + self.geo, converged, gr = integrators.fire_minimization_loop(self.geo, self.set, self.Dofs.Free, g, self.t, + self.numStep) + if converged: + self.iteration_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, + self.numStep, self.t, + self.set.implicit_method) + if not np.isnan(gr) and post_operations: + self.post_newton_raphson(dy, g, gr) + return gr def post_newton_raphson(self, dy, g, gr): @@ -553,9 +559,19 @@ def post_newton_raphson(self, dy, g, gr): :param gr: :return: """ - if ((gr * self.set.dt / self.set.dt0) < self.set.tol and np.all(~np.isnan(g[self.Dofs.Free])) and - np.all(~np.isnan(dy[self.Dofs.Free])) and - (np.max(abs(g[self.Dofs.Free])) * self.set.dt / self.set.dt0) < self.set.tol): + # Standard convergence criterion + if self.set.implicit_method: + converged = ((gr * self.set.dt / self.set.dt0) < self.set.tol and np.all(~np.isnan(g[self.Dofs.Free])) and + np.all(~np.isnan(dy[self.Dofs.Free])) and + (np.max(abs(g[self.Dofs.Free])) * self.set.dt / self.set.dt0) < self.set.tol) + else: + # For explicit methods: Don't reject based on gradient increase + # The adaptive scaling in newton_raphson_iteration_explicit handles stability + # Rejection causes geometry restore and dt reduction, making convergence harder + converged = (np.all(~np.isnan(g[self.Dofs.Free])) and np.all(~np.isnan(dy[self.Dofs.Free])) and + (np.max(abs(g[self.Dofs.Free])) < self.set.tol or self.set.dt <= self.set.dt0 * self.set.dt_tolerance)) + + if converged: self.iteration_converged() else: self.iteration_did_not_converged() @@ -567,20 +583,23 @@ def iteration_did_not_converged(self): If the iteration did not converge, the algorithm will try to relax the value of nu and dt. :return: """ - self.geo, self.geo_n, self.geo_0, self.tr, self.Dofs = load_backup_vars(self.backupVars) - self.relaxingNu = False - if self.set.iter == self.set.MaxIter0 and self.set.implicit_method: - self.set.MaxIter = self.set.MaxIter0 * 3 - self.set.nu = 10 * self.set.nu0 - else: - if (self.set.iter >= self.set.MaxIter and - (self.set.dt / self.set.dt0) > 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 + if self.set.integrator != 'fire': + self.geo, self.geo_n, self.geo_0, self.tr, self.Dofs = load_backup_vars(self.backupVars) + self.relaxingNu = False + if self.set.iter == self.set.MaxIter0 and self.set.implicit_method: + self.set.MaxIter = self.set.MaxIter0 * 3 + self.set.nu = 10 * self.set.nu0 else: - self.didNotConverge = True + if (self.set.iter >= self.set.MaxIter and + (self.set.dt / self.set.dt0) > self.set.dt_tolerance): + 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") def iteration_converged(self): """ @@ -1158,6 +1177,9 @@ def required_purse_string_strength(self, directory, tend=20.1, load_existing=Tru else: run_iteration = True self.set.OutputFolder = output_directory + self.set.integrator = 'euler' + self.set.dt_tolerance = 1e-11 + if (self.set.dt / self.set.dt0) <= 1e-6: return np.inf, np.inf, np.inf, np.inf @@ -1190,6 +1212,9 @@ def required_purse_string_strength(self, directory, tend=20.1, load_existing=Tru for sub_f in os.listdir(sub_dir): shutil.copy(os.path.join(sub_dir, sub_f), os.path.join(dest_sub_dir, sub_f)) + if self.t < tend: + print(f'Could not reach the end time of the simulation for purse string strength exploration. Current time: {self.t}') + return np.inf, np.inf, np.inf, np.inf dy_values, purse_string_strength_values = self.required_purse_string_strength_for_timepoint(directory, timepoint=tend) purse_string_strength_values = purse_string_strength_values[:len(dy_values)] diff --git a/src/pyVertexModel/analysis/analyse_simulation.py b/src/pyVertexModel/analysis/analyse_simulation.py index 94985b18..1dd2697a 100644 --- a/src/pyVertexModel/analysis/analyse_simulation.py +++ b/src/pyVertexModel/analysis/analyse_simulation.py @@ -12,7 +12,6 @@ load_state, load_variables, save_variables, - screenshot, ) @@ -418,6 +417,7 @@ 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/analyse_simulations.py b/src/pyVertexModel/analysis/analyse_simulations.py index e0aed197..15f5ab32 100644 --- a/src/pyVertexModel/analysis/analyse_simulations.py +++ b/src/pyVertexModel/analysis/analyse_simulations.py @@ -3,7 +3,7 @@ import numpy as np import pandas as pd -from pyVertexModel.analysis.analyse_simulation import analyse_simulation, create_video +from pyVertexModel.analysis.analyse_simulation import analyse_simulation from pyVertexModel.util.utils import plot_figure_with_line folder = '/media/pablo/d7c61090-024c-469a-930c-f5ada47fb049/PabloVicenteMunuera/VertexModel/pyVertexModel/Result/same_recoil_wound_area_results/c/' diff --git a/src/pyVertexModel/analysis/find_required_purse_string.py b/src/pyVertexModel/analysis/find_required_purse_string.py index 6996e6ba..abcb5a9c 100644 --- a/src/pyVertexModel/analysis/find_required_purse_string.py +++ b/src/pyVertexModel/analysis/find_required_purse_string.py @@ -2,6 +2,7 @@ import os import sys +import numpy as np import pandas as pd from pyVertexModel import PROJECT_DIRECTORY @@ -14,8 +15,8 @@ plot_figures = True # Folder containing different simulations with different cell shapes -c_folder = os.path.join(PROJECT_DIRECTORY, 'Result/to_calculate_ps_recoil/') -output_csv = os.path.join(PROJECT_DIRECTORY, 'Result/to_calculate_ps_recoil/required_purse_string_strengths.csv') +c_folder = os.path.join(PROJECT_DIRECTORY, 'Result/to_calculate_ps_recoil/') #same_recoil_wound_area_results #to_calculate_ps_recoil +output_csv = os.path.join(c_folder, 'required_purse_string_strengths.csv') if plot_figures and os.path.exists(output_csv): # Load the existing csv file @@ -38,45 +39,71 @@ for timepoint in [0.1, 6.0]: df_time = df[df['time'] == timepoint] - # Keep recoil values that are close to 2.4e-13 +- 4e-14 - # if timepoint == 0.1: - # df_time = df_time[(df_time['Recoil'] > 2.0e-13) & (df_time['Recoil'] < 3e-13)] + if df_time.empty: + print(f"No data found for timepoint {timepoint}. Skipping.") + continue + + if c_folder.__contains__('same_recoil_wound_area_results'): + input_excel_file = 'all_files_features.xlsx' + else: + input_excel_file = 'all_files_features_same_area_ablating.xlsx' + + # Keep only the model and aspect ratio in excel file of 'all_files_features_same_area_ablating' + df_filtered = pd.read_excel(os.path.join(c_folder, input_excel_file), header=0) + df_filtered = df_filtered[df_filtered['last_area_time_top'] > 20.0] + df_filtered['Model_name'] = df_filtered['model'] + '_' + df_filtered['AR'].astype(str) + + # Keep only the df_time rows with model names that are in df_filtered + df_time = df_time[df_time['Model_name'].isin(df_filtered['Model_name'])] + df_time = df_time.reset_index() # Save the cleaned DataFrame with 'filtered' suffix - df_time.to_csv(os.path.join(PROJECT_DIRECTORY, 'Result', 'to_calculate_ps_recoil', 'required_purse_string_strengths_filtered' + f'_time_{timepoint}.csv'), index=False) + df_time.to_csv(os.path.join(c_folder, 'required_purse_string_strengths_filtered' + f'_time_{timepoint}.csv'), index=False) + + # Create a folder for the timepoint if it doesn't exist + timepoint_folder = os.path.join(c_folder, str(timepoint)) + if not os.path.exists(timepoint_folder): + os.makedirs(timepoint_folder) + + # # Normalise purse string strength by its resting line tension + # df_edge_tension = pd.read_excel(os.path.join(c_folder, 'c/all_simulations_metrics.xlsx'), header=0) + # df_edge_tension['Model_name'] = df_edge_tension['model_name'] + '_' + df_edge_tension['AR'].astype(str) + # df_time = df_time.merge(df_edge_tension[['Model_name', 'final_edge_normalised']], on='Model_name', how='left') + # + # df_time['Purse_string_strength'] = df_time['Purse_string_strength'] / df_time['final_edge_normalised'] # Plot recoil over aspect ratio - plot_figure_with_line(df_time, None, os.path.join(PROJECT_DIRECTORY, 'Result', 'to_calculate_ps_recoil', str(timepoint)), + plot_figure_with_line(df_time, None, os.path.join(c_folder, str(timepoint)), x_axis_name='AR', y_axis_name='Recoil', y_axis_label='Recoil velocity (t=' + str(timepoint) + ')') - plot_figure_with_line(df_time, None, os.path.join(PROJECT_DIRECTORY, 'Result', 'to_calculate_ps_recoil', str(timepoint)), + plot_figure_with_line(df_time, None, os.path.join(c_folder, str(timepoint)), x_axis_name='AR', y_axis_name='Purse_string_strength', y_axis_label='Purse string strength (t=' + str(timepoint) + ')') - plot_figure_with_line(df_time, None, os.path.join(PROJECT_DIRECTORY, 'Result', 'to_calculate_ps_recoil', str(timepoint)), + plot_figure_with_line(df_time, None, os.path.join(c_folder, str(timepoint)), x_axis_name='AR', y_axis_name='LambdaS1', y_axis_label=r'$\lambda_{s1}=\lambda_{s3}$') - plot_figure_with_line(df_time, None, os.path.join(PROJECT_DIRECTORY, 'Result', 'to_calculate_ps_recoil', str(timepoint)), + plot_figure_with_line(df_time, None, os.path.join(c_folder, str(timepoint)), x_axis_name='AR', y_axis_name='LambdaS2', y_axis_label=r'$\lambda_{s2}$') - output_csv = output_csv.replace('required_purse_string_strengths.csv', 'all_files_features_filtered.xlsx') - df = pd.read_excel(output_csv) - - plot_figure_with_line(df, None, os.path.join(PROJECT_DIRECTORY, 'Result', 'to_calculate_ps_recoil'), - x_axis_name='AR', - y_axis_name='wound_area_top_extrapolated_60', y_axis_label='Wound area top at 60 min.') - plot_figure_with_line(df, None, os.path.join(PROJECT_DIRECTORY, 'Result', 'to_calculate_ps_recoil'), - x_axis_name='AR', - y_axis_name='wound_area_bottom_extrapolated_60', y_axis_label='Wound area bottom at 60 min.') - plot_figure_with_line(df, None, os.path.join(PROJECT_DIRECTORY, 'Result', 'to_calculate_ps_recoil'), - x_axis_name='AR', - y_axis_name='lS1', y_axis_label=r'$\lambda_{s1}=\lambda_{s3}$') + # output_csv = output_csv.replace('required_purse_string_strengths.csv', 'all_files_features_filtered.xlsx') + # df = pd.read_excel(output_csv) + # + # plot_figure_with_line(df, None, os.path.join(c_folder), + # x_axis_name='AR', + # y_axis_name='wound_area_top_extrapolated_60', y_axis_label='Wound area top at 60 min.') + # plot_figure_with_line(df, None, os.path.join(c_folder), + # x_axis_name='AR', + # y_axis_name='wound_area_bottom_extrapolated_60', y_axis_label='Wound area bottom at 60 min.') + # plot_figure_with_line(df, None, os.path.join(c_folder), + # x_axis_name='AR', + # y_axis_name='lS1', y_axis_label=r'$\lambda_{s1}=\lambda_{s3}$') else: #if not os.path.exists(output_csv): # Get all directories within c_folder @@ -101,11 +128,21 @@ simulations_dirs.sort(reverse=True) directory = simulations_dirs[int(sys.argv[1])] - if directory.startswith('no_'): - continue + #if not directory.startswith('no_Remodelling_ablating_'): + # continue #for directory in simulations_dirs: print(f"Processing directory: {directory}") + # 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_')] + if len(dirs_within) == 0: + print(f"No directories starting with 'no_Remodelling_ablating_' found in {directory}, skipping.") + continue + + directory = dirs_within[0] + print(f"Processing directory: {directory}") + # Get the purse string strength vModel = VertexModelVoronoiFromTimeImage(create_output_folder=False, set_option='wing_disc') @@ -116,13 +153,15 @@ 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 # Run the required purse string strength analysis current_folder = vModel.set.OutputFolder last_folder = current_folder.split('/') vModel.set.OutputFolder = os.path.join(PROJECT_DIRECTORY, 'Result/', last_folder[-1]) _, _, recoiling, purse_string_strength_eq = vModel.required_purse_string_strength( - os.path.join(c_folder, ar_dir, directory), tend=t_ablation + 0.1, load_existing=False) + os.path.join(c_folder, ar_dir, directory), tend=t_ablation + 0.1, load_existing=True) recoilings.append(recoiling) purse_string_strength_0.append(purse_string_strength_eq) @@ -132,29 +171,29 @@ model_name.append(vModel.set.model_name) time_list.append(0.1) - vModel.set.OutputFolder = os.path.join(PROJECT_DIRECTORY, 'Result/', last_folder[-1]) - _, _, recoiling_t_6, purse_string_strength_eq_t_6 = vModel.required_purse_string_strength( - os.path.join(c_folder, ar_dir, directory), tend=t_ablation + 6.0) - - recoilings.append(recoiling_t_6) - purse_string_strength_0.append(purse_string_strength_eq_t_6) - aspect_ratio.append(vModel.set.CellHeight) - lambda_s1_list.append(vModel.set.lambdaS1) - lambda_s2_list.append(vModel.set.lambdaS2) - model_name.append(vModel.set.model_name) - time_list.append(6.0) - analyse_simulation(os.path.join(c_folder, ar_dir, directory)) - - # Append results into an existing csv file - if os.path.exists(output_csv): - df_existing = pd.read_csv(output_csv) - df = pd.DataFrame( - {'Model_name': model_name, 'AR': aspect_ratio, 'LambdaS1': lambda_s1_list, 'LambdaS2': lambda_s2_list, - 'Recoil': recoilings, 'Purse_string_strength': purse_string_strength_0, 'time': time_list}) - df_final = pd.concat([df_existing, df], ignore_index=True) - df_final.to_csv(output_csv, index=False) - else: - df = pd.DataFrame( - {'Model_name': model_name, 'AR': aspect_ratio, 'LambdaS1': lambda_s1_list, 'LambdaS2': lambda_s2_list, - 'Recoil': recoilings, 'Purse_string_strength': purse_string_strength_0, 'time': time_list}) - df.to_csv(output_csv, index=False) \ No newline at end of file + # vModel.set.OutputFolder = os.path.join(PROJECT_DIRECTORY, 'Result/', last_folder[-1]) + # _, _, recoiling_t_6, purse_string_strength_eq_t_6 = vModel.required_purse_string_strength( + # os.path.join(c_folder, ar_dir, directory), tend=t_ablation + 6.0) + # + # recoilings.append(recoiling_t_6) + # purse_string_strength_0.append(purse_string_strength_eq_t_6) + # aspect_ratio.append(vModel.set.CellHeight) + # lambda_s1_list.append(vModel.set.lambdaS1) + # lambda_s2_list.append(vModel.set.lambdaS2) + # model_name.append(vModel.set.model_name) + # time_list.append(6.0) + # analyse_simulation(os.path.join(c_folder, ar_dir, directory)) + + #Append results into an existing csv file + if os.path.exists(output_csv): + df_existing = pd.read_csv(output_csv) + df = pd.DataFrame( + {'Model_name': model_name, 'AR': aspect_ratio, 'LambdaS1': lambda_s1_list, 'LambdaS2': lambda_s2_list, + 'Recoil': recoilings, 'Purse_string_strength': purse_string_strength_0, 'time': time_list}) + df_final = pd.concat([df_existing, df], ignore_index=True) + df_final.to_csv(output_csv, index=False) + else: + df = pd.DataFrame( + {'Model_name': model_name, 'AR': aspect_ratio, 'LambdaS1': lambda_s1_list, 'LambdaS2': lambda_s2_list, + 'Recoil': recoilings, 'Purse_string_strength': purse_string_strength_0, 'time': time_list}) + df.to_csv(output_csv, index=False) \ No newline at end of file diff --git a/src/pyVertexModel/geometry/geo.py b/src/pyVertexModel/geometry/geo.py index 6c458965..815ecd69 100644 --- a/src/pyVertexModel/geometry/geo.py +++ b/src/pyVertexModel/geometry/geo.py @@ -2274,4 +2274,174 @@ def resize_tissue(self, average_volume=0.0003168604676977124): self.update_measures() volumes_after_deformation = np.array([cell.Vol for cell in self.Cells if cell.AliveStatus is not None]) + + def geometry_is_correct(self): + """ + Check if the geometry is structurally valid using vertex model expert criteria. + + This function validates the vertex model geometry structure by checking: + 1. Cell volumes (positive, not degenerate) + 2. Cell Y coordinates (no None, NaN, or Inf) + 3. Face areas (positive, not too small) + 4. Triangle degeneracy (no collapsed triangles) + 5. Face planarity (critical for detecting "spiky cells") + + Dead cells (AliveStatus=None) are ignored as they may legitimately have Y=None. + + Returns: + bool: True if geometry is structurally valid, False if invalid/going wild + """ + + # Get alive cells only + alive_cells = [c for c in self.Cells if hasattr(c, 'AliveStatus') and c.AliveStatus is not None] + + if not alive_cells: + logger.debug("No alive cells found") + 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") + 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") + 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") + return False + + if cell.Vol <= 0: + logger.debug(f"Cell {cell_idx} has non-positive volume: {cell.Vol}") + return False + + if cell.Vol < 1e-10: + logger.debug(f"Cell {cell_idx} has tiny volume: {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'): + 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}") + return False + + if face.Area < 1e-12: + logger.debug(f"Cell {cell_idx}, Face {face_idx} has degenerate area: {face.Area}") + return False + + # 4. Check triangle areas and degeneracy + for cell_idx, cell in enumerate(alive_cells): + if not hasattr(cell, 'Faces'): + continue + + for face_idx, face in enumerate(cell.Faces): + if not hasattr(face, 'Tris'): + continue + + for tri_idx, tri in enumerate(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)") + return False + + # Check triangle area + if hasattr(tri, 'Area') and tri.Area is not None: + if tri.Area < 0: + logger.debug(f"Cell {cell_idx}, Face {face_idx}, 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'): + continue + + for face_idx, face in enumerate(cell.Faces): + if not hasattr(face, 'Tris'): + continue + + for tri in face.Tris: + total_triangles += 1 + + # Count tiny triangles + if hasattr(tri, 'Area') and tri.Area is not None and tri.Area < 1e-10: + tiny_triangle_count += 1 + + # If more than 0.2% of triangles are tiny, geometry is problematic + # Empirically: correct geometries have ~0.05% tiny triangles, going_wild_1 has ~0.29% + if total_triangles > 0: + tiny_ratio = tiny_triangle_count / total_triangles + if tiny_ratio > 0.002: # 0.2% threshold + logger.debug(f"Too many tiny triangles: {tiny_triangle_count}/{total_triangles} ({tiny_ratio*100:.2f}%)") + return False + + # 5. Check face planarity (CRITICAL for detecting "spiky cells") + # In a proper vertex model, faces should be approximately planar. + # 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'): + continue + + for face_idx, face in enumerate(cell.Faces): + if not hasattr(face, 'Tris'): + continue + + # Get all unique vertices in the face + vertices = set() + for tri in 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]) + + # 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})") + + except (IndexError, ValueError) as e: + logger.debug(f"Cell {cell_idx}, Face {face_idx}: Error checking planarity: {e}") + return False + + # Any non-planar faces indicate geometry is going wild + if non_planar_count > 0: + logger.debug(f"Found {non_planar_count} non-planar faces") + return False + + # All checks passed + return True logger.info(f'Volume difference: {np.mean(volumes_after_deformation) - average_volume}') \ No newline at end of file diff --git a/src/pyVertexModel/mesh_remodelling/remodelling.py b/src/pyVertexModel/mesh_remodelling/remodelling.py index 9e8f7cdb..97de0623 100644 --- a/src/pyVertexModel/mesh_remodelling/remodelling.py +++ b/src/pyVertexModel/mesh_remodelling/remodelling.py @@ -3,9 +3,9 @@ import numpy as np import pandas as pd -from pyVertexModel.algorithm.newtonRaphson import ( +from pyVertexModel.algorithm.integrators import ( constrain_bottom_vertices_x_y, - gGlobal, + g_global, newton_raphson_iteration_explicit, ) from pyVertexModel.geometry.geo import ( @@ -323,7 +323,7 @@ def remodel_mesh(self, num_step, remodelled_cells, how_close_to_vertex=0.01): backup_vars = save_backup_vars(self.Geo, self.Geo_n, self.Geo_0, num_step, self.Dofs) # self.Geo.create_vtk_cell(self.Geo_0, self.Set, num_step) - g, energies = gGlobal(self.Geo, self.Geo, self.Geo, self.Set, self.Set.implicit_method) + g, energies = g_global(self.Geo, self.Set, self.Geo, self.Set.implicit_method) gr = np.linalg.norm(g[self.Dofs.Free]) logger.info(f'|gr| before remodelling: {gr}') for key, energy in energies.items(): @@ -506,14 +506,14 @@ def check_if_will_converge(self, best_geo, n_iter_max=20): # Check if the remodelling will converge dy = np.zeros(((best_geo.numY + best_geo.numF + best_geo.nCells) * 3, 1), dtype=np.float64) - g, energies = gGlobal(best_geo, best_geo, best_geo, self.Set, self.Set.implicit_method) + g, energies = g_global(best_geo, self.Set, best_geo, self.Set.implicit_method) previous_gr = np.linalg.norm(g[self.Dofs.Free]) try: for n_iter in range(n_iter_max): best_geo, dy, dyr = newton_raphson_iteration_explicit(best_geo, self.Set, self.Dofs.Free, dy, g) - g, energies = gGlobal(best_geo, best_geo, best_geo, self.Set, self.Set.implicit_method) + g, energies = g_global(best_geo, self.Set, best_geo, self.Set.implicit_method) for key, energy in energies.items(): logger.info(f"{key}: {energy}") g_constrained = constrain_bottom_vertices_x_y(best_geo) @@ -563,7 +563,7 @@ def intercalate_cells(self, segmentFeatures): # Check if the remodelling has improved the gr and the energy # Compute the new energy - g, energies = gGlobal(self.Geo, self.Geo, self.Geo, self.Set, self.Set.implicit_method) + g, energies = g_global(self.Geo, self.Set, self.Geo, self.Set.implicit_method) self.Dofs.get_dofs(self.Geo, self.Set) gr = np.linalg.norm(g[self.Dofs.Free]) logger.info(f'|gr| after remodelling without changes: {gr}') diff --git a/src/pyVertexModel/parameters/set.py b/src/pyVertexModel/parameters/set.py index 3bb4fb61..bda8ad12 100644 --- a/src/pyVertexModel/parameters/set.py +++ b/src/pyVertexModel/parameters/set.py @@ -21,6 +21,7 @@ def __init__(self, mat_file=None): Parameters: mat_file (optional): A MATLAB-style struct or object containing saved Set parameters; when provided, values are read and assigned to this instance. """ + self.dt_tolerance = 1e-6 self.min_3d_neighbours = None self.periodic_boundaries = True self.frozen_face_centres = False @@ -44,6 +45,28 @@ def __init__(self, mat_file=None): self.dt = None self.dt0 = None self.implicit_method = False + self.integrator = 'euler' # Time integrator: 'euler', 'rk2' (midpoint method), or 'fire' (FIRE algorithm) + + # FIRE Algorithm parameters (Bitzek et al., 2006) + # These parameters control the adaptive optimization when integrator='fire' + self.fire_dt_max = None # Maximum timestep (will be set to 10*dt if None) + self.fire_dt_min = None # Minimum timestep (will be set to 0.02*dt if None) + self.fire_N_min = 5 # Steps before acceleration (recommended: 5) + self.fire_f_inc = 1.1 # dt increase factor (recommended: 1.1) + self.fire_f_dec = 0.5 # dt decrease factor (recommended: 0.5) + self.fire_alpha_start = 0.1 # Initial damping coefficient (recommended: 0.1) + self.fire_f_alpha = 0.99 # α decrease factor (recommended: 0.99) + + self.fire_dt_max = 0.5 # Large max dt for fast minimization + self.fire_dt_min = 1e-8 # Very small min dt + + # Convergence tolerances - PRACTICAL FOR VERTEX MODELS + self.fire_force_tol = 1e-6 # Tight for steady-state + self.fire_disp_tol = 1e-8 # Tight displacement + self.fire_vel_tol = 1e-10 # Tight velocity + self.fire_max_iterations = 100 # Allow more iterations for tight convergence + + # Additional parameters self.TypeOfPurseString = None self.Contractility_TimeVariability = None self.Contractility_Variability_LateralCables = None @@ -198,9 +221,8 @@ def check_for_non_used_parameters(self): if not self.brownian_motion: self.brownian_motion_scale = 0 - if self.implicit_method is False: + if not self.implicit_method: self.tol = self.nu - self.tol0 = self.nu/20 if self.Remodelling: self.RemodelStiffness = 0.9 @@ -341,7 +363,7 @@ def bubbles(self): self.check_for_non_used_parameters() def cyst(self): - mat_info = scipy.io.loadmat(os.path.join(PROJECT_DIRECTORY, 'Tests/data/Geo_var_cyst.mat')) + mat_info = scipy.io.loadmat(os.path.join(PROJECT_DIRECTORY, 'Tests_data/Geo_var_cyst.mat')) self.read_mat_file(mat_info['Set']) self.InputGeo = 'Bubbles_Cyst' self.CellHeight = 15 @@ -369,13 +391,14 @@ def wing_disc_apical_constriction(self): self.check_for_non_used_parameters() def wing_disc_equilibrium(self): + self.integrator = 'euler' self.nu_bottom = self.nu self.model_name = 'dWL1' self.initial_filename_state = 'Input/images/' + self.model_name + '.tif' #self.initial_filename_state = 'Input/images/dWP1_150cells_15_scutoids_1.0.pkl' self.percentage_scutoids = 0.0 self.tend = 20 - self.Nincr = self.tend * 100 + self.Nincr = self.tend * 10 self.CellHeight = 1 * 15 #np.array([0.0001, 0.001, 0.01, 0.1, 0.5, 1, 2.0]) * original_wing_disc_height #self.resize_z = 0.01 self.min_3d_neighbours = None # 10 diff --git a/src/pyVertexModel/util/utils.py b/src/pyVertexModel/util/utils.py index 5cbfa7f2..853100b4 100644 --- a/src/pyVertexModel/util/utils.py +++ b/src/pyVertexModel/util/utils.py @@ -688,6 +688,18 @@ def plot_figure_with_line(best_average_values, scutoids, current_path, x_axis_na plt.xticks(fontsize=20, fontweight='bold') plt.yticks(fontsize=20, fontweight='bold') + # Make yticks in scientific notation + plt.ticklabel_format(style='sci', axis='y', scilimits=(0, 0)) + + # From 0 ylim always + plt.ylim(0, None) + + if y_axis_label == 'Purse string strength (t=' + str(0.1) + ')': + plt.ylim(0, 1.2e-3) + + if y_axis_name == 'top_closure_velocity': + plt.ylim(0, 4.0) + # plt.title(f'Boxplot of {param} correlations with {scutoids*100}% scutoids') plt.xlabel(x_axis_label, fontsize=20, fontweight='bold') plt.ylabel(y_axis_label, fontsize=20, fontweight='bold') diff --git a/uv.lock b/uv.lock index db6673da..29908229 100644 --- a/uv.lock +++ b/uv.lock @@ -316,6 +316,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/ff/fa/d3c15189f7c52aaefbaea76fb012119b04b9013f4bf446cb4eb4c26c4e6b/cython-3.2.4-py3-none-any.whl", hash = "sha256:732fc93bc33ae4b14f6afaca663b916c2fdd5dcbfad7114e17fb2434eeaea45c", size = 1257078, upload-time = "2026-01-04T14:14:12.373Z" }, ] +[[package]] +name = "et-xmlfile" +version = "2.0.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/d3/38/af70d7ab1ae9d4da450eeec1fa3918940a5fafb9055e934af8d6eb0c2313/et_xmlfile-2.0.0.tar.gz", hash = "sha256:dab3f4764309081ce75662649be815c4c9081e88f0837825f90fd28317d4da54", size = 17234, upload-time = "2024-10-25T17:25:40.039Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/c1/8b/5fe2cc11fee489817272089c4203e679c63b570a5aaeb18d852ae3cbba6a/et_xmlfile-2.0.0-py3-none-any.whl", hash = "sha256:7a91720bc756843502c3b7504c77b8fe44217c85c537d85037f0f536151b2caa", size = 18059, upload-time = "2024-10-25T17:25:39.051Z" }, +] + [[package]] name = "fonttools" version = "4.61.1" @@ -781,6 +790,37 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/5b/c7/b801bf98514b6ae6475e941ac05c58e6411dd863ea92916bfd6d510b08c1/numpy-2.4.1-pp311-pypy311_pp73-win_amd64.whl", hash = "sha256:4f1b68ff47680c2925f8063402a693ede215f0257f02596b1318ecdfb1d79e33", size = 12492579, upload-time = "2026-01-10T06:44:57.094Z" }, ] +[[package]] +name = "opencv-python" +version = "4.13.0.92" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "numpy", version = "2.4.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, +] +wheels = [ + { url = "https://files.pythonhosted.org/packages/fc/6f/5a28fef4c4a382be06afe3938c64cc168223016fa520c5abaf37e8862aa5/opencv_python-4.13.0.92-cp37-abi3-macosx_13_0_arm64.whl", hash = "sha256:caf60c071ec391ba51ed00a4a920f996d0b64e3e46068aac1f646b5de0326a19", size = 46247052, upload-time = "2026-02-05T07:01:25.046Z" }, + { url = "https://files.pythonhosted.org/packages/08/ac/6c98c44c650b8114a0fb901691351cfb3956d502e8e9b5cd27f4ee7fbf2f/opencv_python-4.13.0.92-cp37-abi3-macosx_14_0_x86_64.whl", hash = "sha256:5868a8c028a0b37561579bfb8ac1875babdc69546d236249fff296a8c010ccf9", size = 32568781, upload-time = "2026-02-05T07:01:41.379Z" }, + { url = "https://files.pythonhosted.org/packages/3e/51/82fed528b45173bf629fa44effb76dff8bc9f4eeaee759038362dfa60237/opencv_python-4.13.0.92-cp37-abi3-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:0bc2596e68f972ca452d80f444bc404e08807d021fbba40df26b61b18e01838a", size = 47685527, upload-time = "2026-02-05T06:59:11.24Z" }, + { url = "https://files.pythonhosted.org/packages/db/07/90b34a8e2cf9c50fe8ed25cac9011cde0676b4d9d9c973751ac7616223a2/opencv_python-4.13.0.92-cp37-abi3-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:402033cddf9d294693094de5ef532339f14ce821da3ad7df7c9f6e8316da32cf", size = 70460872, upload-time = "2026-02-05T06:59:19.162Z" }, + { url = "https://files.pythonhosted.org/packages/02/6d/7a9cc719b3eaf4377b9c2e3edeb7ed3a81de41f96421510c0a169ca3cfd4/opencv_python-4.13.0.92-cp37-abi3-manylinux_2_28_aarch64.whl", hash = "sha256:bccaabf9eb7f897ca61880ce2869dcd9b25b72129c28478e7f2a5e8dee945616", size = 46708208, upload-time = "2026-02-05T06:59:15.419Z" }, + { url = "https://files.pythonhosted.org/packages/fd/55/b3b49a1b97aabcfbbd6c7326df9cb0b6fa0c0aefa8e89d500939e04aa229/opencv_python-4.13.0.92-cp37-abi3-manylinux_2_28_x86_64.whl", hash = "sha256:620d602b8f7d8b8dab5f4b99c6eb353e78d3fb8b0f53db1bd258bb1aa001c1d5", size = 72927042, upload-time = "2026-02-05T06:59:23.389Z" }, + { url = "https://files.pythonhosted.org/packages/fb/17/de5458312bcb07ddf434d7bfcb24bb52c59635ad58c6e7c751b48949b009/opencv_python-4.13.0.92-cp37-abi3-win32.whl", hash = "sha256:372fe164a3148ac1ca51e5f3ad0541a4a276452273f503441d718fab9c5e5f59", size = 30932638, upload-time = "2026-02-05T07:02:14.98Z" }, + { url = "https://files.pythonhosted.org/packages/e9/a5/1be1516390333ff9be3a9cb648c9f33df79d5096e5884b5df71a588af463/opencv_python-4.13.0.92-cp37-abi3-win_amd64.whl", hash = "sha256:423d934c9fafb91aad38edf26efb46da91ffbc05f3f59c4b0c72e699720706f5", size = 40212062, upload-time = "2026-02-05T07:02:12.724Z" }, +] + +[[package]] +name = "openpyxl" +version = "3.1.5" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "et-xmlfile" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/3d/f9/88d94a75de065ea32619465d2f77b29a0469500e99012523b91cc4141cd1/openpyxl-3.1.5.tar.gz", hash = "sha256:cf0e3cf56142039133628b5acffe8ef0c12bc902d2aadd3e0fe5878dc08d1050", size = 186464, upload-time = "2024-06-28T14:03:44.161Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/c0/da/977ded879c29cbd04de313843e76868e6e13408a94ed6b987245dc7c8506/openpyxl-3.1.5-py2.py3-none-any.whl", hash = "sha256:5282c12b107bffeef825f4617dc029afaf41d0ea60823bbb665ef3079dc79de2", size = 250910, upload-time = "2024-06-28T14:03:41.161Z" }, +] + [[package]] name = "packaging" version = "26.0" @@ -1083,6 +1123,8 @@ dependencies = [ { name = "networkx", version = "3.6.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, { name = "numpy", version = "2.4.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, + { name = "opencv-python" }, + { name = "openpyxl" }, { name = "pandas", version = "2.3.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, { name = "pandas", version = "3.0.0", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, { name = "pillow" }, @@ -1093,6 +1135,7 @@ dependencies = [ { name = "scikit-learn", version = "1.8.0", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, { name = "scipy", version = "1.15.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, { name = "scipy", version = "1.17.0", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, + { name = "seaborn" }, { name = "vtk" }, ] @@ -1102,12 +1145,15 @@ requires-dist = [ { name = "matplotlib" }, { name = "networkx" }, { name = "numpy" }, + { name = "opencv-python", specifier = ">=4.13.0.92" }, + { name = "openpyxl", specifier = ">=3.1.5" }, { name = "pandas" }, { name = "pillow" }, { name = "pyvista" }, { name = "scikit-image" }, { name = "scikit-learn" }, { name = "scipy" }, + { name = "seaborn" }, { name = "vtk" }, ] @@ -1513,6 +1559,22 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/a3/bb/bbae36d06c0fd670e8373da67096cd57058b57c9bad7d92969b5e3b730af/scooby-0.11.0-py3-none-any.whl", hash = "sha256:a79663d1a7711eb104e4b2935988ea1ed5f7be6b7288fad23b4fba7462832f9d", size = 19877, upload-time = "2025-11-01T19:22:53.046Z" }, ] +[[package]] +name = "seaborn" +version = "0.13.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "matplotlib" }, + { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "numpy", version = "2.4.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, + { name = "pandas", version = "2.3.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "pandas", version = "3.0.0", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/86/59/a451d7420a77ab0b98f7affa3a1d78a313d2f7281a57afb1a34bae8ab412/seaborn-0.13.2.tar.gz", hash = "sha256:93e60a40988f4d65e9f4885df477e2fdaff6b73a9ded434c1ab356dd57eefff7", size = 1457696, upload-time = "2024-01-25T13:21:52.551Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/83/11/00d3c3dfc25ad54e731d91449895a79e4bf2384dc3ac01809010ba88f6d5/seaborn-0.13.2-py3-none-any.whl", hash = "sha256:636f8336facf092165e27924f223d3c62ca560b1f2bb5dff7ab7fad265361987", size = 294914, upload-time = "2024-01-25T13:21:49.598Z" }, +] + [[package]] name = "six" version = "1.17.0"