diskforge.cpp

#include <algorithm>
#include <array>
#include <atomic>
#include <bit>
#include <chrono>
#include <cfenv>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <sstream>
#include <stdexcept>
#include <string>
#include <string_view>
#include <system_error>
#include <utility>
#include <vector>

#include "diskforge/app.hpp"
#include "diskforge/cell_list.hpp"
#include "diskforge/cli.hpp"
#include "diskforge/foundation.hpp"
#include "diskforge/initialization.hpp"
#include "diskforge/manifest.hpp"
#include "diskforge/overlap.hpp"
#include "diskforge/reference_overlap.hpp"
#include "diskforge/validation.hpp"

#ifdef _OPENMP
#include <omp.h>
#endif

#include <unistd.h>

namespace diskforge {

using Clock = std::chrono::steady_clock;

struct Counters {
    std::uint64_t compression_steps = 0;
    std::uint64_t overlap_particle_scans = 0;
    std::uint64_t representatives = 0;
    std::uint64_t neighbor_lists = 0;
    std::uint64_t neighbor_visits = 0;
    std::uint64_t move_trials = 0;
    std::uint64_t accepted_moves = 0;
    std::uint64_t distance_checks = 0;
    std::uint64_t discovery_distance_checks = 0;
    std::uint64_t relaxation_distance_checks = 0;
    std::uint64_t scheduler_distance_checks = 0;
    std::uint64_t candidate_point_queries = 0;
    std::uint64_t event_candidate_pairs = 0;
    std::uint64_t event_events_scheduled = 0;
    std::uint64_t event_initial_events = 0;
    std::uint64_t event_local_events = 0;
    std::uint64_t event_rebuild_events = 0;
    std::uint64_t event_valid_events = 0;
    std::uint64_t event_stale_events = 0;
    std::uint64_t event_generation_bumps = 0;
    std::uint64_t event_dirty_particles = 0;
    std::uint64_t event_local_refreshes = 0;
    std::uint64_t event_rebuilds = 0;
    std::uint64_t event_rebuild_dropped_events = 0;
    std::uint64_t event_peak_queued = 0;
    std::uint64_t event_queued_end = 0;
    std::uint64_t sqrt_calls = 0;
    std::uint64_t temporary_vectors = 0;
    std::uint64_t cell_removals = 0;
    std::uint64_t cell_additions = 0;
    std::uint64_t radius_writes = 0;
    std::uint64_t anchor_writes = 0;
    std::uint64_t rng_draws = 0;

    bool operator==(const Counters&) const = default;
};

struct SimulationState {
    std::vector<Point> points;
    std::vector<Point> anchors;
    std::vector<double> legacy_radii;
    double radius;
    double box;
    PeriodicCellList cells;

    SimulationState(const InitialConfiguration& initial, double final_radius)
        : points(initial.points),
          anchors(initial.points),
          legacy_radii(initial.points.size(), initial.radius),
          radius(initial.radius),
          box(initial.box),
          cells(initial.box, 2.0 * final_radius, initial.points.size()) {
        if (!(radius <= final_radius)) [[unlikely]] {
            throw std::logic_error(
                "initial radius exceeds the requested final radius");
        }
        for (std::size_t particle = 0; particle < points.size(); ++particle) {
            cells.add(static_cast<int>(particle), points[particle]);
        }
    }
};

struct Workspace {
    static constexpr std::size_t kCandidateNeighborhoods = 9;

    std::vector<int> representatives;
    std::vector<int> neighbors;
    std::vector<int> dirty_particles;
    std::vector<unsigned char> dirty_flags;
    std::vector<unsigned char> overlap_flags;
    std::array<int, kCandidateNeighborhoods> candidate_cells{};
    std::array<std::vector<int>, kCandidateNeighborhoods>
        candidate_neighbors;
    std::vector<int> uncached_candidate_neighbors;
    std::size_t candidate_neighborhood_count = 0;

    explicit Workspace(std::size_t particle_count)
        : dirty_flags(particle_count, 0), overlap_flags(particle_count, 0) {
        representatives.reserve(particle_count);
        neighbors.reserve(particle_count);
        dirty_particles.reserve(particle_count);
    }

    void mark_dirty(int particle) {
        const std::size_t index = static_cast<std::size_t>(particle);
        if (dirty_flags[index] == 0) {
            dirty_flags[index] = 1;
            dirty_particles.push_back(particle);
        }
    }

    void write_dirty_anchors(SimulationState& state, Counters& counters) {
        for (const int particle : dirty_particles) {
            const std::size_t index = static_cast<std::size_t>(particle);
            state.anchors[index] = state.points[index];
            ++counters.anchor_writes;
        }
    }

    void clear_dirty() {
        for (const int particle : dirty_particles) {
            dirty_flags[static_cast<std::size_t>(particle)] = 0;
        }
        dirty_particles.clear();
    }

    void update_dirty_anchors(SimulationState& state, Counters& counters) {
        write_dirty_anchors(state, counters);
        clear_dirty();
    }

    void clear_candidate_neighborhoods() {
        candidate_neighborhood_count = 0;
    }

