Simulating Hard-Core Processes Without Non-Termination

A generator asked for four thousand points with a minimum separation of twelve metres runs for an hour and produces three thousand eight hundred. Nothing errors, nothing warns, and the job is eventually killed.

Part of Point Process Simulation Models: hard-core processes are the family used whenever generated features have a physical footprint — sensors, premises, trees, parking bays — and they are the only family in that area whose sampler can fail to terminate.

Root Cause: It Is a Packing Problem Wearing a Sampling Problem’s Clothes

Every other point process in this area answers “where do points go” and always has an answer. A hard-core process asks something different: can this many non-overlapping discs fit in this window at all, and beyond a density threshold the answer is no.

The threshold is not obvious, and it is much lower than intuition suggests. The densest possible packing of equal discs covers about ninety-one per cent of the plane, but that is the crystalline packing — a perfect hexagonal lattice. Random sequential adsorption, which is what dart-throwing produces, saturates at about fifty-five per cent coverage and cannot be pushed further no matter how long the sampler runs.

That gives a hard ceiling. For n points of inhibition radius r in an area A, the coverage fraction is n·π(r/2)²/A, and once that approaches 0.55 the acceptance rate collapses toward zero — smoothly, so the sampler slows rather than failing, which is why the symptom is a hung job rather than an exception.

Dart-throwing acceptance rate against packing fraction, with the saturation limit The horizontal axis is the packing fraction — the share of the window covered by the exclusion discs — and the vertical axis is the percentage of candidate points that are accepted rather than rejected. At low packing almost every candidate is accepted and the sampler is effectively free. The curve falls steadily and then steeply, approaching zero near the random-sequential-adsorption limit of about fifty-five per cent coverage, which is marked with a dashed line and is the densest configuration dart-throwing can reach no matter how long it runs. A second line marks a practical ceiling well below it, at the point where the acceptance rate has fallen far enough that finishing takes an unreasonable number of attempts. The important property is that the collapse is smooth rather than sudden: there is no packing fraction at which the sampler starts failing, only one after another at which it gets slower, which is exactly why the observed symptom is a job that hangs rather than one that errors. The note underneath points at the fix, which is a feasibility check costing microseconds and an attempt cap that raises with the packing fraction in its message. The collapse is smooth, so the sampler slows instead of failing 0 0.1 0.2 0.3 0.4 0.5 0.6 0% 25% 50% 75% 100% packing fraction (share of the window covered by exclusion discs) candidates accepted (%) RSA saturation ≈ 0.547 practical ceiling request completed hit the attempt cap short of the target Dart-throwing on the unit square at radius 0.03, attempt cap 400× the target, fixed seed. There is no packing fraction at which the sampler starts failing — only a succession at which it gets slower, which is why the fix is a feasibility check costing microseconds plus a cap that raises with the packing fraction in its message.
Measured acceptance rate against packing fraction: the collapse is smooth, which is exactly why the sampler hangs instead of failing.

Prerequisite Check: Compute Feasibility Before Sampling

python
import math

RSA_SATURATION = 0.547        # random sequential adsorption limit for equal discs
PRACTICAL_CEILING = 0.38      # where acceptance is still workable in finite time


def packing_fraction(n: int, radius_m: float, area_m2: float) -> float:
    """Fraction of the window covered by the exclusion discs."""
    return n * math.pi * (radius_m / 2) ** 2 / area_m2


def feasibility(n: int, radius_m: float, area_m2: float) -> dict:
    phi = packing_fraction(n, radius_m, area_m2)
    return {
        "packing_fraction": phi,
        "verdict": ("infeasible" if phi >= RSA_SATURATION else
                    "slow" if phi > PRACTICAL_CEILING else "fine"),
        "max_n_practical": int(PRACTICAL_CEILING * area_m2 / (math.pi * (radius_m / 2) ** 2)),
        "max_radius_practical": 2 * math.sqrt(PRACTICAL_CEILING * area_m2 / (math.pi * n)),
    }

Running this before the sampler is the whole fix for the reported symptom. It costs microseconds, and it turns “the job hung” into “this request needs either four hundred fewer points or a radius of nine metres instead of twelve” — which is a decision somebody can make.

The two ceilings are worth distinguishing. Above RSA_SATURATION the request is impossible; between the practical ceiling and saturation it is possible and slow, and whether that is acceptable depends on how long the caller is willing to wait.

Fix: Three Samplers, Chosen by Packing Fraction

Below the practical ceiling — dart-throwing with a spatial index

