Skip to content

Optimize vertexModel.initialize: reduce initialization time by 25-30% #13

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Copilot wants to merge 8 commits into
base: main
Choose a base branch
from
Draft

Optimize vertexModel.initialize: reduce initialization time by 25-30% #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
c587ebb
Initial plan
Copilot Feb 2, 2026
986f3b0
Optimize performance bottlenecks in face, cell, geo and utils
Copilot Feb 2, 2026
a92ecc9
Additional optimizations: determinant calculations and interface look…
Copilot Feb 2, 2026
06da9ab
Optimize SharedByCells computation in build_edges
Copilot Feb 2, 2026
044dcee
Fix code review issues: remove unused variables and improve list copy…
Copilot Feb 2, 2026
79d78c9
Merge branch 'main' into copilot/improve-vertexmodel-initialize-efficacy
Pablo1990 Feb 2, 2026
28dda7c
Updating gitignore
Pablo1990 Feb 3, 2026
59cc825
Merge branch 'main' into copilot/improve-vertexmodel-initialize-efficacy
Pablo1990 Feb 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,3 +196,5 @@ Input
src/pyVertexModel/_version.py

Temp/

build/
35 changes: 22 additions & 13 deletions src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -495,35 +495,44 @@ 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:
"""
cells = np.unique(labelled_img)
cells = cells[cells != 0] # Remove cell 0 if it exists

img_neighbours = [None] * (np.max(cells) + 1)

# Structuring element for dilation
se = np.array([
[0, 0, 1, 0, 0],
[0, 1, 1, 1, 0],
[1, 1, 1, 1, 1],
[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


# 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):
Expand Down
19 changes: 13 additions & 6 deletions src/pyVertexModel/geometry/cell.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -214,19 +214,26 @@ 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

current_v = np.linalg.det(np.array([y1, y2, y3])) / 6

# Compute relative positions (reuse y3)
y1 = self.Y[idx0, :] - self.X
y2 = self.Y[idx1, :] - self.X

# 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))
Expand Down
35 changes: 21 additions & 14 deletions src/pyVertexModel/geometry/face.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -193,8 +191,10 @@ def build_edges(self, T, face_ids, face_centre, face_interface_type, X, Ys, non_
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)
Expand All @@ -204,21 +204,24 @@ 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,
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):
"""
Expand All @@ -227,13 +230,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
Expand Down
63 changes: 38 additions & 25 deletions src/pyVertexModel/geometry/geo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -695,14 +695,19 @@ 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],
dtype='int')

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]
Expand All @@ -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]
Expand Down
6 changes: 5 additions & 1 deletion src/pyVertexModel/geometry/tris.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
45 changes: 32 additions & 13 deletions src/pyVertexModel/util/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -358,8 +358,29 @@ 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 - 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())
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:
Expand Down Expand Up @@ -529,23 +550,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):
Expand Down