    const std::vector<int>& candidate_neighborhood(
        const PeriodicCellList& cells, int cell) {
        for (std::size_t index = 0;
             index < candidate_neighborhood_count; ++index) {
            if (candidate_cells[index] == cell) {
                return candidate_neighbors[index];
            }
        }
        if (candidate_neighborhood_count >= candidate_cells.size()) [[unlikely]] {
            // Exact boundary rounding can expose more than the mathematical
            // 3x3 set. Preserve correctness with an uncached query rather than
            // turning the optimization's common-case bound into a failure.
            cells.gather_cell_neighbors(cell, uncached_candidate_neighbors);
            return uncached_candidate_neighbors;
        }
        const std::size_t index = candidate_neighborhood_count++;
        candidate_cells[index] = cell;
        cells.gather_cell_neighbors(cell, candidate_neighbors[index]);
        return candidate_neighbors[index];
    }
};

bool has_overlap_legacy(const SimulationState& state,
                        int particle,
                        Counters& counters,
                        std::uint64_t& phase_distance_checks) {
    const std::vector<int> neighbors = state.cells.neighbors_allocating(particle);
    ++counters.temporary_vectors;
    for (const int neighbor : neighbors) {
        if (neighbor == particle) {
            continue;
        }
        ++counters.distance_checks;
        ++phase_distance_checks;
        ++counters.sqrt_calls;
        const double distance =
            std::sqrt(squared_distance(state.points[static_cast<std::size_t>(particle)],
                                       state.points[static_cast<std::size_t>(neighbor)], state.box));
        if (distance <
            state.legacy_radii[static_cast<std::size_t>(particle)] +
                state.legacy_radii[static_cast<std::size_t>(neighbor)]) {
            return true;
        }
    }
    return false;
}

bool has_overlap_at_point_fast(const SimulationState& state,
                               const Point& candidate,
                               const std::vector<int>& neighbors,
                               StrictOverlapCutoff cutoff,
                               Counters& counters) {
    for (const int neighbor : neighbors) {
        ++counters.distance_checks;
        ++counters.relaxation_distance_checks;
        if (strict_overlap_from_squared(
                squared_distance(candidate, state.points[static_cast<std::size_t>(neighbor)],
                                 state.box),
                cutoff)) {
            return true;
        }
    }
    return false;
}

void discover_overlaps_legacy(const SimulationState& state,
                              Workspace& workspace,
                              Counters& counters) {
    workspace.representatives.clear();
    for (std::size_t particle = 0; particle < state.points.size(); ++particle) {
        ++counters.overlap_particle_scans;
        if (has_overlap_legacy(state, static_cast<int>(particle), counters,
                               counters.discovery_distance_checks)) {
            workspace.representatives.push_back(static_cast<int>(particle));
        }
    }
    counters.representatives += workspace.representatives.size();
}

void discover_overlaps_fast(const SimulationState& state,
                            Workspace& workspace,
                            StrictOverlapCutoff cutoff,
                            Counters& counters) {
    workspace.representatives.clear();
    std::fill(workspace.overlap_flags.begin(), workspace.overlap_flags.end(), 0);
    std::uint64_t distance_checks = 0;
    for (std::size_t cell = 0; cell < state.cells.cell_count(); ++cell) {
        state.cells.for_each_unordered_neighbor_pair_owned_by_cell(
            static_cast<int>(cell), [&](int particle, int neighbor) {
                ++distance_checks;
                if (strict_overlap_from_squared(
                        squared_distance(
                            state.points[static_cast<std::size_t>(particle)],
                            state.points[static_cast<std::size_t>(neighbor)],
                            state.box),
                        cutoff)) [[unlikely]] {
                    workspace.overlap_flags[
                        static_cast<std::size_t>(particle)] = 1;
                    workspace.overlap_flags[static_cast<std::size_t>(neighbor)] = 1;
                }
            });
    }
    counters.overlap_particle_scans += state.points.size();
    counters.distance_checks += distance_checks;
    counters.discovery_distance_checks += distance_checks;
    for (std::size_t particle = 0; particle < workspace.overlap_flags.size(); ++particle) {
        if (workspace.overlap_flags[particle] != 0) {
            workspace.representatives.push_back(static_cast<int>(particle));
        }
    }
    counters.representatives += workspace.representatives.size();
}

int discover_overlaps_omp(const SimulationState& state,
                          Workspace& workspace,
                          StrictOverlapCutoff cutoff,
                          Counters& counters,
                          int requested_threads) {
    workspace.representatives.clear();
    std::fill(workspace.overlap_flags.begin(), workspace.overlap_flags.end(), 0);
    std::uint64_t distance_checks = 0;
    int threads_used = 1;

#ifdef _OPENMP
    static const int runtime_default_team_size = omp_get_max_threads();
    const int previous_dynamic = omp_get_dynamic();
    const int team_size =
        requested_threads > 0 ? requested_threads : runtime_default_team_size;
    omp_set_dynamic(0);
#pragma omp parallel num_threads(team_size) reduction(+ : distance_checks)
    {
#pragma omp single
        { threads_used = omp_get_num_threads(); }
#pragma omp for schedule(static)
        for (std::ptrdiff_t cell = 0;
             cell < static_cast<std::ptrdiff_t>(state.cells.cell_count()); ++cell) {
            state.cells.for_each_unordered_neighbor_pair_owned_by_cell(
                static_cast<int>(cell), [&](int particle, int neighbor) {
                    ++distance_checks;
                    if (strict_overlap_from_squared(
                            squared_distance(
                                state.points[static_cast<std::size_t>(particle)],
                                state.points[static_cast<std::size_t>(neighbor)],
                                state.box),
                            cutoff)) [[unlikely]] {
                        std::atomic_ref<unsigned char>(
                            workspace.overlap_flags[
                                static_cast<std::size_t>(particle)])
                            .store(1, std::memory_order_relaxed);
                        std::atomic_ref<unsigned char>(
                            workspace.overlap_flags[
                                static_cast<std::size_t>(neighbor)])
                            .store(1, std::memory_order_relaxed);
                    }
                });
        }
    }
    omp_set_dynamic(previous_dynamic);
#else
    (void)requested_threads;
    for (std::size_t cell = 0; cell < state.cells.cell_count(); ++cell) {
        state.cells.for_each_unordered_neighbor_pair_owned_by_cell(
            static_cast<int>(cell), [&](int particle, int neighbor) {
                ++distance_checks;
                if (strict_overlap_from_squared(
                        squared_distance(
                            state.points[static_cast<std::size_t>(particle)],
                            state.points[static_cast<std::size_t>(neighbor)],
                            state.box),
                        cutoff)) [[unlikely]] {
                    workspace.overlap_flags[
                        static_cast<std::size_t>(particle)] = 1;
                    workspace.overlap_flags[static_cast<std::size_t>(neighbor)] = 1;
                }
            });
    }
#endif

    counters.overlap_particle_scans += state.points.size();
    counters.distance_checks += distance_checks;
    counters.discovery_distance_checks += distance_checks;
    // Serial compaction is the reproducibility boundary: output is always index ordered.
    for (std::size_t particle = 0; particle < workspace.overlap_flags.size(); ++particle) {
        if (workspace.overlap_flags[particle] != 0) {
            workspace.representatives.push_back(static_cast<int>(particle));
        }
    }
    counters.representatives += workspace.representatives.size();
    return threads_used;
}

void relocate_until_clear_legacy(SimulationState& state,
                                 int particle,
                                 ReproducibleRng& rng,
                                 const Parameters& parameters,
                                 Counters& counters) {
    const std::size_t index = static_cast<std::size_t>(particle);
    const double displacement_extent = 2.0 * state.radius;
    for (std::size_t trial = 0; trial < parameters.max_trials_per_move; ++trial) {
        ++counters.move_trials;
        state.cells.remove(particle);
        ++counters.cell_removals;
        const Point displacement =
            rng.symmetric_pair(displacement_extent);
        state.points[index] = wrap_point(
            {state.anchors[index].x + displacement.x,
             state.anchors[index].y + displacement.y},
            state.box);
        counters.rng_draws += 2;
        state.cells.add(particle, state.points[index]);
        ++counters.cell_additions;

        if (!has_overlap_legacy(state, particle, counters,
                                counters.relaxation_distance_checks)) {
            ++counters.accepted_moves;
            return;
        }
    }
    throw std::runtime_error("relocation trial limit exhausted for particle " +
                             std::to_string(particle));
}

void relocate_until_clear_fast(SimulationState& state,
                               int particle,
                               ReproducibleRng& rng,
                               const Parameters& parameters,
                               Workspace& workspace,
                               StrictOverlapCutoff cutoff,
                               Counters& counters) {
    const std::size_t index = static_cast<std::size_t>(particle);
    const double displacement_extent = 2.0 * state.radius;
    const Point previous_position = state.points[index];
    state.cells.remove(particle);
    ++counters.cell_removals;
    workspace.clear_candidate_neighborhoods();

    for (std::size_t trial = 0; trial < parameters.max_trials_per_move; ++trial) {
        ++counters.move_trials;
        const Point displacement =
            rng.symmetric_pair(displacement_extent);
        const Point candidate = wrap_point(
            {state.anchors[index].x + displacement.x,
             state.anchors[index].y + displacement.y},
            state.box);
        counters.rng_draws += 2;
        ++counters.candidate_point_queries;

        const int candidate_cell = state.cells.cell_for_point(candidate);
        const std::vector<int>& candidate_neighbors =
            workspace.candidate_neighborhood(state.cells, candidate_cell);
        if (!has_overlap_at_point_fast(
                state, candidate, candidate_neighbors, cutoff, counters)) {
            state.points[index] = candidate;
            state.cells.add(particle, candidate);
            ++counters.cell_additions;
            ++counters.accepted_moves;
            return;
        }
    }

    state.points[index] = previous_position;
    state.cells.add(particle, previous_position);
    ++counters.cell_additions;
    throw std::runtime_error("relocation trial limit exhausted for particle " +
                             std::to_string(particle));
}

void relax_representatives(SimulationState& state,
                           Workspace& workspace,
                           ReproducibleRng& rng,
                           const Parameters& parameters,
                           Engine engine,
                           StrictOverlapCutoff cutoff,
                           Counters& counters) {
    for (const int representative : workspace.representatives) {
        if (engine == Engine::scan_legacy) {
            std::vector<int> neighbors = state.cells.neighbors_allocating(representative);
            ++counters.temporary_vectors;
            ++counters.neighbor_lists;
            counters.neighbor_visits += neighbors.size();
            for (const int particle : neighbors) {
                relocate_until_clear_legacy(state, particle, rng, parameters, counters);
            }
        } else {
            state.cells.gather_neighbors(representative, workspace.neighbors);
            ++counters.neighbor_lists;
            counters.neighbor_visits += workspace.neighbors.size();
            for (const int particle : workspace.neighbors) {
                workspace.mark_dirty(particle);
                relocate_until_clear_fast(
                    state, particle, rng, parameters, workspace, cutoff,
                    counters);
            }
        }
    }
}

enum class StepVerification {
    none,
    discovery,
    discovery_and_clear,
};

struct RunControls {
    StepVerification verification = StepVerification::none;
    bool force_event_rebuilds = false;
};

struct VerificationStats {
    std::uint64_t discovery_scans = 0;
    std::uint64_t post_relaxation_scans = 0;
    std::uint64_t pairs_checked = 0;
};

struct RunResult {
    Engine engine = Engine::scan_fast;
    std::vector<Point> points;
    double radius = 0.0;
    double box = 0.0;
    double seconds = 0.0;
    double setup_seconds = 0.0;
    double discovery_seconds = 0.0;
    double relaxation_seconds = 0.0;
    double scheduler_seconds = 0.0;
    int threads_used = 1;
    Counters counters;
    OracleResult oracle;
    VerificationStats verification;
    std::uint64_t state_checksum = 0;
};

enum class EventScheduleOrigin {
    initial,
    local,
    rebuild,
};

struct PairEvent {
    std::uint32_t a = 0;
    std::uint32_t b = 0;
    std::uint64_t scheduled_after_step = 0;
};

static_assert(sizeof(PairEvent) == 16);

std::size_t first_strict_event_step(const std::vector<double>& diameters,
                                    double distance,
                                    std::size_t current_step) {
    if (current_step + 1U >= diameters.size()) {
        return diameters.size();
    }
    const auto search_begin =
        diameters.begin() + static_cast<std::ptrdiff_t>(current_step + 1U);
    const auto event_position =
        std::upper_bound(search_begin, diameters.end(), distance);
    return static_cast<std::size_t>(event_position - diameters.begin());
}

class EventCalendar {
  public:
    EventCalendar(const SimulationState& state,
                  const Parameters& parameters,
                  Counters& counters)
        : radii_(build_radii(state.radius, parameters)),
          diameters_(radii_.size()),
          events_at_step_(radii_.size()),
          last_dirty_step_(state.points.size(), 0),
          endpoint_bits_((state.points.size() + 63U) / 64U, 0) {
        const double final_diameter = 2.0 * parameters.final_radius;
        if (state.cells.cell_size() < final_diameter) {
            throw std::logic_error(
                "event scheduler requires cell width at least the final diameter");
        }
        std::transform(radii_.begin(), radii_.end(), diameters_.begin(),
                       [](double radius) { return radius + radius; });
        enumerate_all_pairs(state, 0, EventScheduleOrigin::initial, counters);
        live_count_at_last_rebuild_ = queued_events_;
        counters.event_queued_end = queued_events_;
    }

    std::size_t step_count() const {
        return radii_.size() - 1U;
    }

    double radius_at(std::size_t step) const {
        return radii_.at(step);
    }

