VCAD

VCAD is a certified adaptive decomposition method for implicitly defined sets. It combines interval enclosures of first- and second-order derivatives, the Vertex Sufficiency Theorem (TVS), and Hilbert or Morton space-filling orderings to classify dyadic cells and estimate area or volume.

The repository contains the C++17 computational core, Python bindings, unit tests, benchmark programs, reproducibility scripts, and the data and figures used in the accompanying paper.

The Python import remains elipsus for compatibility with the current API and the experiments reported in the paper.

Method overview

For a region represented as

[ \Omega = {x \in D : G(x) \leq 0}, ]

the certified pipeline processes cells in two tracks at each refinement level:

1. Gate: cells on the ordinary refinement track are tested using G at the centre and an interval enclosure of the gradient over the complete cell.
2. TVS-VT: cells left undecided by the Gate use a Hessian enclosure to test whether every gradient component keeps a strict sign on that cell.
3. Track transition: cells passing TVS-VT enter the TVS-certified track. At the end of the level, both undecided tracks are subdivided into 2^n children.
4. TVS classification at the next level: children of a TVS-certified cell inherit its gradient-sign certificate. Their extrema occur at two opposite vertices, so the children are classified using those two values. Children still crossing the level set remain on this track; ordinary undecided cells return to the Gate.

The certified mode never silently falls back to the historical centre-Hessian approximation. A function without interval derivative extensions is rejected.

Here “certified mode” refers to the mathematical cell-wide derivative bounds used by the Gate and TVS-VT. The current implementation uses ordinary double arithmetic for point values, point gradients, and scalar reductions in the certificate tests. Consequently, it does not claim a fully end-to-end floating-point proof for the exceptional case in which a decision quantity lies within a few rounding units of zero. Closing that last numerical gap requires interval point evaluations and outward-rounded reductions, not an arbitrary tolerance.

Requirements

• A C++17 compiler
• CMake 3.16 or newer
• OpenMP
• Python 3.10 or newer for the Python interface
• NumPy and Matplotlib

Eigen, pybind11, and GoogleTest are found locally when available or fetched by CMake during configuration.

Python installation

From the repository root:

python3 -m pip install -e ".[dev,bench,viz3d]"

The bench and viz3d extras install the optional dependencies used by all benchmark utilities and 3-D visualizations. For the core library and tests, python3 -m pip install -e ".[dev]" is sufficient.

Minimal certified example:

import numpy as np
import elipsus

params = elipsus.DecomposerParams()
params.dimension = 2
params.initial_order = 3
params.refinement_depth = 8
params.domain_width = 2.5
params.epsilon = 0.0
params.integration_mode = True
params.certification = elipsus.CertificationMode.CERTIFIED
params.curve = elipsus.SpaceCurve.HILBERT

center = np.full(2, params.domain_width / 2)
disk = elipsus.ShiftedUnitSphere(2, center)
result = elipsus.Decomposer(params).run(disk)

lower = result.inside_volume
estimate = result.total_volume
upper = result.inside_volume + result.boundary_volume

print(f"volume bracket: [{lower:.8f}, {upper:.8f}]")
print(f"midpoint estimate: {estimate:.8f}")

For volume-only experiments, integration_mode=True avoids retaining IDs for cells already classified as inside and substantially reduces memory use.

C++ build and tests

cmake -S . -B build \
  -DCMAKE_BUILD_TYPE=Release \
  -DELIPSUS_BUILD_PYTHON=OFF \
  -DELIPSUS_BUILD_TESTS=ON \
  -DELIPSUS_BUILD_BENCHMARKS=ON
cmake --build build --parallel
ctest --test-dir build --output-on-failure

The benchmark build produces two main executables:

• bench_baselines: this work, its Morton ablation, and the uniform grid;
• bench_octree: the independent interval-octree baseline.

The octree implementation needed by the benchmark is vendored under benchmarks/octree, so no directory outside this repository is required.

Supported benchmark geometries

The built-in functions include shifted balls, a Bernoulli lemniscate, Rastrigin sublevel sets, a polar snowflake-shaped domain, intersecting spheres, a torus, and a three-dimensional snowflake tube. The Rastrigin and snowflake examples are general implicit sets rather than semi-algebraic sets.

Unit and shifted balls, Rastrigin, the torus, the polar snowflake, and the Bernoulli lemniscate provide the interval derivative extensions required by certified mode. IntersectingSpheres and SnowflakeTube3D currently provide only pointwise derivatives and must be used in legacy mode; certified mode rejects them explicitly.

Custom C++ functions implement FunctionInterface. Certified execution requires point evaluations of G and its gradient, together with interval extensions of the gradient and Hessian. The point Hessian remains part of the interface for the legacy mode and diagnostic use. Python subclasses can provide the pointwise interface, but certified interval extensions currently need a C++ implementation.

Boolean set operations

Decomposer.run_set supports unions and intersections of regions. Each defining function is certified independently before the classifications are combined. See the scripts under examples/ for complete uses of the API.

Reproducing the paper results

After building with ELIPSUS_BUILD_BENCHMARKS=ON:

# Main benchmark tables (computationally expensive)
python3 benchmarks/bench_sphere_table.py
python3 benchmarks/bench_2d_suite_accuracy.py

# Figures generated from the checked-in current measurements
python3 benchmarks/make_sphere_convergence_plot.py --from-csv
python3 benchmarks/make_sphere_volume_convergence_plot.py
python3 benchmarks/make_sphere_cells_processed_plot.py
python3 benchmarks/make_snowflake_volume_convergence_plot.py --from-csv

Omit --from-csv to rerun the certified decompositions before plotting. The exact benchmark drivers use Linux process-affinity and address-space controls; Linux is therefore required to reproduce the published timing protocol.

The full table runs use four pinned OpenMP threads and explicit per-process memory limits. The high-resolution uniform-grid cases can still take several minutes. Canonical result files and a mapping from each artifact to its generator are in paper/.

Repository layout

include/elipsus/   Public C++ headers and built-in implicit functions
src/               C++ implementation and Python bindings
python/elipsus/    Python package and visualization helpers
tests/             C++ and Python tests
examples/          Usage and visualization examples
benchmarks/        Benchmark executables and reproducibility scripts
paper/figures/     Current tables, measurements, and publication figures

Limitations

• Refinement creates 2^n children per cell, so time and memory grow rapidly with dimension.
• Certified functions currently require explicit analytical derivative and interval-extension routines.
• Natural interval extensions can overestimate derivative ranges and cause additional refinement.
• The certified derivative tests have a greater cost per processed cell than direct interval evaluation of G in the octree baseline.

Citation

Citation metadata is provided in CITATION.cff. For a paper, prefer the DOI of the archived release once it is deposited in Zenodo rather than citing the changing default branch.

License

This project is released under the MIT License.