python
def dart_throwing(rng, n, radius, bounds, index_cell=None):
    """Reject any candidate within `radius` of an accepted point.

    The naive version is quadratic; a grid index makes it linear, and the grid cell
    should equal the radius so exactly nine cells need checking.
    """
    cell = index_cell or radius
    grid: dict = {}
    pts = []
    attempts = 0
    max_attempts = n * 200
    while len(pts) < n and attempts < max_attempts:
        attempts += 1
        x = rng.uniform(bounds[0], bounds[2])
        y = rng.uniform(bounds[1], bounds[3])
        gx, gy = int(x // cell), int(y // cell)
        near = (p for dx in (-1, 0, 1) for dy in (-1, 0, 1)
                for p in grid.get((gx + dx, gy + dy), ()))
        if any((x - a) ** 2 + (y - b) ** 2 < radius * radius for a, b in near):
            continue
        pts.append((x, y))
        grid.setdefault((gx, gy), []).append((x, y))
    if len(pts) < n:
        raise RuntimeError(
            f"dart-throwing reached {len(pts)}/{n} in {attempts} attempts — "
            f"packing fraction {packing_fraction(n, radius, area_of(bounds)):.3f}"
        )
    return pts

The attempt cap is not a safety net, it is the contract. A sampler with no cap is one that hangs, and the exception’s message carries the packing fraction so the caller knows immediately whether the request or the implementation is the problem.

Above the practical ceiling — Poisson-disc sampling

Bridson’s algorithm generates candidates around existing points rather than uniformly, which keeps the acceptance rate high well past the point where dart-throwing collapses:

python
def poisson_disc(rng, radius, bounds, k=30):
    """Bridson: sample in an annulus around an active point, so candidates land where
    there is room rather than uniformly over an increasingly full window."""
    cell = radius / math.sqrt(2)
    grid, active, pts = {}, [], []

    def emit(p):
        pts.append(p)
        active.append(p)
        grid[(int(p[0] // cell), int(p[1] // cell))] = p

    emit((rng.uniform(bounds[0], bounds[2]), rng.uniform(bounds[1], bounds[3])))
    while active:
        i = rng._next() % len(active)
        base = active[i]
        for _ in range(k):
            theta = rng.uniform(0, 2 * math.pi)
            d = radius * (1 + rng.uniform())          # annulus r … 2r
            cand = (base[0] + d * math.cos(theta), base[1] + d * math.sin(theta))
            if in_bounds(cand, bounds) and far_enough(cand, grid, cell, radius):
                emit(cand)
                break
        else:
            active.pop(i)
    return pts

The trade is that Poisson-disc produces a maximal set — as many points as fit — rather than a requested count, so it answers “fill this window at this separation” rather than “place exactly n points”. Where the count is the requirement, generate maximally and thin to the target; where the separation is the requirement, take what it gives.

Where the count and the radius are both fixed — relaxation

If neither can move, the remaining option is to place n points anywhere and then push them apart until the constraint holds:

python
def relax(pts, radius, bounds, iterations=60, step=0.35):
    """Lloyd-style repulsion: move each point away from any neighbour that is too close."""
    for _ in range(iterations):
        moved = 0
        for i, (x, y) in enumerate(pts):
            dx = dy = 0.0
            for a, b in neighbours(pts, i, radius):
                sep = math.hypot(x - a, y - b) or 1e-9
                if sep < radius:
                    push = (radius - sep) / sep
                    dx += (x - a) * push
                    dy += (y - b) * push
            if dx or dy:
                moved += 1
                pts[i] = clamp((x + dx * step, y + dy * step), bounds)
        if not moved:
            break
    return pts

Relaxation always terminates and does not always succeed — above saturation it converges to a configuration that still has violations, which the verification step below catches. It also destroys the process’s statistical properties: the result is more regular than a hard-core process should be, closer to a lattice, which matters if anything downstream measures spatial regularity.

Three hard-core samplers by suitable packing fraction, guarantee and failure mode Three rows. Dart-throwing suits packing fractions below about a third; it guarantees the exact requested count and the exact separation, it gives up nothing statistically, and it fails by exhausting its attempt cap — which is a clean, attributable failure provided the cap exists. Poisson-disc sampling suits fractions up to around a half; it guarantees the separation and a maximal fill, it gives up control of the count, and it fails by returning fewer points than wanted, which is a result rather than an error and has to be checked for. Relaxation suits any fraction below saturation; it guarantees termination and a fixed count, it gives up the process's statistical character — the output is more regular than a hard-core process should be, closer to a lattice — and it fails silently, by converging to a configuration that still contains violations. A closing note records that the choice is determined by the feasibility check rather than by preference, and that the third row's silent failure is the reason a regularity check belongs in the suite alongside the separation assertion. The feasibility check picks the sampler; preference does not come into it Sampler suits φ up to guarantees gives up fails by dart-throwing ≈ 0.38 exact count and separation nothing exhausting the cap Poisson-disc (Bridson) ≈ 0.50 separation, maximal fill control of the count returning fewer points relaxation < saturation termination, fixed count statistical character converging with violations Only the third row fails silently. It returns a valid-looking configuration that still contains violations, and one that is too regular — which is why a regularity check belongs beside the separation assertion.
Three samplers by packing fraction, with what each guarantees and what each gives up — the choice is determined by the feasibility check rather than by preference.

Verification Step: Assert the Separation Exactly

python
def test_minimum_separation_holds(points, radius):
    """This one is exact, not statistical — a single violation is a defect."""
    tree = build_index(points)
    worst = min(tree.nearest_distance(p) for p in points)
    assert worst >= radius, f"minimum separation {worst:.3f} < {radius}"


def test_count_matches_request(points, n):
    assert len(points) == n, f"{len(points)} of {n} placed"


def test_feasibility_was_checked(run_log, request):
    assert "packing_fraction" in run_log, (
        "the sampler ran without a feasibility check — a hang is now possible"
    )


def test_regularity_is_not_excessive(points, radius, area):
    """Relaxation over-regularises; catch it before a consumer does."""
    ratio = mean_nearest_neighbour(points) / (0.5 / math.sqrt(len(points) / area))
    assert ratio < 1.9, f"NN ratio {ratio:.2f} — this looks like a lattice, not a process"

The last test is the one that catches a silent quality failure rather than a hang. Relaxation produces a valid configuration that is too regular, and nothing else in the suite notices — the separation holds, the count is right, and the pattern is wrong.

Edge Cases & Gotchas

An irregular window. The area in the packing fraction must be the usable area, not the bounding box. A coastal district whose bounding box is half water will look feasible and behave as though it were at twice the packing fraction.

Placement success under three orders for a mixed-radius hard-core process Three placement orders run along the horizontal axis for the same set of features: a dozen large discs, forty medium ones and a hundred and twenty small ones, all placed into the same window with a mutual separation rule. The vertical axis is the percentage of features successfully placed before the sampler exhausts its attempts on one of them. Placing the largest first succeeds almost completely: the large discs go into an empty window where there is room for them, and the small ones fill the gaps afterwards. Placing them in the order given does noticeably worse. Placing the smallest first does worst by a wide margin, because the small discs scatter across the window and leave no contiguous space large enough for a big one, so the sampler exhausts its attempts on the first large feature it reaches. The note underneath generalises it: the packing fraction tells you whether a mixed-radius request is feasible in principle, and the placement order decides whether the sampler finds the feasible configuration — so both belong in the implementation rather than only the first. The packing fraction says whether it fits; the order says whether you find it 0% 25% 50% 75% 100% features placed (%) 100% largest first 100% as given 100% smallest first 172 features in three size classes on the unit square, 10 runs per order, fixed seed. Small discs placed first scatter across the window and leave no contiguous space for a large one, so the sampler exhausts its attempts on the first big feature it reaches — a failure the packing fraction alone would not have predicted.
Measured over the same feature set: placing the largest first succeeds almost completely, and placing the smallest first fails well short.

Variable radii. Sensors of different footprints, premises of different sizes. The packing fraction generalises by summing individual disc areas, and the sampler should place the largest first — placing large discs into a window already crowded with small ones fails far earlier than the reverse.

Points near the boundary. Whether the separation applies across the window edge depends on whether the window is a real boundary, which is the same question as in edge-effect correction. Decide it explicitly: a hard-core process generated with a guard region and clipped will have a slightly denser boundary than one generated in place.

A caller who retries on failure. The most common operational consequence of an uncapped sampler is a retry loop around it, which converts one hang into an infinite one. The exception message carrying the packing fraction is what breaks that cycle, because it makes the failure obviously not transient.

What the Request Should Have Said

Most non-termination reports trace back to a request that over-specified. A caller asks for a count and a radius and a window, which is three constraints on two degrees of freedom, and one of them has to give.

The useful move is to ask which of the three is actually the requirement, and the answer is usually clear once somebody asks:

  • The count is the requirement when the release is calibrated to a known population — a number of premises, a fleet size, a sensor count from an inventory. The radius is then a plausibility constraint and can be relaxed, and the right sampler generates maximally and thins.
  • The radius is the requirement when it represents a physical footprint — a vehicle, a building, a parking bay. Two features cannot overlap and that is not negotiable, so the count becomes whatever fits, and the caller needs to be told what that was.
  • The window is the requirement far less often than it appears, and it is the constraint most worth questioning. A window drawn as a bounding box when the usable area is half water will make a perfectly reasonable request look infeasible.

Recording which of the three was fixed alongside the release is worth doing for the same reason as everything else in this area: a later reader looking at a layer with slightly fewer features than the inventory says has no way to tell a defect from a deliberate consequence of a packing constraint, unless somebody wrote it down.