    void discover_representatives(const SimulationState& state,
                                  std::size_t step,
                                  Workspace& workspace,
                                  StrictOverlapCutoff cutoff,
                                  Counters& counters) {
        if (step == 0 || step >= events_at_step_.size()) {
            throw std::out_of_range("event calendar step is out of range");
        }

        std::vector<PairEvent> bucket =
            std::move(events_at_step_[step]);
        if (queued_events_ < bucket.size()) {
            throw std::logic_error("event queue accounting underflow");
        }
        queued_events_ -= bucket.size();

        for (const PairEvent& event : bucket) {
            const std::size_t a = event.a;
            const std::size_t b = event.b;
            if (last_dirty_step_[a] > event.scheduled_after_step ||
                last_dirty_step_[b] > event.scheduled_after_step) {
                ++stale_since_rebuild_;
                ++counters.event_stale_events;
                continue;
            }
#ifndef NDEBUG
            if (!strict_overlap_from_squared(
                    squared_distance(state.points[a], state.points[b], state.box),
                    cutoff)) {
                throw std::logic_error("live event is not overlapping at its scheduled step");
            }
#else
            (void)cutoff;
#endif
            endpoint_bits_[a / 64U] |= 1ULL << (a % 64U);
            endpoint_bits_[b / 64U] |= 1ULL << (b % 64U);
            ++valid_since_rebuild_;
            ++counters.event_valid_events;
        }

        workspace.representatives.clear();
        for (std::size_t block_index = 0; block_index < endpoint_bits_.size(); ++block_index) {
            std::uint64_t block =
                std::exchange(endpoint_bits_[block_index], 0);
            while (block != 0) {
                const unsigned offset = std::countr_zero(block);
                const std::size_t particle = block_index * 64U + offset;
                if (particle < state.points.size()) {
                    workspace.representatives.push_back(static_cast<int>(particle));
                }
                block &= block - 1U;
            }
        }
        counters.representatives += workspace.representatives.size();
        counters.event_queued_end = queued_events_;
    }

    void finish_step(SimulationState& state,
                     std::size_t step,
                     Workspace& workspace,
                     Counters& counters,
                     bool force_rebuild) {
        const std::size_t dirty_count = workspace.dirty_particles.size();
        counters.event_dirty_particles += dirty_count;
        for (const int particle : workspace.dirty_particles) {
            last_dirty_step_[static_cast<std::size_t>(particle)] = step;
            ++counters.event_generation_bumps;
        }
        workspace.write_dirty_anchors(state, counters);

        if (step < step_count()) {
            const bool rebuild =
                force_rebuild || dirty_count > state.points.size() / 3U ||
                queued_events_ >
                    std::max(4ULL * static_cast<std::uint64_t>(state.points.size()),
                             3ULL * live_count_at_last_rebuild_) ||
                stale_since_rebuild_ >
                    std::max(static_cast<std::uint64_t>(state.points.size()),
                             std::uint64_t{2} * valid_since_rebuild_);
            if (rebuild) {
                rebuild_future(state, step, counters);
            } else if (dirty_count != 0) {
                ++counters.event_local_refreshes;
                enumerate_dirty_pairs(state, step, workspace, counters);
            }
        }

        workspace.clear_dirty();
        counters.event_queued_end = queued_events_;
    }

  private:
    static std::vector<double> build_radii(double initial_radius,
                                           const Parameters& parameters) {
        std::vector<double> radii;
        radii.reserve(std::min<std::size_t>(parameters.max_compression_steps, 4095U) + 1U);
        radii.push_back(initial_radius);
        while (radii.back() < parameters.final_radius) {
            if (radii.size() - 1U >= parameters.max_compression_steps) {
                throw std::runtime_error("compression step limit exhausted");
            }
            radii.push_back(next_compression_radius(radii.back(), parameters));
        }
        return radii;
    }

    void schedule_pair(const SimulationState& state,
                       int first,
                       int second,
                       std::size_t current_step,
                       EventScheduleOrigin origin,
                       Counters& counters) {
        if (first == second || current_step >= step_count()) {
            return;
        }
        const int a = std::min(first, second);
        const int b = std::max(first, second);
        ++counters.event_candidate_pairs;
        ++counters.distance_checks;
        ++counters.scheduler_distance_checks;
        const double distance =
            std::sqrt(squared_distance(state.points[static_cast<std::size_t>(a)],
                                       state.points[static_cast<std::size_t>(b)], state.box));
        ++counters.sqrt_calls;
        const std::size_t event_step =
            first_strict_event_step(diameters_, distance, current_step);
        if (event_step == diameters_.size()) {
            return;
        }

        events_at_step_[event_step].push_back(
            {static_cast<std::uint32_t>(a), static_cast<std::uint32_t>(b),
             static_cast<std::uint64_t>(current_step)});
        ++queued_events_;
        ++counters.event_events_scheduled;
        switch (origin) {
            case EventScheduleOrigin::initial:
                ++counters.event_initial_events;
                break;
            case EventScheduleOrigin::local:
                ++counters.event_local_events;
                break;
            case EventScheduleOrigin::rebuild:
                ++counters.event_rebuild_events;
                break;
        }
        counters.event_peak_queued = std::max(counters.event_peak_queued, queued_events_);
    }

    void enumerate_all_pairs(const SimulationState& state,
                             std::size_t current_step,
                             EventScheduleOrigin origin,
                             Counters& counters) {
        for (std::size_t particle = 0; particle < state.points.size(); ++particle) {
            state.cells.for_each_neighbor_until(static_cast<int>(particle), [&](int neighbor) {
                if (neighbor > static_cast<int>(particle)) {
                    schedule_pair(state, static_cast<int>(particle), neighbor, current_step,
                                  origin, counters);
                }
                return true;
            });
        }
    }

    void enumerate_dirty_pairs(const SimulationState& state,
                               std::size_t current_step,
                               const Workspace& workspace,
                               Counters& counters) {
        for (const int particle : workspace.dirty_particles) {
            state.cells.for_each_neighbor_until(particle, [&](int neighbor) {
                if (neighbor == particle) {
                    return true;
                }
                if (workspace.dirty_flags[static_cast<std::size_t>(neighbor)] != 0 &&
                    neighbor < particle) {
                    return true;
                }
                schedule_pair(state, particle, neighbor, current_step,
                              EventScheduleOrigin::local, counters);
                return true;
            });
        }
    }

    void rebuild_future(const SimulationState& state,
                        std::size_t current_step,
                        Counters& counters) {
        ++counters.event_rebuilds;
        counters.event_rebuild_dropped_events += queued_events_;
        for (std::size_t step = current_step + 1U; step < events_at_step_.size(); ++step) {
            events_at_step_[step].clear();
        }
        queued_events_ = 0;
        enumerate_all_pairs(state, current_step, EventScheduleOrigin::rebuild, counters);
        live_count_at_last_rebuild_ = queued_events_;
        stale_since_rebuild_ = 0;
        valid_since_rebuild_ = 0;
    }

