diff --git a/CHANGELOG.md b/CHANGELOG.md index 7ee76c8..34f360c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,25 @@ All notable changes to this project are documented in this file. The format is based on *Keep a Changelog*, and this project follows *Semantic Versioning*. +## [0.4.2] - 2026-03-04 + +### Changed + +- Vendored Voro++: updated the vendored snapshot to include the upstream numeric robustness fix for + power/Laguerre (radical) pruning (fixes rare cross-platform edge cases in fully periodic power + tessellations). +- Removed the previously vendored local `nextafter`-based `max_radius` inflation patch (no longer needed). + +## [0.4.1] - 2026-02-16 + +### Fixed + +- Vendored Voro++: inflate the stored global `max_radius` by 1 ULP (via `nextafter`) in + power/Laguerre mode to make radical pruning robust across platforms. +- Removed a Python-side workaround that recomputed fully periodic orthorhombic power tessellations + via the periodic (triclinic) backend when periodic face-shift assignment failed. +- Updated documentation to reflect the patched vendored Voro++ snapshot. + ## [0.4.0] - 2026-02-15 Initial public release. diff --git a/NOTICE.md b/NOTICE.md index b7016ec..dc7f14d 100644 --- a/NOTICE.md +++ b/NOTICE.md @@ -9,3 +9,7 @@ This project vendors and links against the following third-party software: - License: See `vendor/voro++/LICENSE` pyvoro2 is licensed under the MIT License (see `LICENSE`). The included Voro++ code remains under its original license. + +Local modifications: + +- None. The vendored Voro++ snapshot is kept unmodified (aside from being vendored into this repository). diff --git a/README.md b/README.md index 3b5e454..f3a2ac7 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,6 @@ # pyvoro2 [![CI](https://github.com/IvanChernyshov/pyvoro2/actions/workflows/ci.yml/badge.svg)](https://github.com/IvanChernyshov/pyvoro2/actions/workflows/ci.yml) [![Docs](https://github.com/IvanChernyshov/pyvoro2/actions/workflows/docs.yml/badge.svg)](https://github.com/IvanChernyshov/pyvoro2/actions/workflows/docs.yml) [![PyPI](https://img.shields.io/pypi/v/pyvoro2.svg)](https://pypi.org/project/pyvoro2/) [![Python Versions](https://img.shields.io/pypi/pyversions/pyvoro2.svg)](https://pypi.org/project/pyvoro2/) [![License](https://img.shields.io/pypi/l/pyvoro2.svg)](https://github.com/IvanChernyshov/pyvoro2/blob/main/LICENSE) - **Documentation:** https://IvanChernyshov.github.io/pyvoro2/ @@ -19,7 +18,7 @@ matter physics — especially **periodic boundary conditions** and **neighbor gr pyvoro2 is designed to be **honest and predictable**: -- it vendors and wraps **unmodified** upstream Voro++; +- it vendors and wraps an upstream Voro++ snapshot (with a small numeric robustness patch for power/Laguerre diagrams); - the core tessellation modes are **standard Voronoi** and **power/Laguerre**. ## Quickstart @@ -100,19 +99,10 @@ For stricter post-hoc checks, see: - `pyvoro2.validate_tessellation(..., level='strict')` - `pyvoro2.validate_normalized_topology(..., level='strict')` -## Platform note (macOS) - -On some macOS builds, fully periodic **power/Laguerre** tessellations can -occasionally produce a non-reciprocal face/neighbor graph. In that case, -requesting `return_face_shifts=True` may raise a `ValueError` because pyvoro2 -cannot assign consistent periodic shifts. - -Workarounds: - -- Avoid requesting shifts (`return_face_shifts=False`), or -- disable shift validation (`validate_face_shifts=False`, shifts may be - unreliable), or -- run strict periodic power workflows on Linux/Windows. +Note: pyvoro2 vendors a Voro++ snapshot that includes the upstream numeric robustness fix for +*power/Laguerre* mode (radical pruning). This avoids rare cross-platform edge cases where fully +periodic power tessellations could yield a non-reciprocal face/neighbor graph under aggressive +floating-point codegen. ## Why use pyvoro2 diff --git a/docs/guide/domains.md b/docs/guide/domains.md index 4458816..8fb5cab 100644 --- a/docs/guide/domains.md +++ b/docs/guide/domains.md @@ -60,18 +60,13 @@ Notes: - `return_face_shifts=True` is supported whenever at least one axis is periodic. - For non-periodic axes, query points outside the bounds are treated as out-of-domain. -Platform note (macOS): +Robustness note: -- On some macOS builds (compiler/CPU-dependent), fully periodic **power/Laguerre** - tessellations can occasionally produce a non-reciprocal face/neighbor graph. - In that case, requesting `return_face_shifts=True` may raise a `ValueError` - because pyvoro2 cannot assign consistent periodic shifts. - - Workarounds (choose what matches your workflow): - - - Disable shift validation: `validate_face_shifts=False` (shifts may be unreliable). - - Avoid requesting shifts: `return_face_shifts=False`. - - Use a non-macOS build for strict periodic power workflows. +- pyvoro2 vendors a Voro++ snapshot that includes the upstream numeric robustness fix for fully + periodic **power/Laguerre** tessellations (radical pruning). +- If `return_face_shifts=True` fails with a `ValueError`, it is typically due to a + genuine geometric degeneracy (for example, nearly co-spherical sites) rather than + a platform-specific issue. ## `PeriodicCell` (triclinic, fully periodic) diff --git a/docs/index.md b/docs/index.md index a4e516e..831aa49 100644 --- a/docs/index.md +++ b/docs/index.md @@ -12,7 +12,7 @@ matter physics — especially **periodic boundary conditions** and **neighbor gr pyvoro2 is designed to be **honest and predictable**: -- it vendors and wraps **unmodified** upstream Voro++; +- it vendors and wraps an upstream Voro++ snapshot (with a small numeric robustness patch for power/Laguerre diagrams); - the core tessellation modes are **standard Voronoi** and **power/Laguerre**. ## Quickstart @@ -93,6 +93,11 @@ For stricter post-hoc checks, see: - `pyvoro2.validate_tessellation(..., level='strict')` - `pyvoro2.validate_normalized_topology(..., level='strict')` +Note: pyvoro2 vendors a Voro++ snapshot that includes the upstream numeric robustness fix for +*power/Laguerre* mode (radical pruning). This avoids rare cross-platform edge cases where fully +periodic power tessellations could yield a non-reciprocal face/neighbor graph under aggressive +floating-point codegen. + ## Why use pyvoro2 Voro++ is fast and feature-rich, but it is a C++ library with a low-level API. diff --git a/docs/project/about.md b/docs/project/about.md index 633e41e..5cea2d0 100644 --- a/docs/project/about.md +++ b/docs/project/about.md @@ -20,8 +20,9 @@ At the core, pyvoro2 exposes only two mathematically standard tessellations: Voro++ is a widely used C++ library for computing Voronoi cells efficiently in 3D. It implements robust algorithms and is commonly used in computational physics and materials science. -pyvoro2 vendors an **unmodified** snapshot of upstream Voro++ and builds its Python extension against it. -This keeps the core geometry trustworthy and makes it easy to track upstream behavior. +pyvoro2 vendors a snapshot of upstream Voro++ and builds its Python extension against it. +The vendored snapshot includes the upstream numeric robustness fix for *power/Laguerre* (radical) pruning, +which avoids rare cross-platform edge cases in fully periodic power tessellations. ## When should you use pyvoro2? diff --git a/src/pyvoro2/__about__.py b/src/pyvoro2/__about__.py index db97f50..eb6a4a6 100644 --- a/src/pyvoro2/__about__.py +++ b/src/pyvoro2/__about__.py @@ -5,4 +5,4 @@ """ # Keep this as a simple assignment so scikit-build-core can extract it via regex. -__version__ = '0.4.0' +__version__ = '0.4.2' diff --git a/src/pyvoro2/api.py b/src/pyvoro2/api.py index 0fa3b32..4714983 100644 --- a/src/pyvoro2/api.py +++ b/src/pyvoro2/api.py @@ -342,85 +342,6 @@ def compute( else tuple(bool(x) for x in domain.periodic) ) is_periodic = isinstance(domain, OrthorhombicCell) and any(periodic_flags) - full_periodic_ortho = ( - isinstance(domain, OrthorhombicCell) - and periodic_flags == (True, True, True) - ) - - def _compute_full_periodic_ortho_via_periodic_backend( - *, - mode: Literal['standard', 'power'], - rr: np.ndarray | None, - ) -> list[dict[str, Any]]: - """Compute fully periodic OrthorhombicCell via Voro++ periodic containers. - - This is a fallback path for rare platform-dependent issues in the - rectangular container backend (or its face bookkeeping). - - It remaps points into the primary cell, shifts the origin to (0,0,0), - calls the periodic container backend, then shifts outputs back. - """ - assert isinstance(domain, OrthorhombicCell) - (xmin, xmax), (ymin, ymax), (zmin, zmax) = bounds - origin = np.asarray((xmin, ymin, zmin), dtype=np.float64) - bx = float(xmax - xmin) - by = float(ymax - ymin) - bz = float(zmax - zmin) - - # Wrap into the primary domain (half-open in periodic axes) - pts_w = domain.remap_cart(pts, return_shifts=False) - pts_i = pts_w - origin - cell_params = (bx, 0.0, by, 0.0, 0.0, bz) - - if mode == 'standard': - cells_i = core.compute_periodic_standard( - pts_i, - ids_internal, - cell_params, - (nx, ny, nz), - init_mem, - opts, - ) - elif mode == 'power': - if rr is None: - raise ValueError('internal error: rr is required for power mode') - cells_i = core.compute_periodic_power( - pts_i, - ids_internal, - rr, - cell_params, - (nx, ny, nz), - init_mem, - opts, - ) - else: - raise ValueError(f'unknown mode: {mode}') - - if include_empty: - # Voro++ remaps inserted points into the primary cell; mirror that - # for any empty-cell records we inject. - _add_empty_cells_inplace(cells_i, n=n, sites=pts_i, opts=opts) - - # Shift outputs back to the original bounds origin. - if float(xmin) != 0.0 or float(ymin) != 0.0 or float(zmin) != 0.0: - ox, oy, oz = float(origin[0]), float(origin[1]), float(origin[2]) - for cc in cells_i: - s = cc.get('site') - if s is not None and len(s) == 3: - cc['site'] = [ - float(s[0]) + ox, - float(s[1]) + oy, - float(s[2]) + oz, - ] - if return_vertices: - vv = cc.get('vertices') - if vv: - cc['vertices'] = [ - [float(x) + ox, float(y) + oy, float(z) + oz] - for (x, y, z) in vv - ] - return cells_i - if return_face_shifts: if not is_periodic: raise ValueError( @@ -479,50 +400,21 @@ def _compute_full_periodic_ortho_via_periodic_backend( if return_face_shifts: assert isinstance(domain, OrthorhombicCell) a, b, cvec = domain.lattice_vectors - try: - _add_periodic_face_shifts_inplace( - cells, - lattice_vectors=(a, b, cvec), - periodic_mask=periodic_flags, - mode=mode, - radii=( - np.asarray(radii, dtype=np.float64) - if radii is not None - else None - ), - search=int(face_shift_search), - tol=face_shift_tol, - validate=bool(validate_face_shifts), - repair=bool(repair_face_shifts), - ) - except ValueError: - # Fallback for rare platform-dependent issues (observed on macOS) - # where the rectangular periodic backend can produce a face/neighbor - # graph that prevents consistent shift assignment. - if full_periodic_ortho and mode == 'power': - cells = _compute_full_periodic_ortho_via_periodic_backend( - mode=mode, - rr=rr, - ) - # Re-run shift assignment on fallback output. - _add_periodic_face_shifts_inplace( - cells, - lattice_vectors=(a, b, cvec), - periodic_mask=(True, True, True), - mode=mode, - radii=( - np.asarray(radii, dtype=np.float64) - if radii is not None - else None - ), - search=int(face_shift_search), - tol=face_shift_tol, - validate=bool(validate_face_shifts), - repair=bool(repair_face_shifts), - ) - else: - raise - + _add_periodic_face_shifts_inplace( + cells, + lattice_vectors=(a, b, cvec), + periodic_mask=periodic_flags, + mode=mode, + radii=( + np.asarray(radii, dtype=np.float64) + if radii is not None + else None + ), + search=int(face_shift_search), + tol=face_shift_tol, + validate=bool(validate_face_shifts), + repair=bool(repair_face_shifts), + ) if ids_user is not None: _remap_ids_inplace(cells, ids_user) diff --git a/tests/test_power_periodic_pruning_regression.py b/tests/test_power_periodic_pruning_regression.py new file mode 100644 index 0000000..1fcd6e7 --- /dev/null +++ b/tests/test_power_periodic_pruning_regression.py @@ -0,0 +1,119 @@ +import numpy as np + +from pyvoro2 import OrthorhombicCell, compute + + +def test_power_periodic_pruning_regression_issue43(): + """Regression for Voro++ power/Laguerre pruning edge case (upstream issue #43). + + The invariant checks are intentionally simple: + + - Sum of cell volumes should equal the periodic box volume. + - Neighbor reciprocity between particles 0 and 9 should hold. + + This dataset was found to produce cross-platform inconsistencies under + sufficiently aggressive floating-point codegen in older Voro++ snapshots. + """ + + # Fully periodic orthorhombic domain [0, L) in all axes. + L = 8.0 + domain = OrthorhombicCell( + ((0.0, L), (0.0, L), (0.0, L)), periodic=(True, True, True) + ) + + # Input points and radii (power diagram / radical Voronoi). + # + # These values are copied from the upstream reproduction in issue #43. + # They were generated via: + # rng = np.random.default_rng(3) + # pts = rng.uniform(0.0, 8.0, size=(25, 3)) + # radii = rng.uniform(0.1, 0.5, size=(25,)) + pts = np.array( + [ + [0.68519333714899489, 1.8944840527687976, 6.4101957216511751], + [4.6572962885149423, 0.75302913792319348, 3.4650155218917904], + [3.8324103851266722, 1.2779113170966285, 5.8766172112737163], + [0.90937615937122729, 3.1298255239652963, 4.1339214609709094], + [3.4450241633134224, 4.6943885715051259, 5.9027022983372817], + [7.6501380386887883, 2.2736093099903316, 5.1883776566386004], + [5.5697279733612435, 2.341765992099897, 0.011920668070689366], + [7.7876821981313018, 2.3872097841350053, 2.5118880162746944], + [7.1336885635612575, 4.6813035191272645, 3.770477321454651], + [6.1862160771905312, 0.24276806129976958, 5.6557207652449879], + [2.9939506678277663, 0.72682170803406265, 5.2840005394231581], + [7.451710837930836, 1.65752934464801, 5.0407215982827442], + [2.385304725259398, 5.9340534405546434, 5.7773184651369397], + [1.7497233965504364, 6.6390949941944982, 5.2612176869859457], + [5.4623912628828019, 6.5606060013642797, 3.4285832343876956], + [6.069643689239352, 7.0278414773300308, 0.81855937537659516], + [6.79814669972923, 3.151418661058681, 3.8374713880981997], + [1.1706765580655878, 5.5874107595767479, 2.335828927902808], + [6.9691131983487129, 2.2029950155846167, 4.4944777498471194], + [3.197249769043622, 4.9032759352195133, 1.5731139181769898], + [1.4423003267447614, 5.9748830851987027, 6.0177873469542187], + [4.5358229956232979, 7.3686374182634395, 1.6462003858321506], + [6.8072089967332117, 1.3518984906166898, 7.7148617671542366], + [4.989541824849214, 4.8550703040801178, 7.7644701050609902], + [6.2962617103101679, 6.3193398087084347, 0.4327500629677221], + ], + dtype=float, + ) + + radii = np.array( + [ + 0.24771452243776798, + 0.13395790886618722, + 0.17741103258620577, + 0.18554679627685888, + 0.44345677286636176, + 0.15070199188732217, + 0.2187031073557798, + 0.29713879156973305, + 0.43978416447581475, + 0.4860925683349502, + 0.3832579082748967, + 0.18547488777221927, + 0.31799307258320164, + 0.3823836171947548, + 0.12075301109065006, + 0.3719536668896286, + 0.24731273405563245, + 0.33588011471168333, + 0.36781331245666604, + 0.36765099845837435, + 0.30921991021034256, + 0.3218949729860357, + 0.1792598854454132, + 0.2980749566141609, + 0.15016378905810127, + ], + dtype=float, + ) + + cells = compute( + pts, + domain=domain, + mode='power', + radii=radii, + # Explicitly match the upstream repro harness (and the wrapper default). + blocks=(1, 1, 1), + init_mem=8, + return_vertices=False, + return_adjacency=False, + return_faces=True, + ) + + # Invariant 1: volume sum equals box volume. + vol_sum = float(sum(c["volume"] for c in cells)) + assert np.isfinite(vol_sum) + assert abs(vol_sum - (L * L * L)) < 1e-6 + + # Invariant 2: reciprocity for the known problematic pair. + by_id = {int(c["id"]): c for c in cells} + assert 0 in by_id and 9 in by_id + + neigh0 = {int(f["adjacent_cell"]) for f in by_id[0]["faces"]} + neigh9 = {int(f["adjacent_cell"]) for f in by_id[9]["faces"]} + + assert 9 in neigh0 + assert 0 in neigh9 diff --git a/tests/test_strict_validation.py b/tests/test_strict_validation.py index dc010b1..f24fdee 100644 --- a/tests/test_strict_validation.py +++ b/tests/test_strict_validation.py @@ -2,7 +2,6 @@ import numpy as np import pytest -import sys import pyvoro2 @@ -79,29 +78,16 @@ def test_validate_tessellation_strict_passes_for_power_mode_periodic_case() -> N pts = rng.uniform(0.0, 8.0, size=(25, 3)) radii = rng.uniform(0.1, 0.5, size=(25,)) - try: - cells = pyvoro2.compute( - pts, - domain=dom, - mode='power', - radii=radii, - return_vertices=True, - return_faces=True, - return_face_shifts=True, - repair_face_shifts=True, - ) - except ValueError as e: - # Known upstream/platform-specific instability: on some macOS builds, - # Voro++ can return a non-reciprocal face/neighbor graph for fully - # periodic *power* diagrams, which prevents pyvoro2 from assigning - # consistent periodic face shifts. - # - # We accept this as an expected failure on macOS for now and document - # workarounds in the README/docs. - if sys.platform == 'darwin': - pytest.xfail(f'macOS power-periodic instability: {e}') - raise - + cells = pyvoro2.compute( + pts, + domain=dom, + mode='power', + radii=radii, + return_vertices=True, + return_faces=True, + return_face_shifts=True, + repair_face_shifts=True, + ) diag = pyvoro2.validate_tessellation(cells, dom, level='strict') assert diag.ok_volume diff --git a/vendor/voro++/2d/src/rad_option.hh b/vendor/voro++/2d/src/rad_option.hh index c53a9f8..935f34e 100644 --- a/vendor/voro++/2d/src/rad_option.hh +++ b/vendor/voro++/2d/src/rad_option.hh @@ -94,8 +94,9 @@ class radius_poly { * \param[in] ijk the block that the particle is within. * \param[in] s the index of the particle within the block. */ inline void r_init(int ijk,int s) { - r_rad=ppr[ijk][3*s+2]*ppr[ijk][3*s+2]; - r_mul=r_rad-max_radius*max_radius; + double &r=ppr[ijk][3*s+2]; + r_rad=r*r; + r_mul=(r-max_radius)*(r+max_radius); } /** Sets a required constant to be used when carrying out a * plane bounds check. */ @@ -106,7 +107,10 @@ class radius_poly { * vertex multiplied by two. * \return True if particles at this radius could not possibly * cut the cell, false otherwise. */ - inline bool r_ctest(double crs,double mrs) {return crs+r_mul>sqrt(mrs*crs);} + inline bool r_ctest(double crs,double mrs) { + double cc=crs+r_mul; + return cc>0&&cc*cc>mrs*crs; + } /** Scales a plane displacement during a plane bounds check. * \param[in] lrs the plane displacement. * \return The scaled value. */ @@ -146,8 +150,8 @@ class radius_poly { inline bool r_scale_check(double &rs,double mrs,int ijk,int q) { double trs=rs; rs+=r_rad-ppr[ijk][3*q+2]*ppr[ijk][3*q+2]; - return rssqrt(mrs*crs);} + inline bool r_ctest(double crs,double mrs) { + double cc=crs+r_mul; + return cc>0&&cc*cc>mrs*crs; + } /** Scales a plane displacement during a plane bounds check. * \param[in] lrs the plane displacement. * \return The scaled value. */ @@ -148,7 +152,7 @@ class radius_poly { inline bool r_scale_check(double &rs,double mrs,int ijk,int q) { double trs=rs; rs+=r_rad-ppr[ijk][4*q+3]*ppr[ijk][4*q+3]; - return rs