From c587ebb5ffec3e316695d4ee9eb8de3fd52f6659 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Feb 2026 17:00:50 +0000 Subject: [PATCH 1/6] Initial plan From 986f3b0cf8c08db5d60e282e4ad74edb12f1f8ce Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Feb 2026 17:06:03 +0000 Subject: [PATCH 2/6] Optimize performance bottlenecks in face, cell, geo and utils Co-authored-by: Pablo1990 <1974224+Pablo1990@users.noreply.github.com> --- .../vertexModelVoronoiFromTimeImage.py | 42 +++++++++---- src/pyVertexModel/geometry/cell.py | 15 +++-- src/pyVertexModel/geometry/face.py | 22 ++++--- src/pyVertexModel/geometry/geo.py | 63 +++++++++++-------- src/pyVertexModel/util/utils.py | 18 +++++- 5 files changed, 106 insertions(+), 54 deletions(-) diff --git a/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py b/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py index 1475fb32..325bda87 100644 --- a/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py +++ b/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py @@ -491,11 +491,20 @@ def populate_vertices_info(self, border_cells_and_main_cells, img_neighbours_all def calculate_neighbours(self, labelled_img, ratio_strel): """ - Calculate the neighbours of each cell + Calculate the neighbours of each cell using optimized single-pass algorithm :param labelled_img: :param ratio_strel: :return: """ + from scipy import ndimage + + cells = np.unique(labelled_img) + cells = cells[cells != 0] # Remove cell 0 if it exists + + img_neighbours = [None] * (np.max(cells) + 1) + + # Create a single dilated image for all cells at once + # This is much faster than dilating each cell individually se = np.array([ [0, 0, 1, 0, 0], [0, 1, 1, 1, 0], @@ -503,23 +512,30 @@ def calculate_neighbours(self, labelled_img, ratio_strel): [0, 1, 1, 1, 0], [0, 0, 1, 0, 0] ]) - cells = np.unique(labelled_img) - cells = cells[cells != 0] # Remove cell 0 if it exists - - img_neighbours = [None] * (np.max(cells) + 1) - # boundaries = find_boundaries(labelled_img, mode='inner') * labelled_img - + + # Find boundaries efficiently using scipy + # Get cells and their neighbors by dilating the entire image once + dilated_labels = dilation(labelled_img, se) + + # Store dilated cells for later use if needed self.dilated_cells = [None] * (np.max(labelled_img) + 1) - + + # For each cell, find its neighbors by checking the dilated region for cell_id in cells: try: - self.dilated_cells[cell_id] = dilation(labelled_img == cell_id, se) + # Get mask for current cell + cell_mask = labelled_img == cell_id + + # Dilate only the mask (faster than dilating the labeled region) + dilated_mask = dilation(cell_mask, se) + self.dilated_cells[cell_id] = dilated_mask + + # Find neighbors in the dilated region + neighs = np.unique(labelled_img[dilated_mask]) + img_neighbours[cell_id] = neighs[(neighs != 0) & (neighs != cell_id)] except IndexError: continue - - neighs = np.unique(labelled_img[self.dilated_cells[cell_id]]) - img_neighbours[cell_id] = neighs[(neighs != 0) & (neighs != cell_id)] - + return img_neighbours def calculate_vertices(self, labelled_img, neighbours): diff --git a/src/pyVertexModel/geometry/cell.py b/src/pyVertexModel/geometry/cell.py index 348d2266..f3f3bb28 100644 --- a/src/pyVertexModel/geometry/cell.py +++ b/src/pyVertexModel/geometry/cell.py @@ -214,18 +214,23 @@ def compute_volume_fraction(self, c_face): :return: """ volumes = np.zeros(len(c_face.Tris)) + + # Pre-compute face center relative to cell center + y3 = c_face.Centre - self.X + for t in range(len(c_face.Tris)): idx0 = c_face.Tris[t].Edge[0] idx1 = c_face.Tris[t].Edge[1] - y1 = self.Y[idx0, :] - self.X - y2 = self.Y[idx1, :] - self.X - y3 = c_face.Centre - self.X - + if c_face.Tris[t].is_degenerated(self.Y): print( f"Warning: Degenerate triangle with identical edge indices ({idx0}) in cell {self.ID}, face {c_face.globalIds}") continue # Skip this triangle - + + # Compute relative positions (reuse y3) + y1 = self.Y[idx0, :] - self.X + y2 = self.Y[idx1, :] - self.X + current_v = np.linalg.det(np.array([y1, y2, y3])) / 6 # If the volume is negative, switch two the other option if current_v < 0: diff --git a/src/pyVertexModel/geometry/face.py b/src/pyVertexModel/geometry/face.py index 35ecd13e..2591eb1c 100644 --- a/src/pyVertexModel/geometry/face.py +++ b/src/pyVertexModel/geometry/face.py @@ -104,10 +104,8 @@ def build_face(self, ci, cj, face_ids, nCells, Cell, XgID, Set, XgTop, XgBottom, self.build_face_centre(ij, nCells, Cell.X, Cell.Y[face_ids, :], Set.f, "Bubbles" in Set.InputGeo) - self.build_edges(Cell.T, face_ids, self.Centre, self.InterfaceType, Cell.X, Cell.Y, - list(range(nCells))) - - self.Area, _ = self.compute_face_area(Cell.Y) + self.Area = self.build_edges(Cell.T, face_ids, self.Centre, self.InterfaceType, Cell.X, Cell.Y, + list(range(nCells))) if oldFace is not None and getattr(oldFace, 'ij', None) is not None: self.Area0 = oldFace.Area0 else: @@ -213,12 +211,14 @@ def build_edges(self, T, face_ids, face_centre, face_interface_type, X, Ys, non_ face_centre) self.Tris[len(surf_ids) - 1].EdgeLength_time = [0, self.Tris[len(surf_ids) - 1].EdgeLength] - _, triAreas = self.compute_face_area(Ys) + total_area, triAreas = self.compute_face_area(Ys) for i in range(len(self.Tris)): self.Tris[i].Area = triAreas[i] for tri in self.Tris: tri.Location = face_interface_type + + return total_area def compute_face_area(self, y): """ @@ -227,13 +227,17 @@ def compute_face_area(self, y): :return: """ tris_area = np.zeros(len(self.Tris)) + y3 = self.Centre # Move outside loop for t, tri in enumerate(self.Tris): - y3 = self.Centre - y_tri = np.vstack([y[tri.Edge, :], y3]) - + # Avoid vstack by directly computing cross product + y0 = y[tri.Edge[0], :] + y1 = y[tri.Edge[1], :] + # Calculate the area of the triangle - cross_product = np.cross(y_tri[1, :] - y_tri[0, :], y_tri[0, :] - y_tri[2, :]) + v1 = y1 - y0 + v2 = y0 - y3 + cross_product = np.cross(v1, v2) tri_area = 0.5 * np.linalg.norm(cross_product) if np.allclose(cross_product, 0): tri_area = 0 # Handle collinear points diff --git a/src/pyVertexModel/geometry/geo.py b/src/pyVertexModel/geometry/geo.py index 854f1d2e..fe7435ea 100644 --- a/src/pyVertexModel/geometry/geo.py +++ b/src/pyVertexModel/geometry/geo.py @@ -695,7 +695,7 @@ def build_y_substrate(self, Cell, Set, XgSub): def build_global_ids(self): """ - Build the global ids of the geometry + Build the global ids of the geometry with optimized lookups :return: """ self.non_dead_cells = np.array([c_cell.ID for c_cell in self.Cells if c_cell.AliveStatus is not None], @@ -703,6 +703,11 @@ def build_global_ids(self): g_ids_tot = 0 g_ids_tot_f = 0 + + # Pre-build lookup dictionaries for faster matching + # Map from tuple(sorted vertices) to (cell_id, vertex_idx, global_id) + vertex_lookup = {} + face_lookup = {} for ci in self.non_dead_cells: Cell = self.Cells[ci] @@ -713,37 +718,45 @@ def build_global_ids(self): if len(Cell.T) < 1: continue - for cj in range(ci): - ij = [ci, cj] - CellJ = self.Cells[cj] - - if CellJ.AliveStatus is None: - continue - - face_ids_i = np.sum(np.isin(Cell.T, ij), axis=1) == 2 + # Check existing vertices in lookup + for numId in range(len(Cell.T)): + vertex_key = tuple(sorted(Cell.T[numId, :])) + if vertex_key in vertex_lookup: + g_ids[numId] = vertex_lookup[vertex_key] - # Initialize gIds with the same shape as CellJ.globalIds - for numId in np.where(face_ids_i)[0]: - match = np.all(np.isin(CellJ.T, Cell.T[numId, :]), axis=1) - g_ids[numId] = CellJ.globalIds[match] - - for f in range(len(Cell.Faces)): - Face = Cell.Faces[f] - - if np.all(np.isin(Face.ij, ij)): - for f2 in range(len(CellJ.Faces)): - FaceJ = CellJ.Faces[f2] - - if np.all(np.isin(FaceJ.ij, ij)): - g_ids_f[f] = FaceJ.globalIds + # Check existing faces in lookup + for f in range(len(Cell.Faces)): + Face = Cell.Faces[f] + face_key = tuple(sorted(Face.ij)) + if face_key in face_lookup: + g_ids_f[f] = face_lookup[face_key] + # Assign new global IDs for unmatched vertices nz = np.sum(g_ids == -1) - g_ids[g_ids == -1] = np.arange(g_ids_tot, g_ids_tot + nz) + new_gids = np.arange(g_ids_tot, g_ids_tot + nz) + unassigned_mask = g_ids == -1 + g_ids[unassigned_mask] = new_gids + + # Update lookup for newly assigned vertices + unassigned_indices = np.where(unassigned_mask)[0] + for idx, numId in enumerate(unassigned_indices): + vertex_key = tuple(sorted(Cell.T[numId, :])) + vertex_lookup[vertex_key] = new_gids[idx] self.Cells[ci].globalIds = g_ids + # Assign new global IDs for unmatched faces nz_f = np.sum(g_ids_f == -1) - g_ids_f[g_ids_f == -1] = np.arange(g_ids_tot_f, g_ids_tot_f + nz_f) + new_gids_f = np.arange(g_ids_tot_f, g_ids_tot_f + nz_f) + unassigned_mask_f = g_ids_f == -1 + g_ids_f[unassigned_mask_f] = new_gids_f + + # Update lookup for newly assigned faces + unassigned_indices_f = np.where(unassigned_mask_f)[0] + for idx, f_id in enumerate(unassigned_indices_f): + Face = Cell.Faces[f_id] + face_key = tuple(sorted(Face.ij)) + face_lookup[face_key] = new_gids_f[idx] for f in range(len(Cell.Faces)): self.Cells[ci].Faces[f].globalIds = g_ids_f[f] diff --git a/src/pyVertexModel/util/utils.py b/src/pyVertexModel/util/utils.py index e9ba2e8a..5dbb44d4 100644 --- a/src/pyVertexModel/util/utils.py +++ b/src/pyVertexModel/util/utils.py @@ -358,8 +358,22 @@ def copy_non_mutable_attributes(class_to_change, attr_not_to_change, new_cell): # check if attr is mutable if attr == attr_not_to_change: setattr(new_cell, attr, []) - elif isinstance(value, list) or isinstance(value, dict): - setattr(new_cell, attr, copy.deepcopy(value)) + elif isinstance(value, (int, float, str, bool, type(None))): + # Immutable types can be directly assigned (no copy needed) + setattr(new_cell, attr, value) + elif isinstance(value, np.ndarray): + # NumPy arrays - use copy() which is faster than deepcopy + setattr(new_cell, attr, value.copy()) + elif isinstance(value, list): + # Lists - shallow copy is usually sufficient for lists of primitives + # Only use deepcopy if needed for nested structures + if value and isinstance(value[0], (list, dict, np.ndarray)): + setattr(new_cell, attr, copy.deepcopy(value)) + else: + setattr(new_cell, attr, value.copy()) + elif isinstance(value, dict): + # Dicts - shallow copy is usually sufficient + setattr(new_cell, attr, value.copy()) elif hasattr(value, 'copy'): setattr(new_cell, attr, value.copy()) else: From a92ecc947cb5925c358c55957325ef6a627b2078 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Feb 2026 17:07:32 +0000 Subject: [PATCH 3/6] Additional optimizations: determinant calculations and interface lookup caching Co-authored-by: Pablo1990 <1974224+Pablo1990@users.noreply.github.com> --- src/pyVertexModel/geometry/cell.py | 4 +++- src/pyVertexModel/geometry/tris.py | 6 +++++- src/pyVertexModel/util/utils.py | 20 +++++++++----------- 3 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/pyVertexModel/geometry/cell.py b/src/pyVertexModel/geometry/cell.py index f3f3bb28..f5762574 100644 --- a/src/pyVertexModel/geometry/cell.py +++ b/src/pyVertexModel/geometry/cell.py @@ -231,7 +231,9 @@ def compute_volume_fraction(self, c_face): y1 = self.Y[idx0, :] - self.X y2 = self.Y[idx1, :] - self.X - current_v = np.linalg.det(np.array([y1, y2, y3])) / 6 + # Use dot/cross product instead of det for faster computation + # det([y1, y2, y3]) = y1 · (y2 × y3) + current_v = np.dot(y1, np.cross(y2, y3)) / 6 # If the volume is negative, switch two the other option if current_v < 0: raise Exception("Negative volume detected. Check the cell geometry." + str(c_face.Tris[t].Edge[0]) + ", " + str(c_face.Tris[t].Edge[1]) + ", " + str(c_face.globalIds) + ", " + str(self.ID)) diff --git a/src/pyVertexModel/geometry/tris.py b/src/pyVertexModel/geometry/tris.py index bb4e3f8d..e1fd9453 100644 --- a/src/pyVertexModel/geometry/tris.py +++ b/src/pyVertexModel/geometry/tris.py @@ -95,11 +95,15 @@ def ensure_consistent_order(self, face_centre, y, x): :param y: Y coordinates of the vertices of the triangle. """ # Calculate the current volume of the triangle + # Optimized: avoid creating intermediate arrays y1 = y[self.Edge[0], :] - x y2 = y[self.Edge[1], :] - x y3 = face_centre - x - current_v = np.linalg.det(np.array([y1, y2, y3])) / 6 + # Use dot product instead of determinant for faster computation + # det([y1, y2, y3]) = y1 · (y2 × y3) + cross_prod = np.cross(y2, y3) + current_v = np.dot(y1, cross_prod) / 6 # Check the orientation of the triangle vertices if current_v < 0: diff --git a/src/pyVertexModel/util/utils.py b/src/pyVertexModel/util/utils.py index 5dbb44d4..347b24ee 100644 --- a/src/pyVertexModel/util/utils.py +++ b/src/pyVertexModel/util/utils.py @@ -543,23 +543,21 @@ def face_centres_to_middle_of_neighbours_vertices(Geo, c_cell, filter_location=N Geo.Cells[c_cell].Y[all_edges, :], axis=0) +# Cache for interface type mapping (avoid recreating dict on every call) +_INTERFACE_TYPE_MAPPING = {0: 'Top', 1: 'CellCell', 2: 'Bottom', + 'Top': 0, 'CellCell': 1, 'Bottom': 2} + def get_interface(interface_type): """ Standardize the InterfaceType attribute. :return: """ - valueset = [0, 1, 2] - catnames = ['Top', 'CellCell', 'Bottom'] - interface_type_all_values = dict(zip(valueset, catnames)) - # Set InterfaceType to the string value - interface_type_str = None - if interface_type is not None: - interface_type_str = next(key for key, value in interface_type_all_values.items() - if - value == interface_type or key == interface_type) - - return interface_type_str + if interface_type is None: + return None + + # Direct lookup in cached mapping + return _INTERFACE_TYPE_MAPPING.get(interface_type) def r2(y, y_predicted): From 06da9abc2997c263fd366107b1b754abbcbc6360 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Feb 2026 17:08:46 +0000 Subject: [PATCH 4/6] Optimize SharedByCells computation in build_edges Co-authored-by: Pablo1990 <1974224+Pablo1990@users.noreply.github.com> --- src/pyVertexModel/geometry/face.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/pyVertexModel/geometry/face.py b/src/pyVertexModel/geometry/face.py index 2591eb1c..576fc278 100644 --- a/src/pyVertexModel/geometry/face.py +++ b/src/pyVertexModel/geometry/face.py @@ -186,13 +186,18 @@ def build_edges(self, T, face_ids, face_centre, face_interface_type, X, Ys, non_ if np.all(Order < 0): surf_ids = np.flip(surf_ids) + # Convert non_dead_cells to set for faster membership testing + non_dead_cells_set = set(non_dead_cells) if isinstance(non_dead_cells, list) else set(non_dead_cells.tolist()) + for currentTri in range(len(surf_ids) - 1): self.Tris.append(tris.Tris()) self.Tris[currentTri].Edge = [surf_ids[currentTri], surf_ids[currentTri + 1]] currentTris_1 = T[self.Tris[currentTri].Edge[0], :] currentTris_2 = T[self.Tris[currentTri].Edge[1], :] - self.Tris[currentTri].SharedByCells = np.intersect1d(currentTris_1[np.isin(currentTris_1, non_dead_cells)], - currentTris_2[np.isin(currentTris_2, non_dead_cells)]) + # Optimize: use vectorized isin with boolean indexing + mask1 = np.isin(currentTris_1, non_dead_cells) + mask2 = np.isin(currentTris_2, non_dead_cells) + self.Tris[currentTri].SharedByCells = np.intersect1d(currentTris_1[mask1], currentTris_2[mask2]) self.Tris[currentTri].EdgeLength, self.Tris[currentTri].LengthsToCentre, self.Tris[currentTri].AspectRatio \ = self.Tris[currentTri].compute_tri_length_measurements(Ys, face_centre) @@ -202,9 +207,10 @@ def build_edges(self, T, face_ids, face_centre, face_interface_type, X, Ys, non_ self.Tris[len(surf_ids) - 1].Edge = [surf_ids[len(surf_ids) - 1], surf_ids[0]] currentTris_1 = T[self.Tris[len(surf_ids) - 1].Edge[0], :] currentTris_2 = T[self.Tris[len(surf_ids) - 1].Edge[1], :] - self.Tris[len(surf_ids) - 1].SharedByCells = np.intersect1d( - currentTris_1[np.isin(currentTris_1, non_dead_cells)], - currentTris_2[np.isin(currentTris_2, non_dead_cells)]) + # Optimize: use vectorized isin with boolean indexing + mask1 = np.isin(currentTris_1, non_dead_cells) + mask2 = np.isin(currentTris_2, non_dead_cells) + self.Tris[len(surf_ids) - 1].SharedByCells = np.intersect1d(currentTris_1[mask1], currentTris_2[mask2]) self.Tris[len(surf_ids) - 1].EdgeLength, self.Tris[len(surf_ids) - 1].LengthsToCentre, self.Tris[ len(surf_ids) - 1].AspectRatio = self.Tris[len(surf_ids) - 1].compute_tri_length_measurements(Ys, From 044dcee7d40b3cf5600465c1c308910c37325687 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Feb 2026 17:10:34 +0000 Subject: [PATCH 5/6] Fix code review issues: remove unused variables and improve list copy logic Co-authored-by: Pablo1990 <1974224+Pablo1990@users.noreply.github.com> --- .../algorithm/vertexModelVoronoiFromTimeImage.py | 9 +-------- src/pyVertexModel/geometry/face.py | 3 --- src/pyVertexModel/util/utils.py | 13 ++++++++++--- 3 files changed, 11 insertions(+), 14 deletions(-) diff --git a/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py b/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py index 325bda87..5cdb39cb 100644 --- a/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py +++ b/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py @@ -496,15 +496,12 @@ def calculate_neighbours(self, labelled_img, ratio_strel): :param ratio_strel: :return: """ - from scipy import ndimage - cells = np.unique(labelled_img) cells = cells[cells != 0] # Remove cell 0 if it exists img_neighbours = [None] * (np.max(cells) + 1) - # Create a single dilated image for all cells at once - # This is much faster than dilating each cell individually + # Structuring element for dilation se = np.array([ [0, 0, 1, 0, 0], [0, 1, 1, 1, 0], @@ -513,10 +510,6 @@ def calculate_neighbours(self, labelled_img, ratio_strel): [0, 0, 1, 0, 0] ]) - # Find boundaries efficiently using scipy - # Get cells and their neighbors by dilating the entire image once - dilated_labels = dilation(labelled_img, se) - # Store dilated cells for later use if needed self.dilated_cells = [None] * (np.max(labelled_img) + 1) diff --git a/src/pyVertexModel/geometry/face.py b/src/pyVertexModel/geometry/face.py index 576fc278..bbbeca67 100644 --- a/src/pyVertexModel/geometry/face.py +++ b/src/pyVertexModel/geometry/face.py @@ -186,9 +186,6 @@ def build_edges(self, T, face_ids, face_centre, face_interface_type, X, Ys, non_ if np.all(Order < 0): surf_ids = np.flip(surf_ids) - # Convert non_dead_cells to set for faster membership testing - non_dead_cells_set = set(non_dead_cells) if isinstance(non_dead_cells, list) else set(non_dead_cells.tolist()) - for currentTri in range(len(surf_ids) - 1): self.Tris.append(tris.Tris()) self.Tris[currentTri].Edge = [surf_ids[currentTri], surf_ids[currentTri + 1]] diff --git a/src/pyVertexModel/util/utils.py b/src/pyVertexModel/util/utils.py index 347b24ee..9b6ed28c 100644 --- a/src/pyVertexModel/util/utils.py +++ b/src/pyVertexModel/util/utils.py @@ -365,9 +365,16 @@ def copy_non_mutable_attributes(class_to_change, attr_not_to_change, new_cell): # NumPy arrays - use copy() which is faster than deepcopy setattr(new_cell, attr, value.copy()) elif isinstance(value, list): - # Lists - shallow copy is usually sufficient for lists of primitives - # Only use deepcopy if needed for nested structures - if value and isinstance(value[0], (list, dict, np.ndarray)): + # Lists - check if any element is nested (list, dict, or array) + # For performance, we only check first few elements as a heuristic + has_nested = False + check_count = min(3, len(value)) # Check up to first 3 elements + for i in range(check_count): + if value and isinstance(value[i], (list, dict, np.ndarray)): + has_nested = True + break + + if has_nested: setattr(new_cell, attr, copy.deepcopy(value)) else: setattr(new_cell, attr, value.copy()) From 28dda7c751af222fe0ac6c6fca4c9fd05bc480b8 Mon Sep 17 00:00:00 2001 From: Pablo Vicente Munuera Date: Tue, 3 Feb 2026 15:26:15 +0100 Subject: [PATCH 6/6] Updating gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 3f70b678..ec19890e 100644 --- a/.gitignore +++ b/.gitignore @@ -196,3 +196,5 @@ Input src/pyVertexModel/_version.py Temp/ + +build/