    std::vector<double> radii_;
    std::vector<double> diameters_;
    std::vector<std::vector<PairEvent>> events_at_step_;
    std::vector<std::uint64_t> last_dirty_step_;
    std::vector<std::uint64_t> endpoint_bits_;
    std::uint64_t queued_events_ = 0;
    std::uint64_t live_count_at_last_rebuild_ = 0;
    std::uint64_t stale_since_rebuild_ = 0;
    std::uint64_t valid_since_rebuild_ = 0;
};

void verify_discovery_with_reference(const SimulationState& state,
                                     const std::vector<int>& representatives,
                                     std::size_t step,
                                     Engine engine,
                                     VerificationStats& stats) {
    const ReferenceOverlapScan reference =
        scan_overlaps_reference(state.points, state.radius, state.box);
    ++stats.discovery_scans;
    stats.pairs_checked += reference.pairs_checked;
    if (reference.representatives != representatives) {
        throw std::logic_error(
            std::string(engine_name(engine)) +
            "/independent-reference representative mismatch at compression step " +
            std::to_string(step));
    }
}

void verify_clear_with_reference(const SimulationState& state,
                                 std::size_t step,
                                 Engine engine,
                                 VerificationStats& stats) {
    const ReferenceOverlapScan reference =
        scan_overlaps_reference(state.points, state.radius, state.box);
    ++stats.post_relaxation_scans;
    stats.pairs_checked += reference.pairs_checked;
    if (reference.overlapping_pairs != 0) {
        throw std::logic_error(
            std::string(engine_name(engine)) +
            "/independent-reference found overlaps after relaxation at compression step " +
            std::to_string(step));
    }
}

void validate_final_state_fast(const SimulationState& state) {
    if (!(state.radius > 0.0) || !std::isfinite(state.radius) ||
        !(state.box > 0.0) || !std::isfinite(state.box) ||
        state.radius + state.radius > state.box) {
        throw std::logic_error(
            "final state violates finite periodic geometry invariants");
    }
    if (!state.cells.validate(state.points)) {
        throw std::logic_error(
            "final state violates cell-list membership invariants");
    }
    const StrictOverlapCutoff cutoff =
        make_strict_overlap_cutoff(state.radius + state.radius);
    std::uint64_t distance_checks = 0;
    for (std::size_t particle = 0; particle < state.points.size(); ++particle) {
        if (has_strict_overlap(
                state.points, state.cells, static_cast<int>(particle),
                cutoff, state.box, distance_checks)) {
            throw std::logic_error(
                "final state contains an overlap under the source predicate");
        }
    }
}

RunResult run_engine(const InitialConfiguration& initial,
                     const Parameters& parameters,
                     Engine engine,
                     RunControls controls = {}) {
    const auto start = Clock::now();
    SimulationState state(initial, parameters.final_radius);
    ReproducibleRng rng = initial.rng_after_initialization;
    Workspace workspace(state.points.size());
    Counters counters;
    VerificationStats verification_stats;
    StepVerification verification = controls.verification;
#ifndef NDEBUG
    if (verification == StepVerification::none && state.points.size() <= 128U) {
        verification = StepVerification::discovery_and_clear;
    }
#endif
    const double setup_seconds =
        std::chrono::duration<double>(Clock::now() - start).count();
    int threads_used = 1;
    double discovery_seconds = 0.0;
    double relaxation_seconds = 0.0;
    double scheduler_seconds = 0.0;

    if (engine == Engine::event) {
        const auto scheduler_start = Clock::now();
        EventCalendar calendar(state, parameters, counters);
        scheduler_seconds +=
            std::chrono::duration<double>(Clock::now() - scheduler_start).count();
        for (std::size_t step = 1; step <= calendar.step_count(); ++step) {
            state.radius = calendar.radius_at(step);
            ++counters.compression_steps;

            const auto discovery_start = Clock::now();
            const StrictOverlapCutoff cutoff =
                make_strict_overlap_cutoff(
                    state.radius + state.radius, &counters.sqrt_calls);
            calendar.discover_representatives(
                state, step, workspace, cutoff, counters);
            discovery_seconds +=
                std::chrono::duration<double>(Clock::now() - discovery_start).count();
            if (verification != StepVerification::none) [[unlikely]] {
                verify_discovery_with_reference(
                    state, workspace.representatives, step, engine, verification_stats);
            }

            const auto relaxation_start = Clock::now();
            relax_representatives(
                state, workspace, rng, parameters, engine, cutoff, counters);
            relaxation_seconds +=
                std::chrono::duration<double>(Clock::now() - relaxation_start).count();
            if (verification == StepVerification::discovery_and_clear) [[unlikely]] {
                verify_clear_with_reference(state, step, engine, verification_stats);
            }

            const auto refresh_start = Clock::now();
            calendar.finish_step(
                state, step, workspace, counters, controls.force_event_rebuilds);
            scheduler_seconds +=
                std::chrono::duration<double>(Clock::now() - refresh_start).count();

#ifndef NDEBUG
            if (!state.cells.validate(state.points)) {
                throw std::logic_error("cell-list invariant failed during compression");
            }
#endif
        }
    } else {
        while (state.radius < parameters.final_radius) {
            if (counters.compression_steps >= parameters.max_compression_steps) {
                throw std::runtime_error("compression step limit exhausted");
            }
            state.radius =
                next_compression_radius(state.radius, parameters);
            ++counters.compression_steps;

            const auto discovery_start = Clock::now();
            StrictOverlapCutoff cutoff;
            if (engine == Engine::scan_legacy) {
                std::fill(state.legacy_radii.begin(), state.legacy_radii.end(), state.radius);
                counters.radius_writes += state.legacy_radii.size();
                discover_overlaps_legacy(state, workspace, counters);
            } else {
                cutoff = make_strict_overlap_cutoff(
                    state.radius + state.radius, &counters.sqrt_calls);
                if (engine == Engine::scan_fast) {
                    discover_overlaps_fast(
                        state, workspace, cutoff, counters);
                } else {
                    threads_used = discover_overlaps_omp(
                        state, workspace, cutoff, counters,
                        parameters.threads);
                }
            }
            discovery_seconds +=
                std::chrono::duration<double>(Clock::now() - discovery_start).count();
            if (verification != StepVerification::none) [[unlikely]] {
                verify_discovery_with_reference(
                    state, workspace.representatives,
                    static_cast<std::size_t>(counters.compression_steps), engine,
                    verification_stats);
            }

            const auto relaxation_start = Clock::now();
            relax_representatives(
                state, workspace, rng, parameters, engine, cutoff, counters);

            if (engine == Engine::scan_legacy) {
                state.anchors = state.points;
                counters.anchor_writes += state.anchors.size();
            } else {
                workspace.update_dirty_anchors(state, counters);
            }
            relaxation_seconds +=
                std::chrono::duration<double>(Clock::now() - relaxation_start).count();
            if (verification == StepVerification::discovery_and_clear) [[unlikely]] {
                verify_clear_with_reference(
                    state, static_cast<std::size_t>(counters.compression_steps), engine,
                    verification_stats);
            }

#ifndef NDEBUG
            if (!state.cells.validate(state.points)) {
                throw std::logic_error("cell-list invariant failed during compression");
            }
#endif
        }
    }

    const auto compression_end = Clock::now();
    validate_final_state_fast(state);

    RunResult result;
    result.engine = engine;
    result.seconds =
        std::chrono::duration<double>(compression_end - start).count();
    result.setup_seconds = setup_seconds;
    result.discovery_seconds = discovery_seconds;
    result.relaxation_seconds = relaxation_seconds;
    result.scheduler_seconds = scheduler_seconds;
    result.points = std::move(state.points);
    result.radius = state.radius;
    result.box = state.box;
    result.threads_used = threads_used;
    result.counters = counters;
    result.verification = verification_stats;
    if (!parameters.skip_oracle) {
        result.oracle = final_oracle(
            result.points, result.radius, result.box,
            parameters.oracle_threads);
    }
    result.state_checksum = checksum(result.points, result.radius, result.box);
    return result;
}

double maximum_coordinate_difference(const RunResult& lhs, const RunResult& rhs) {
    if (lhs.points.size() != rhs.points.size()) {
        return std::numeric_limits<double>::infinity();
    }
    double maximum = 0.0;
    for (std::size_t particle = 0; particle < lhs.points.size(); ++particle) {
        maximum = std::max(maximum, std::fabs(lhs.points[particle].x - rhs.points[particle].x));
        maximum = std::max(maximum, std::fabs(lhs.points[particle].y - rhs.points[particle].y));
    }
    return maximum;
}

bool positions_bitwise_equal(const RunResult& lhs, const RunResult& rhs) {
    if (lhs.points.size() != rhs.points.size() ||
        std::bit_cast<std::uint64_t>(lhs.radius) !=
            std::bit_cast<std::uint64_t>(rhs.radius) ||
        std::bit_cast<std::uint64_t>(lhs.box) !=
            std::bit_cast<std::uint64_t>(rhs.box)) {
        return false;
    }
    for (std::size_t particle = 0; particle < lhs.points.size(); ++particle) {
        if (std::bit_cast<std::uint64_t>(lhs.points[particle].x) !=
                std::bit_cast<std::uint64_t>(rhs.points[particle].x) ||
            std::bit_cast<std::uint64_t>(lhs.points[particle].y) !=
                std::bit_cast<std::uint64_t>(rhs.points[particle].y)) {
            return false;
        }
    }
    return true;
}

std::string render_configuration(const RunResult& result,
                                 const Parameters& parameters) {
    std::ostringstream output;
    output << "# diskforge C++20 configuration\n"
           << "# engine " << engine_name(result.engine) << " seed " << parameters.seed << " N "
           << result.points.size() << " box " << std::setprecision(17) << result.box << " radius "
           << result.radius << "\n"
           << "# columns: particle x y radius\n";
    for (std::size_t particle = 0; particle < result.points.size(); ++particle) {
        output << particle << ' ' << std::setprecision(17) << result.points[particle].x << ' '
               << result.points[particle].y << ' ' << result.radius << '\n';
    }
    if (!output) {
        throw std::runtime_error("failed to render configuration output");
    }
    return output.str();
}

struct StagedConfiguration {
    std::filesystem::path target;
    std::filesystem::path temporary;
    std::string artifact = "output";
};

StagedConfiguration stage_configuration(const RunResult& result,
                                         const Parameters& parameters,
                                         const std::string& path) {
    static std::atomic<std::uint64_t> sequence{0};
    const std::filesystem::path target(path);
    std::error_code status_error;
    const std::filesystem::file_status target_status =
        std::filesystem::status(target, status_error);
    if (status_error &&
        status_error != std::errc::no_such_file_or_directory) {
        throw std::runtime_error(
            "cannot inspect output path " + path + ": " +
            status_error.message());
    }
    if (!status_error && std::filesystem::exists(target_status) &&
        !std::filesystem::is_regular_file(target_status)) {
        throw std::runtime_error(
            "output path exists and is not a regular file: " + path);
    }
    std::filesystem::path temporary = target;
    temporary += ".tmp." + std::to_string(static_cast<long long>(::getpid())) +
                 "." + std::to_string(sequence.fetch_add(1));

    const std::string content = render_configuration(result, parameters);
    if (content.size() >
        static_cast<std::size_t>(
            std::numeric_limits<std::streamsize>::max())) {
        throw std::runtime_error(
            "configuration output exceeds stream size limits");
    }
    std::ofstream output(
        temporary, std::ios::binary | std::ios::out | std::ios::trunc);
    if (!output) {
        throw std::runtime_error(
            "cannot stage output file " + path + " via " +
            temporary.string());
    }
    output.write(content.data(), static_cast<std::streamsize>(content.size()));
    output.flush();
    if (!output) {
        output.close();
        std::error_code ignored;
        std::filesystem::remove(temporary, ignored);
        throw std::runtime_error(
            "failed while staging output file: " + path);
    }
    output.close();
    if (output.fail()) {
        std::error_code ignored;
        std::filesystem::remove(temporary, ignored);
        throw std::runtime_error(
            "failed to close staged output file: " + path);
    }
    return {target, temporary};
}

StagedConfiguration stage_manifest(const std::string& content,
                                    const std::string& path) {
    static std::atomic<std::uint64_t> sequence{0};
    const std::filesystem::path target(path);
    std::error_code status_error;
    const std::filesystem::file_status target_status =
        std::filesystem::status(target, status_error);
    if (status_error &&
        status_error != std::errc::no_such_file_or_directory) {
        throw std::runtime_error(
            "cannot inspect manifest path " + path + ": " +
            status_error.message());
    }
    if (!status_error && std::filesystem::exists(target_status) &&
        !std::filesystem::is_regular_file(target_status)) {
        throw std::runtime_error(
            "manifest path exists and is not a regular file: " + path);
    }
    std::filesystem::path temporary = target;
    temporary += ".tmp." + std::to_string(static_cast<long long>(::getpid())) +
                 "." + std::to_string(sequence.fetch_add(1));
    if (content.size() >
        static_cast<std::size_t>(
            std::numeric_limits<std::streamsize>::max())) {
        throw std::runtime_error(
            "run manifest exceeds stream size limits");
    }
    std::ofstream output(
        temporary, std::ios::binary | std::ios::out | std::ios::trunc);
    if (!output) {
        throw std::runtime_error(
            "cannot stage manifest file " + path + " via " +
            temporary.string());
    }
    output.write(content.data(), static_cast<std::streamsize>(content.size()));
    output.flush();
    if (!output) {
        output.close();
        std::error_code ignored;
        std::filesystem::remove(temporary, ignored);
        throw std::runtime_error(
            "failed while staging manifest file: " + path);
    }
    output.close();
    if (output.fail()) {
        std::error_code ignored;
        std::filesystem::remove(temporary, ignored);
        throw std::runtime_error(
            "failed to close staged manifest file: " + path);
    }
    return {target, temporary, "manifest"};
}

void discard_staged_configurations(
    const std::vector<StagedConfiguration>& staged) noexcept {
    for (const StagedConfiguration& output : staged) {
        std::error_code ignored;
        std::filesystem::remove(output.temporary, ignored);
    }
}

void publish_staged_configurations(
    std::vector<StagedConfiguration>& staged) {
    for (StagedConfiguration& output : staged) {
        std::error_code error;
        std::filesystem::rename(output.temporary, output.target, error);
        if (error) {
            throw std::runtime_error(
                "cannot publish " + output.artifact + " file " +
                output.target.string() +
                ": " + error.message());
        }
        output.temporary.clear();
    }
}

std::filesystem::path normalized_artifact_path(
    const std::string& path) {
    std::error_code error;
    std::filesystem::path absolute =
        std::filesystem::absolute(path, error);
    if (error) {
        throw std::runtime_error(
            "cannot resolve artifact path " + path + ": " +
            error.message());
    }
    std::filesystem::path normalized =
        std::filesystem::weakly_canonical(absolute, error);
    if (error) {
        normalized = absolute.lexically_normal();
    }
    return normalized;
}

std::filesystem::path resolved_executable_path(
    const char* invoked_path) {
    std::error_code error;
    const std::filesystem::path proc_target =
        std::filesystem::read_symlink("/proc/self/exe", error);
    if (!error && !proc_target.empty()) {
        return normalized_artifact_path(proc_target.string());
    }
    return normalized_artifact_path(invoked_path);
}

std::vector<std::pair<std::string, std::string>>
relevant_runtime_environment() {
    static constexpr std::array names{
        "DISKFORGE_VERIFY_REFERENCE",
        "OMP_NUM_THREADS",
        "OMP_DYNAMIC",
        "OMP_PROC_BIND",
        "OMP_PLACES",
        "OMP_WAIT_POLICY",
        "OMP_THREAD_LIMIT",
        "OMP_SCHEDULE",
        "GOMP_CPU_AFFINITY",
        "KMP_AFFINITY",
        "LD_LIBRARY_PATH",
        "LANG",
        "LC_ALL",
        "LC_NUMERIC",
        "TZ",
    };
    std::vector<std::pair<std::string, std::string>> values;
    values.reserve(names.size());
    for (const std::string_view name : names) {
        const std::string owned_name(name);
        if (const char* value = std::getenv(owned_name.c_str());
            value != nullptr) {
            values.emplace_back(owned_name, value);
        }
    }
    return values;
}

std::string floating_point_rounding_mode() {
    switch (std::fegetround()) {
        case FE_TONEAREST:
            return "to-nearest";
        case FE_DOWNWARD:
            return "downward";
        case FE_UPWARD:
            return "upward";
        case FE_TOWARDZERO:
            return "toward-zero";
        default:
            return "unknown";
    }
}

void require_distinct_artifact_paths(
    const std::vector<std::string>& paths) {
    std::vector<std::filesystem::path> normalized;
    normalized.reserve(paths.size());
    for (const std::string& path : paths) {
        const std::filesystem::path candidate =
            normalized_artifact_path(path);
        for (std::size_t previous = 0;
             previous < normalized.size(); ++previous) {
            std::error_code equivalent_error;
            const bool equivalent =
                std::filesystem::equivalent(
                    candidate, normalized[previous], equivalent_error);
            if (candidate == normalized[previous] ||
                (!equivalent_error && equivalent)) {
                throw std::invalid_argument(
                    "artifact paths alias each other: " +
                    paths[previous] + " and " + path);
            }
        }
        normalized.push_back(candidate);
    }
}

void print_initialization(const InitialConfiguration& initial, const Parameters& parameters) {
    std::cout << std::setprecision(10) << "INIT"
              << " initializer=" << initializer_name(parameters.initializer)
              << " n=" << initial.points.size() << " seed=" << parameters.seed
              << " box=" << initial.box << " radius=" << initial.radius
              << " pf=" << parameters.initial_packing_fraction
              << " seconds=" << initial.stats.seconds << " trials=" << initial.stats.trials
              << " accepted=" << initial.stats.accepted
              << " distance_checks=" << initial.stats.distance_checks
              << " rng_draws=" << initial.stats.rng_draws << '\n';
}

void print_result(const RunResult& result, const Parameters& parameters) {
    const double packing_fraction =
        static_cast<double>(result.points.size()) * kPi * result.radius * result.radius /
        (result.box * result.box);
    std::cout << std::setprecision(10) << "RESULT"
              << " engine=" << engine_name(result.engine) << " n=" << result.points.size()
              << " seed=" << parameters.seed << " threads=" << result.threads_used
              << " radius=" << result.radius << " pf=" << packing_fraction
              << " seconds=" << result.seconds << " oracle_seconds=" << result.oracle.seconds
              << " engine_setup_seconds=" << result.setup_seconds
              << " discovery_seconds=" << result.discovery_seconds
              << " relaxation_seconds=" << result.relaxation_seconds
              << " scheduler_seconds=" << result.scheduler_seconds
              << " checksum=0x" << std::hex << result.state_checksum << std::dec
              << " oracle=" << (result.oracle.performed ? "performed" : "skipped");
    if (result.oracle.performed) {
        std::cout << " oracle_pairs=" << result.oracle.pairs_checked
                  << " oracle_overlaps=" << result.oracle.overlapping_pairs
                  << " min_distance=" << result.oracle.minimum_distance;
    }
    std::cout << '\n';
    const Counters& counters = result.counters;
    std::cout << "COUNTERS"
              << " engine=" << engine_name(result.engine)
              << " steps=" << counters.compression_steps
              << " overlap_particle_scans=" << counters.overlap_particle_scans
              << " representatives=" << counters.representatives
              << " neighbor_lists=" << counters.neighbor_lists
              << " neighbor_visits=" << counters.neighbor_visits
              << " move_trials=" << counters.move_trials
              << " accepted_moves=" << counters.accepted_moves
              << " distance_checks=" << counters.distance_checks
              << " discovery_distance_checks=" << counters.discovery_distance_checks
              << " relaxation_distance_checks=" << counters.relaxation_distance_checks
              << " scheduler_distance_checks=" << counters.scheduler_distance_checks
              << " candidate_point_queries=" << counters.candidate_point_queries
              << " event_candidate_pairs=" << counters.event_candidate_pairs
              << " event_events_scheduled=" << counters.event_events_scheduled
              << " event_initial_events=" << counters.event_initial_events
              << " event_local_events=" << counters.event_local_events
              << " event_rebuild_events=" << counters.event_rebuild_events
              << " event_valid_events=" << counters.event_valid_events
              << " event_stale_events=" << counters.event_stale_events
              << " event_generation_bumps=" << counters.event_generation_bumps
              << " event_dirty_particles=" << counters.event_dirty_particles
              << " event_local_refreshes=" << counters.event_local_refreshes
              << " event_rebuilds=" << counters.event_rebuilds
              << " event_rebuild_dropped_events="
              << counters.event_rebuild_dropped_events
              << " event_peak_queued=" << counters.event_peak_queued
              << " event_queued_end=" << counters.event_queued_end
              << " sqrt_calls=" << counters.sqrt_calls
              << " temporary_vectors=" << counters.temporary_vectors
              << " cell_removals=" << counters.cell_removals
              << " cell_additions=" << counters.cell_additions
              << " radius_writes=" << counters.radius_writes
              << " anchor_writes=" << counters.anchor_writes
              << " rng_draws=" << counters.rng_draws << '\n';
}

std::string render_json_result(const RunResult& result,
                               const InitialConfiguration& initial,
                               const OracleResult& initial_oracle,
                               const Parameters& parameters,
                               const RunResult* baseline) {
    const bool validation_performed = initial_oracle.performed && result.oracle.performed;
    const double validation_seconds =
        validation_performed ? initial_oracle.seconds + result.oracle.seconds : 0.0;
    const double production_seconds = initial.stats.seconds + result.seconds;
    const double total_seconds = production_seconds + validation_seconds;
    const double packing_fraction =
        static_cast<double>(result.points.size()) * kPi * result.radius * result.radius /
        (result.box * result.box);
    const Counters& counters = result.counters;
    std::ostringstream output;

    output << std::setprecision(17) << '{'
              << "\"engine\":\"" << engine_name(result.engine) << "\","
              << "\"initializer\":\"" << initializer_name(parameters.initializer) << "\","
              << "\"n\":" << result.points.size() << ','
              << "\"seed\":" << parameters.seed << ','
              << "\"threads\":" << result.threads_used << ','
              << "\"timing_seconds\":{"
              << "\"initialization\":" << initial.stats.seconds << ','
              << "\"compression\":" << result.seconds << ','
              << "\"engine_setup\":" << result.setup_seconds << ','
              << "\"discovery\":" << result.discovery_seconds << ','
              << "\"relaxation\":" << result.relaxation_seconds << ','
              << "\"scheduler\":" << result.scheduler_seconds << ','
              << "\"production\":" << production_seconds << ','
              << "\"validation\":" << validation_seconds << ','
              << "\"total\":" << total_seconds << "},"
              << "\"state\":{"
              << "\"box\":" << result.box << ','
              << "\"radius\":" << result.radius << ','
              << "\"packing_fraction\":" << packing_fraction << ','
              << "\"checksum\":\"0x" << std::hex << result.state_checksum << std::dec << "\"},"
              << "\"valid\":";
    if (validation_performed) {
        output << (result.oracle.overlapping_pairs == 0 ? "true" : "false");
    } else {
        output << "null";
    }
    output << ",\"oracle\":";
    if (validation_performed) {
        output << '{'
                  << "\"pairs\":" << result.oracle.pairs_checked << ','
                  << "\"overlaps\":" << result.oracle.overlapping_pairs << ','
                  << "\"minimum_distance\":";
        if (std::isfinite(result.oracle.minimum_distance)) {
            output << result.oracle.minimum_distance;
        } else {
            output << "null";
        }
        output << '}';
    } else {
        output << "null";
    }
    output << ",\"counters\":{"
              << "\"steps\":" << counters.compression_steps << ','
              << "\"overlap_particle_scans\":" << counters.overlap_particle_scans << ','
              << "\"representatives\":" << counters.representatives << ','
              << "\"neighbor_lists\":" << counters.neighbor_lists << ','
              << "\"neighbor_visits\":" << counters.neighbor_visits << ','
              << "\"move_trials\":" << counters.move_trials << ','
              << "\"accepted_moves\":" << counters.accepted_moves << ','
              << "\"distance_checks\":" << counters.distance_checks << ','
              << "\"discovery_distance_checks\":" << counters.discovery_distance_checks << ','
              << "\"relaxation_distance_checks\":" << counters.relaxation_distance_checks << ','
              << "\"scheduler_distance_checks\":" << counters.scheduler_distance_checks << ','
              << "\"candidate_point_queries\":" << counters.candidate_point_queries << ','
              << "\"event_candidate_pairs\":" << counters.event_candidate_pairs << ','
              << "\"event_events_scheduled\":" << counters.event_events_scheduled << ','
              << "\"event_initial_events\":" << counters.event_initial_events << ','
              << "\"event_local_events\":" << counters.event_local_events << ','
              << "\"event_rebuild_events\":" << counters.event_rebuild_events << ','
              << "\"event_valid_events\":" << counters.event_valid_events << ','
              << "\"event_stale_events\":" << counters.event_stale_events << ','
              << "\"event_generation_bumps\":" << counters.event_generation_bumps << ','
              << "\"event_dirty_particles\":" << counters.event_dirty_particles << ','
              << "\"event_local_refreshes\":" << counters.event_local_refreshes << ','
              << "\"event_rebuilds\":" << counters.event_rebuilds << ','
              << "\"event_rebuild_dropped_events\":"
              << counters.event_rebuild_dropped_events << ','
              << "\"event_peak_queued\":" << counters.event_peak_queued << ','
              << "\"event_queued_end\":" << counters.event_queued_end << ','
              << "\"sqrt_calls\":" << counters.sqrt_calls << ','
              << "\"temporary_vectors\":" << counters.temporary_vectors << ','
              << "\"cell_removals\":" << counters.cell_removals << ','
              << "\"cell_additions\":" << counters.cell_additions << ','
              << "\"radius_writes\":" << counters.radius_writes << ','
              << "\"anchor_writes\":" << counters.anchor_writes << ','
              << "\"rng_draws\":" << counters.rng_draws << '}';
    if (baseline != nullptr) {
        output << ",\"comparison\":{"
                  << "\"baseline\":\"" << engine_name(baseline->engine) << "\","
                  << "\"speedup\":" << baseline->seconds / result.seconds << ','
                  << "\"bitwise_equal\":"
                  << (positions_bitwise_equal(*baseline, result) ? "true" : "false") << ','
                  << "\"max_coordinate_difference\":"
                  << maximum_coordinate_difference(*baseline, result) << '}';
    }
    output << '}';
    if (!output) {
        throw std::runtime_error("failed to render JSON result");
    }
    return output.str();
}

void require(bool condition, std::string_view message) {
    if (!condition) {
        throw std::runtime_error("self-test failure: " + std::string(message));
    }
}

bool relaxation_counters_equal(const Counters& lhs, const Counters& rhs) {
    return lhs.compression_steps == rhs.compression_steps &&
           lhs.representatives == rhs.representatives &&
           lhs.neighbor_lists == rhs.neighbor_lists &&
           lhs.neighbor_visits == rhs.neighbor_visits &&
           lhs.move_trials == rhs.move_trials &&
           lhs.accepted_moves == rhs.accepted_moves &&
           lhs.relaxation_distance_checks == rhs.relaxation_distance_checks &&
           lhs.candidate_point_queries == rhs.candidate_point_queries &&
           lhs.cell_removals == rhs.cell_removals &&
           lhs.cell_additions == rhs.cell_additions &&
           lhs.anchor_writes == rhs.anchor_writes &&
           lhs.rng_draws == rhs.rng_draws;
}

int self_test() {
    const auto reference_wrap = [](double value, double box) {
        value -= box * std::floor(value / box);
        if (value >= box) {
            value = 0.0;
        }
        return value;
    };
    constexpr double wrap_box = 7.25;
    std::vector<double> wrap_values{
        0.0,
        -0.0,
        std::nextafter(0.0, -std::numeric_limits<double>::infinity()),
        std::nextafter(0.0, std::numeric_limits<double>::infinity()),
        wrap_box,
        -wrap_box,
        std::nextafter(wrap_box, 0.0),
        std::nextafter(wrap_box, std::numeric_limits<double>::infinity()),
        std::nextafter(-wrap_box, -std::numeric_limits<double>::infinity()),
        std::nextafter(-wrap_box, 0.0),
        1234.56789,
        -1234.56789,
    };
    for (int multiple = -12; multiple <= 12; ++multiple) {
        wrap_values.push_back(static_cast<double>(multiple) * wrap_box);
    }
    for (const double value : wrap_values) {
        require(std::bit_cast<std::uint64_t>(wrap_coordinate(value, wrap_box)) ==
                    std::bit_cast<std::uint64_t>(reference_wrap(value, wrap_box)),
                "wrap fast path differs from the floor reference");
    }

    Parameters parameters;
    parameters.particle_count = 36;
    parameters.seed = 424242;
    parameters.initial_packing_fraction = 0.30;
    parameters.final_packing_fraction = 0.36;
    parameters.final_radius = 0.5;
    parameters.compression_a = 100.0;
    parameters.epsilon_floor = 0.004;
    parameters.mc_sweeps = 3;
    parameters.mc_step_scale = 0.75;
    parameters.max_trials_per_move = 100000;
    parameters.threads = 4;
    parameters.oracle_threads = 1;
    validate_parameters(parameters);

    const InitialConfiguration initial = initialize(parameters);
    const OracleResult initial_oracle =
        final_oracle(initial.points, initial.radius, initial.box);
    require(initial_oracle.overlapping_pairs == 0, "lattice-MC initializer has overlaps");

#ifdef _OPENMP
    const int omp_max_threads_before = omp_get_max_threads();
    const int omp_dynamic_before = omp_get_dynamic();
    const int omp_thread_limit = omp_get_thread_limit();
#endif
    const RunControls verified_steps{StepVerification::discovery_and_clear, false};
    const RunControls verified_rebuilds{StepVerification::discovery_and_clear, true};
    const RunResult legacy =
        run_engine(initial, parameters, Engine::scan_legacy, verified_steps);
    const RunResult fast =
        run_engine(initial, parameters, Engine::scan_fast, verified_steps);
    const RunResult omp =
        run_engine(initial, parameters, Engine::scan_omp, verified_steps);
    const RunResult event =
        run_engine(initial, parameters, Engine::event, verified_steps);
    const RunResult event_rebuild =
        run_engine(initial, parameters, Engine::event, verified_rebuilds);
    const RunResult fast_repeat = run_engine(initial, parameters, Engine::scan_fast);
    std::vector<RunResult> omp_matrix;
    for (const int thread_count : std::vector<int>{1, 2, 3, 0, 3}) {
        Parameters omp_parameters = parameters;
        omp_parameters.threads = thread_count;
        omp_matrix.push_back(
            run_engine(initial, omp_parameters, Engine::scan_omp, verified_steps));
    }

    const auto require_complete_step_verification =
        [&](const RunResult& result, std::string_view label) {
            const std::uint64_t pairs_per_scan =
                static_cast<std::uint64_t>(result.points.size()) *
                static_cast<std::uint64_t>(result.points.size() - 1U) / 2U;
            require(result.verification.discovery_scans ==
                        result.counters.compression_steps &&
                        result.verification.post_relaxation_scans ==
                            result.counters.compression_steps,
                    std::string(label) +
                        " did not run both independent scans at every compression step");
            require(result.verification.pairs_checked ==
                        2U * result.counters.compression_steps * pairs_per_scan,
                    std::string(label) +
                        " independent verification pair accounting differs");
        };
    require_complete_step_verification(legacy, "legacy");
    require_complete_step_verification(fast, "fast");
    require_complete_step_verification(omp, "OpenMP");
    require_complete_step_verification(event, "event");
    require_complete_step_verification(event_rebuild, "rebuilding event");
    for (const RunResult& matrix_result : omp_matrix) {
        require_complete_step_verification(matrix_result, "OpenMP matrix");
    }
    require(legacy.oracle.overlapping_pairs == 0, "legacy final oracle found an overlap");
    require(fast.oracle.overlapping_pairs == 0, "fast final oracle found an overlap");
    require(omp.oracle.overlapping_pairs == 0, "OpenMP final oracle found an overlap");
    require(event.oracle.overlapping_pairs == 0, "event final oracle found an overlap");
    require(event_rebuild.oracle.overlapping_pairs == 0,
            "rebuilding event final oracle found an overlap");
    require(positions_bitwise_equal(legacy, fast),
            "legacy and fast positions are not bitwise identical");
    require(positions_bitwise_equal(fast, omp),
            "fast and OpenMP positions are not bitwise identical");
    for (const RunResult& matrix_result : omp_matrix) {
        require(positions_bitwise_equal(fast, matrix_result),
                "OpenMP thread matrix changed the trajectory");
        require(matrix_result.counters == fast.counters,
                "OpenMP thread matrix changed deterministic counters");
        require(matrix_result.oracle.pairs_checked == fast.oracle.pairs_checked &&
                    matrix_result.oracle.overlapping_pairs ==
                        fast.oracle.overlapping_pairs &&
                    matrix_result.oracle.minimum_distance ==
                        fast.oracle.minimum_distance,
                "OpenMP thread matrix changed oracle results");
    }
    require(positions_bitwise_equal(omp_matrix[2], omp_matrix[4]) &&
                omp_matrix[2].counters == omp_matrix[4].counters,
            "repeated three-thread OpenMP run is not deterministic");
#ifdef _OPENMP
    const std::vector<int> requested_thread_counts{1, 2, 3, 0, 3};
    for (std::size_t index = 0; index < omp_matrix.size(); ++index) {
        const int requested = requested_thread_counts[index];
        const int configured =
            requested > 0 ? requested : omp_max_threads_before;
        require(omp_matrix[index].threads_used ==
                    std::min(configured, omp_thread_limit),
                "OpenMP run did not use the requested bounded team size: requested=" +
                    std::to_string(requested) +
                    " configured=" + std::to_string(configured) +
                    " limit=" + std::to_string(omp_thread_limit) +
                    " actual=" +
                    std::to_string(omp_matrix[index].threads_used));
    }
    require(omp_get_max_threads() == omp_max_threads_before &&
                omp_get_dynamic() == omp_dynamic_before,
            "OpenMP discovery leaked runtime settings into the caller");
#endif
    require(positions_bitwise_equal(fast, event),
            "fast and event positions are not bitwise identical");
    require(positions_bitwise_equal(fast, event_rebuild),
            "fast and rebuilding-event positions are not bitwise identical");
    RunResult altered_box = fast;
    altered_box.box =
        std::nextafter(altered_box.box, std::numeric_limits<double>::infinity());
    require(!positions_bitwise_equal(fast, altered_box),
            "bitwise state comparison ignored the periodic box");
    require(relaxation_counters_equal(fast.counters, event.counters),
            "fast and event RNG-sensitive counters differ");
    require(relaxation_counters_equal(fast.counters, event_rebuild.counters),
            "fast and rebuilding-event RNG-sensitive counters differ");
    require(fast.state_checksum == fast_repeat.state_checksum,
            "repeated fast run is not deterministic");
    require(legacy.counters.sqrt_calls > 0 && legacy.counters.temporary_vectors > 0,
            "legacy path did not exercise sqrt/allocation baseline");
    require(fast.counters.sqrt_calls >= fast.counters.compression_steps &&
                fast.counters.sqrt_calls <=
                    2U * fast.counters.compression_steps &&
                fast.counters.sqrt_calls < legacy.counters.sqrt_calls &&
                fast.counters.temporary_vectors == 0,
            "fast path did not cache one exact contact cutoff per step");
    require(event.counters.sqrt_calls >=
                    event.counters.scheduler_distance_checks +
                        event.counters.compression_steps &&
                event.counters.sqrt_calls <=
                    event.counters.scheduler_distance_checks +
                        2U * event.counters.compression_steps,
            "event path did not separate scheduler sqrt calls from cached cutoffs");
    require(legacy.counters.anchor_writes ==
                legacy.counters.compression_steps * parameters.particle_count,
            "legacy did not refresh every anchor");
    require(fast.counters.anchor_writes <= legacy.counters.anchor_writes,
            "dirty-anchor optimization wrote more anchors than legacy");
    require(legacy.counters.cell_removals == legacy.counters.move_trials &&
                legacy.counters.cell_additions == legacy.counters.move_trials,
            "legacy path no longer mutates the cell list on every trial");
    require(fast.counters.cell_removals == fast.counters.accepted_moves &&
                fast.counters.cell_additions == fast.counters.accepted_moves,
            "fast path did not reduce cell mutations to accepted moves");
    require(fast.counters.candidate_point_queries == fast.counters.move_trials,
            "fast path did not query every candidate point directly");
    require(legacy.counters.rng_draws ==
                    2U * legacy.counters.move_trials &&
                fast.counters.rng_draws ==
                    2U * fast.counters.move_trials &&
                event.counters.rng_draws ==
                    2U * event.counters.move_trials,
            "paired relocation draws diverged from the trial count");
    require(fast.counters.distance_checks ==
                fast.counters.discovery_distance_checks +
                    fast.counters.relaxation_distance_checks +
                    fast.counters.scheduler_distance_checks,
            "fast distance-check phase counters do not sum to the aggregate");
    require(event.counters.distance_checks ==
                event.counters.discovery_distance_checks +
                    event.counters.relaxation_distance_checks +
                    event.counters.scheduler_distance_checks,
            "event distance-check phase counters do not sum to the aggregate");
    require(event.counters.event_events_scheduled > 0 &&
                event.counters.event_queued_end == 0,
            "event path did not exercise and drain its calendar");
    require(event.counters.event_generation_bumps == event.counters.anchor_writes,
            "event generations and dirty-anchor updates diverged");
    require(event_rebuild.counters.event_rebuilds > 0 &&
                event_rebuild.counters.event_rebuild_events > 0,
            "forced event rebuild path was not exercised");

    const std::vector<double> strict_thresholds{1.0, 4.0, 9.0};
    require(first_strict_event_step(strict_thresholds, 4.0, 0) == 2,
            "exact tangency was scheduled as an overlap");
    require(first_strict_event_step(
                strict_thresholds,
                std::nextafter(4.0, -std::numeric_limits<double>::infinity()), 0) == 1,
            "value immediately below tangency was not scheduled at the contact step");
    require(first_strict_event_step(
                strict_thresholds,
                std::nextafter(4.0, std::numeric_limits<double>::infinity()), 0) == 2,
            "value immediately above tangency was scheduled too early");
    require(first_strict_event_step(strict_thresholds, 9.0, 0) ==
                strict_thresholds.size(),
            "final-step tangency incorrectly produced an event");

    const double contact_radius = 0x1.d391075c3a5f8p-2;
    const double contact_dx = 0x1.d37ac47146d6cp-1;
    const double contact_dy = 0x1.208d140f826dap-6;
    const double contact_d2 = contact_dx * contact_dx + contact_dy * contact_dy;
    const double contact_diameter = contact_radius + contact_radius;
    require(contact_d2 < contact_diameter * contact_diameter &&
                std::sqrt(contact_d2) == contact_diameter,
            "contact regression no longer exercises squared-predicate disagreement");
    std::uint64_t contact_sqrt_calls = 0;
    require(!strict_overlap_from_squared(contact_d2, contact_diameter,
                                         &contact_sqrt_calls) &&
                contact_sqrt_calls == 1,
            "guarded contact predicate differs from the legacy strict predicate");

    const double subnormal_diameter = 0x1.0000000000001p-537;
    const double subnormal_d2 = subnormal_diameter * subnormal_diameter;
    require(strict_overlap_from_squared(subnormal_d2, subnormal_diameter),
            "guarded contact predicate mishandled subnormal product rounding");
    require(!strict_overlap_from_squared(
                std::nextafter(contact_diameter * contact_diameter,
                               std::numeric_limits<double>::infinity()),
                contact_diameter),
            "guarded contact predicate accepted a distance above the squared guard");
    const double overflow_diameter =
        std::nextafter(std::sqrt(std::numeric_limits<double>::max()),
                       std::numeric_limits<double>::infinity());
    require(std::isinf(overflow_diameter * overflow_diameter) &&
                strict_overlap_from_squared(std::numeric_limits<double>::max(),
                                            overflow_diameter),
            "guarded contact predicate mishandled squared-threshold overflow");

    bool non_finite_oracle_rejected = false;
    try {
        const std::vector<Point> non_finite_points{
            {std::numeric_limits<double>::quiet_NaN(), 0.0}};
        (void)final_oracle(non_finite_points, 0.5, 2.0);
    } catch (const std::invalid_argument&) {
        non_finite_oracle_rejected = true;
    }
    require(non_finite_oracle_rejected, "oracle accepted a non-finite coordinate");

    PeriodicCellList incomplete_cells(4.0, 1.0, 2);
    const std::vector<Point> incomplete_points{{0.5, 0.5}, {2.5, 2.5}};
    incomplete_cells.add(0, incomplete_points[0]);
    require(!incomplete_cells.validate(incomplete_points),
            "cell-list validation accepted a missing particle");
    bool invalid_cell_coordinate_rejected = false;
    try {
        incomplete_cells.add(
            1, {std::numeric_limits<double>::quiet_NaN(), 2.5});
    } catch (const std::invalid_argument&) {
        invalid_cell_coordinate_rejected = true;
    }
    require(invalid_cell_coordinate_rejected,
            "cell list accepted a non-finite coordinate");

    constexpr double rounded_box = 1.6;
    constexpr double rounded_cutoff = 0.1;
    PeriodicCellList rounded_cells(rounded_box, rounded_cutoff, 16);
    std::vector<Point> rounded_points;
    rounded_points.reserve(16);
    for (int y = 12; y <= 15; ++y) {
        for (int x = 12; x <= 15; ++x) {
            rounded_points.push_back(
                {(static_cast<double>(x) + 0.5) * rounded_cutoff,
                 (static_cast<double>(y) + 0.5) * rounded_cutoff});
            rounded_cells.add(
                static_cast<int>(rounded_points.size() - 1),
                rounded_points.back());
        }
    }
    Workspace rounded_workspace(rounded_points.size());
    const double rounded_anchor =
        std::nextafter(14.0 * rounded_cutoff, 0.0);
    const std::array<double, 4> rounded_displacements{
        -rounded_cutoff,
        0.0,
        0.5 * rounded_cutoff,
        (1.0 - 0x1.0p-52) * rounded_cutoff,
    };
    std::vector<int> rounded_candidate_cells;
    for (const double dx : rounded_displacements) {
        for (const double dy : rounded_displacements) {
            const int cell = rounded_cells.cell_for_point(
                {rounded_anchor + dx, rounded_anchor + dy});
            if (std::find(
                    rounded_candidate_cells.begin(),
                    rounded_candidate_cells.end(), cell) !=
                rounded_candidate_cells.end()) {
                continue;
            }
            rounded_candidate_cells.push_back(cell);
            const std::vector<int>& cached =
                rounded_workspace.candidate_neighborhood(
                    rounded_cells, cell);
            std::vector<int> direct;
            rounded_cells.gather_cell_neighbors(cell, direct);
            require(cached == direct,
                    "overflow candidate-neighborhood cache changed order");
        }
    }
    require(rounded_candidate_cells.size() == 16,
            "cache overflow self-test did not visit 16 rounded cells");

    Parameters periodic_parameters = parameters;
    periodic_parameters.particle_count = 2;
    periodic_parameters.compression_a = 100.0;
    periodic_parameters.epsilon_floor = 0.002;
    periodic_parameters.max_compression_steps = 1000;
    ReproducibleRng periodic_rng(991);
    InitStats periodic_stats;
    const InitialConfiguration edge_initial(
        {{0.10, 5.0}, {9.15, 5.0}}, periodic_rng, 10.0, 0.45, periodic_stats);
    const RunResult edge_fast =
        run_engine(edge_initial, periodic_parameters, Engine::scan_fast,
                   verified_steps);
    const RunResult edge_event =
        run_engine(edge_initial, periodic_parameters, Engine::event,
                   verified_rebuilds);
    require(positions_bitwise_equal(edge_fast, edge_event),
            "periodic-edge event trajectory differs from scan-fast");

    ReproducibleRng corner_rng(992);
    const InitialConfiguration corner_initial(
        {{0.10, 0.10}, {9.40, 9.40}}, corner_rng, 10.0, 0.45, periodic_stats);
    const RunResult corner_fast =
        run_engine(corner_initial, periodic_parameters, Engine::scan_fast,
                   verified_steps);
    const RunResult corner_event =
        run_engine(corner_initial, periodic_parameters, Engine::event,
                   verified_rebuilds);
    require(positions_bitwise_equal(corner_fast, corner_event),
            "periodic-corner event trajectory differs from scan-fast");

    Parameters rsa_parameters = parameters;
    rsa_parameters.particle_count = 16;
    rsa_parameters.seed = 17;
    rsa_parameters.initial_packing_fraction = 0.10;
    rsa_parameters.final_packing_fraction = 0.10;
    rsa_parameters.initializer = Initializer::rsa;
    rsa_parameters.rsa_attempts_per_particle = 10000;
    const InitialConfiguration rsa = initialize(rsa_parameters);
    require(final_oracle(rsa.points, rsa.radius, rsa.box).overlapping_pairs == 0,
            "RSA initializer has overlaps");

    std::cout << "SELF_TEST PASS"
              << " checksum=0x" << std::hex << fast.state_checksum << std::dec
              << " steps=" << fast.counters.compression_steps
              << " legacy_seconds=" << legacy.seconds << " fast_seconds=" << fast.seconds
              << " omp_seconds=" << omp.seconds << " event_seconds=" << event.seconds
              << " event_rebuilds=" << event_rebuild.counters.event_rebuilds << '\n';
    return 0;
}

int run_command(const Parameters& parameters,
                const std::string& manifest_path,
                int argc,
                char** argv) {
    const InitialConfiguration initial = initialize(parameters);
    OracleResult initial_oracle;
    if (!parameters.skip_oracle) {
        initial_oracle = final_oracle(
            initial.points, initial.radius, initial.box,
            parameters.oracle_threads);
    }
    if (initial_oracle.performed && initial_oracle.overlapping_pairs != 0) {
        throw std::runtime_error("initializer failed the independent O(N^2) oracle");
    }

    std::vector<RunResult> results;
    if (parameters.engine == "all") {
        results.push_back(run_engine(initial, parameters, Engine::scan_legacy));
        results.push_back(run_engine(initial, parameters, Engine::scan_fast));
        results.push_back(run_engine(initial, parameters, Engine::scan_omp));
        results.push_back(run_engine(initial, parameters, Engine::event));
    } else {
        results.push_back(run_engine(initial, parameters, parse_engine(parameters.engine)));
    }

    bool oracle_failure = false;
    for (const RunResult& result : results) {
        oracle_failure = oracle_failure ||
                         (result.oracle.performed &&
                          result.oracle.overlapping_pairs != 0);
    }

    std::vector<std::string> rendered_results;
    if (parameters.json || !manifest_path.empty()) {
        rendered_results.reserve(results.size());
        for (std::size_t index = 0; index < results.size(); ++index) {
            const RunResult* baseline =
                results.size() > 1 && index > 0
                    ? &results.front()
                    : nullptr;
            rendered_results.push_back(render_json_result(
                results[index], initial, initial_oracle,
                parameters, baseline));
        }
    }

    std::vector<std::string> coordinate_paths;
    if (!parameters.output_path.empty()) {
        coordinate_paths.reserve(results.size());
        for (const RunResult& result : results) {
            std::string path = parameters.output_path;
            if (results.size() > 1) {
                path += "." + std::string(engine_name(result.engine));
            }
            coordinate_paths.push_back(std::move(path));
        }
    }

    std::vector<StagedConfiguration> staged;
    if (!oracle_failure &&
        (!coordinate_paths.empty() ||
         !manifest_path.empty())) {
        try {
            std::vector<std::string> artifact_paths =
                coordinate_paths;
            if (!manifest_path.empty()) {
                artifact_paths.push_back(manifest_path);
            }
            require_distinct_artifact_paths(artifact_paths);

            staged.reserve(artifact_paths.size());
            for (std::size_t index = 0;
                 index < coordinate_paths.size(); ++index) {
                staged.push_back(stage_configuration(
                    results[index], parameters,
                    coordinate_paths[index]));
            }
            if (!manifest_path.empty()) {
                RunManifestData manifest;
                manifest.argv.reserve(static_cast<std::size_t>(argc));
                for (int index = 0; index < argc; ++index) {
                    manifest.argv.emplace_back(argv[index]);
                }
                manifest.working_directory =
                    std::filesystem::current_path().string();
                manifest.executable =
                    resolved_executable_path(argv[0]).string();
                manifest.manifest_path = manifest_path;
                manifest.floating_point_rounding_mode =
                    floating_point_rounding_mode();
                manifest.environment =
                    relevant_runtime_environment();
                manifest.parameters = parameters;
                manifest.initial_particle_count =
                    initial.points.size();
                manifest.initial_box = initial.box;
                manifest.initial_radius = initial.radius;
                manifest.initial_checksum = checksum(
                    initial.points, initial.radius, initial.box);
                ReproducibleRng continuation =
                    initial.rng_after_initialization;
                manifest.next_rng_raw = continuation.raw();
                manifest.initialization = initial.stats;
                manifest.initial_oracle = initial_oracle;
                manifest.result_json_objects = rendered_results;
                manifest.coordinate_paths = coordinate_paths;
                staged.push_back(stage_manifest(
                    render_run_manifest(manifest),
                    manifest_path));
            }
            publish_staged_configurations(staged);
        } catch (...) {
            discard_staged_configurations(staged);
            throw;
        }
    }

    if (!parameters.json) {
        print_initialization(initial, parameters);
    }
    for (std::size_t index = 0; index < results.size(); ++index) {
        const RunResult& result = results[index];
        if (parameters.json) {
            std::cout << rendered_results[index] << '\n';
        } else {
            print_result(result, parameters);
        }
    }

    if (results.size() > 1 && !parameters.json) {
        const RunResult& legacy = results[0];
        for (std::size_t index = 1; index < results.size(); ++index) {
            const RunResult& candidate = results[index];
            std::cout << std::setprecision(10) << "COMPARISON"
                      << " baseline=" << engine_name(legacy.engine)
                      << " candidate=" << engine_name(candidate.engine)
                      << " speedup=" << legacy.seconds / candidate.seconds
                      << " bitwise_equal="
                      << (positions_bitwise_equal(legacy, candidate) ? "yes" : "no")
                      << " max_coordinate_difference="
                      << maximum_coordinate_difference(legacy, candidate) << '\n';
        }
    }

    if (oracle_failure) {
        std::cerr << "ERROR: final O(N^2) oracle found overlapping pairs\n";
        return 3;
    }
    return 0;
}

int app_main(int argc, char** argv) {
    try {
        if (argc >= 2) {
            const std::string_view command(argv[1]);
            if (command == "self-test" || command == "--self-test") {
                return self_test();
            }
            if (command == "--help" || command == "-h" || command == "help") {
                print_usage(std::cout);
                return 0;
            }
            if (command == "run") {
                const ParsedRunArguments parsed =
                    parse_run_arguments(argc, argv, 2);
                return run_command(
                    parsed.parameters, parsed.manifest_path, argc, argv);
            }
        }
        const ParsedRunArguments parsed =
            parse_run_arguments(argc, argv, 1);
        return run_command(
            parsed.parameters, parsed.manifest_path, argc, argv);
    } catch (const std::exception& error) {
        std::cerr << "ERROR: " << error.what() << '\n';
        return 2;
    }
}

}  // namespace diskforge