Handle degenerate 4-fold Voronoi vertices causing unpacking errors in _assemble_forces

[wwang721]

Root Cause

Degenerate 4-fold Voronoi vertices cause a crash in FiniteVoronoiSimulator._assemble_forces when more than three seed-point indices are found for a Voronoi vertex, observed in large systems (N=1000+). The code (v0.4.7) assumes each inner vertex connects to three cells and unpacks accordingly, leading to ValueError when this is not true:

Traceback (most recent call last):
  File "xxx.py", line xxx, in <module>
    diag = sim.build()
           ^^^^^^^^^^^
  File "/xxx/lib/python3.11/site-packages/pyafv/finite_voronoi.py", line 758, in build
    F = self._assemble_forces(
        ^^^^^^^^^^^^^^^^^^^^^^
  File "/xxx/lib/python3.11/site-packages/pyafv/finite_voronoi.py", line 528, in _assemble_forces
    I[t], J[t], K[t] = vertex_points[h]
    ^^^^^^^^^^^^^^^^
ValueError: too many values to unpack (expected 3)

Detailed Explanation

• Code in question:

  I[t], J[t], K[t] = vertex_points[h]   # blows up if len > 3

• When does this happen?
  • In rare cases, Qhull returns a single vertex with 4 incident ridges — typically when 4 (or more) points are co-circular (e.g., corners of a square or crystalline patches at high N).
  • This increases in likelihood at larger N due to combinatorics and relaxed crystalline formation.

Minimal Example to Reproduce

import numpy as np
import pyafv as afv

square_pts = np.array([[ 1., 1.], [-1., 1.], [-1., -1.], [ 1., -1.]]) * 0.5
params = afv.PhysicalParams(r=1.0)
sim = afv.FiniteVoronoiSimulator(square_pts, params)
diag = sim.build()   # raises: ValueError: too many values to unpack (expected 3)

[wwang721]

This issue arose during large-batch simulations with N=1000 using random initial positions. This is strange, as there should be no co-circular degeneracies (four points lying exactly on a common circumcircle) given the small floating-point tolerances of SciPy QHull. For example, consider the four co-circular points above. If we apply a small perturbation to the first point,

import numpy as np
import pyafv as afv

eps = 1e-16  # small perturbation
square_pts = np.array([[ 1., 1. + eps], [-1., 1.], [-1., -1.], [ 1., -1.]]) * 0.5
params = afv.PhysicalParams(r=1.0)
sim = afv.FiniteVoronoiSimulator(square_pts, params)
diag = sim.build()   # raises: ValueError: too many values to unpack (expected 3)

the ValueError reappears only when eps is on the order of $10^{-16}$, whereas $10^{-15}$ works without issue.

[wwang721]

The source of the issue might be attributed to the fact that Qhull’s floating-point tolerance scales with the overall coordinate span of the input points. To illustrate this, consider adding 100 random points around the four nearly co-circular points:

import numpy as np
import pyafv as afv

eps = 1e-10
square_pts = np.array([[1., 1. + eps], [-1., 1.], [-1., -1.], [1., -1.]]) * 0.5

np.random.seed(42)                        # For reproducibility
N = 100                                   # number of additional points
span = 100.0                              # span of the random points
pts_add = np.random.rand(N, 2) * span
pts_total = np.vstack((square_pts, pts_add))

params = afv.PhysicalParams(r=1.0)
sim = afv.FiniteVoronoiSimulator(pts_total, params)
diag = sim.build()   # raises: ValueError: too many values to unpack (expected 3)

with an added point cloud spanning ~100, the critical eps shifts to $10^{-10}$. Increasing the span further to 1000 raises this threshold to around $10^{-8}$.

[wwang721]

See scipy/scipy#10464