//#define USE_CAIRO
#include "util.h"
#include "solver.h"
#include "wide_real.h"
#include "GlobalCheckpoint.h"
#include "SearchArchive.h"
#include "cuda_search.h"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <Eigen/Dense>
#include <filesystem>
#include <algorithm>
#include <array>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <limits>
#include <random>
#include <atomic>
#include <chrono>
#include <cstdint>
#include <mutex>
#include <memory>
#include <map>
#include <numeric>
#include <thread>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#ifdef _WIN32
#ifndef NOMINMAX
#define NOMINMAX
#endif
#include <windows.h>
#endif
#ifdef USE_CAIRO
#include <cairo.h>
#endif

#define NUM_TOPOLOGIES 59
#define DUAL_PROBLEM 0
constexpr int GLOBAL_PLANE_VALUE_COUNT = 36;
constexpr std::uint32_t GLOBAL_OBJECTIVE_VERSION = 4;
constexpr std::uint32_t MIN_COMPATIBLE_ARCHIVE_OBJECTIVE_VERSION = 3;
constexpr int CUDA_SESSION_CACHE_LIMIT = NUM_TOPOLOGIES;

struct StudyOptions {
    std::string obj_path = "data/shape_c2_i0_0.obj";
    std::string out_prefix = "runtime/candidates/top4_study";
    std::string report_path = "runtime/reports/02_study_cpp.md";
    std::string objective = "cross-zint";
    int topology = 4;
    int seed = 30000157;
    int max_iters = 360;
    int clusters = 10000;
    double sigma = 1e-2;
    double beta = 0.9;
    bool use_symmetry = false;
};

struct LocalRepairOptions {
    std::string obj_path = "data/shape_c2_i0_0.obj";
    std::string start_planes_path;
    std::string stop_file_path;
    std::string global_dir = "results/search";
    std::string out_prefix = "runtime/candidates/top4_local";
    std::string report_path = "runtime/reports/03_repair_local_cpp.md";
    int topology = 4;
    int seed = 30000157;
    int iterations = 50000;
    int report_every = 1000;
    double step = 0.5;
    double beta = 0.9995;
    double temperature = 0.02;
    std::vector<int> movable_vertices = {32, 33, 27, 4, 12, 23};
    int restarts = 256;
    int stagnation = 2000;
    int trials = 2000;
    int threads = 0;
    int time_limit_seconds = 0;
    double jump_chance = 0.08;
    double min_step_ratio = 1e-5;
    int topology_from = 0;
    int topology_to = NUM_TOPOLOGIES - 1;
    int cuda_chains = 0;
    int cuda_iterations = 64;
    int checkpoint_seconds = 30;
    double degeneracy_weight = 0.01;
    bool prioritize_worst = false;
    bool use_cuda = false;
};

bool save_plane_state(const std::filesystem::path& path, const VectorXd& x);
bool load_plane_state(const std::filesystem::path& path, VectorXd& x);

double (*select_objective(const std::string& name))(const VectorXd&) {
    if (name == "sum") { return objective_sum; }
    if (name == "sum-q") { return objective_sum_q; }
    if (name == "cross") { return objective_cross; }
    if (name == "cross-int") { return objective_cross_int; }
    if (name == "cross-int-q") { return objective_cross_int_q; }
    if (name == "cross-zint") { return objective_cross_zint; }
    if (name == "int-cross") { return objective_int_cross; }
    if (name == "int-zcross") { return objective_int_zcross; }
    if (name == "wsum") { return objective_wsum; }
    if (name == "wsum-q") { return objective_wsum_q; }
    return nullptr;
}

void print_usage(const char* exe_name) {
    std::cout
        << "Usage:\n"
        << "  " << exe_name << "\n"
        << "      Interactive random solver. Prompts for Seed and Topology.\n\n"
        << "  " << exe_name << " --study <obj> [options]\n"
        << "      Continue optimization from an existing OBJ candidate.\n\n"
        << "  " << exe_name << " --repair-local <obj> [options]\n"
        << "      Move local vertices [33,34,28,5,13,24] with a smooth crossing surrogate.\n\n"
        << "  " << exe_name << " --hunt-local <obj> [options]\n"
        << "      Fast focused search around the two known crossings on edge 33-34.\n\n"
        << "  " << exe_name << " --batch-hunt <obj> [options]\n"
        << "      Run many focused hunts with varied seeds and parameters.\n\n"
        << "  " << exe_name << " --global-search [options]\n"
        << "      Independent breadth/depth search over all 59 topologies.\n\n"
        << "Study options:\n"
        << "  --topology <n>       Default: 4\n"
        << "  --seed <n>           Default: 30000157\n"
        << "  --iters <n>          Default: 360\n"
        << "  --clusters <n>       Default: 10000\n"
        << "  --sigma <x>          Default: 0.01\n"
        << "  --beta <x>           Default: 0.9\n"
        << "  --objective <name>   Default: cross-zint\n"
        << "                       Names: sum, sum-q, cross, cross-int, cross-int-q,\n"
        << "                              cross-zint, int-cross, int-zcross, wsum, wsum-q\n"
        << "  --symmetry           Apply original O4 symmetry transform during study.\n"
        << "  --out <prefix>       Default: runtime/candidates/top4_study\n"
        << "  --report <path>      Default: runtime/reports/02_study_cpp.md\n\n"
        << "Local repair options:\n"
        << "  --iters <n>          Default: 50000\n"
        << "  --sigma <x>          Initial coordinate step. Default: 0.5\n"
        << "  --beta <x>           Step cooling. Default: 0.9995\n"
        << "  --temperature <x>    Annealing temperature. Default: 0.02\n"
        << "  --report-every <n>   Default: 1000\n"
        << "  --trials <n>         batch-hunt attempts. Default: 2000\n"
        << "  --threads <n>        batch-hunt worker threads. 0 = all cores.\n"
        << "  --minutes <x>        Clean time limit for batch-hunt. 0 = unlimited.\n"
        << "  --stop-file <path>   Stop cleanly when this file appears.\n"
        << "  --global-dir <path>  Checkpoints for --global-search. Default: results/search\n"
        << "  --topology-from <n>  First topology to search, inclusive. Default: 0\n"
        << "  --topology-to <n>    Last topology to search, inclusive. Default: 58\n"
        << "  --cuda                Require the FP32 CUDA search backend.\n"
        << "  --cuda-chains <n>     Parallel GPU chains. 0 = automatic (default).\n"
        << "  --cuda-iters <n>      Iterations per short GPU batch. Default: 64\n"
        << "  --checkpoint-seconds <n>  Durable checkpoint period. Default: 30\n"
        << "  --degeneracy-weight <x>   Worst-barrier degeneracy penalty. Default: 0.01\n"
        << "  --prioritize-worst    Give depth priority to the worst current topologies.\n"
        << "  --start-planes <p>   Continue batch-hunt from a saved .planes sidecar.\n"
        << "  --restarts <n>       hunt-local restarts. Default: 256\n"
        << "  --stagnation <n>     Iterations before hunt-local reheat. Default: 2000\n"
        << "  --jump-chance <x>    Large local jump probability. Default: 0.08\n"
        << "  --min-step-ratio <x> Smallest step/base-step ratio. Default: 1e-5\n";
}

#ifdef USE_CAIRO
void render_cutout(const Verts3D& v3ds, const Planes& planes, int width) {
    //Convert to 2D faces
    std::vector<Verts2D> faces(planes.size());
    double cur_x = 0.0;
    double cur_y = 0.0;
    double max_x = 0.0;
    double max_y = 0.0;
    for (size_t i = 0; i < planes.size(); ++i) {
        //Project points
        make_2d_projection(v3ds, g_polys[i], planes[i], faces[i]);

        //Figure out a bounding box
        Vector2d minCoord(1e9, 1e9);
        Vector2d maxCoord(-1e9, -1e9);
        for (Vector2d& v : faces[i]) {
            minCoord = minCoord.cwiseMin(v);
            maxCoord = maxCoord.cwiseMax(v);
        }

        //Transform coordinates
        for (Vector2d& v : faces[i]) {
            v -= minCoord;
            v.x() += cur_x;
            v.y() += cur_y;
        }

        //Advance height to next slot
        cur_y += maxCoord.y() - minCoord.y();
        max_x = std::max(max_x, maxCoord.x() - minCoord.x());

        if (i % 3 == 2) {
            cur_x += max_x;
            max_y = std::max(max_y, cur_y);
            cur_y = 0.0;
            max_x = 0.0;
        }
    }

    //Compute the scale factor
    const double padding = 4.0;
    const double scale = double(width - padding*2.0) / cur_x;

    //Create surface to draw on
    cairo_surface_t* surface = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, width, int(max_y * scale + padding * 2.0));
    cairo_t* cr = cairo_create(surface);

    cairo_set_line_width(cr, 2.0);
    cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);

    for (const Verts2D& face : faces) {
      const Vector2d& start_pt = face[face.size() - 1];
      cairo_move_to(cr, double(start_pt.x() * scale + padding), double(start_pt.y() * scale + padding));
      for (const Vector2d& p : face) {
        cairo_line_to(cr, double(p.x() * scale + padding), double(p.y() * scale + padding));
        cairo_stroke(cr);
        cairo_move_to(cr, double(p.x() * scale + padding), double(p.y() * scale + padding));
      }
    }

    //Save the image and free all memory
    cairo_surface_write_to_png(surface, "face1.png");
    cairo_destroy(cr);
    cairo_surface_destroy(surface);
}
#endif

void validate_files(bool change_files=true) {
    for (int i = 0; i < NUM_TOPOLOGIES; ++i) {
#if DUAL_PROBLEM != 0
        const std::filesystem::path dir("results/topologies/Dual/topology_" + std::to_string(i));
#else
        const std::filesystem::path dir("results/topologies/topology_" + std::to_string(i));
#endif
        for (const auto& dir_entry : std::filesystem::directory_iterator(dir)) {
            //Get intersections and crossings reported from file
            const auto& path = dir_entry.path();
            if (path.extension() != ".obj") { continue; }
            std::vector<std::string> name_split = split(path.filename().string(), '_');
            const int f_crossings = std::stoi(name_split[1].substr(1));
            const int f_intersections = std::stoi(name_split[2].substr(1));

            //Check for nan or inf in file
            Verts3D obj_verts;
            g_topology = i;
            import_obj(path.string().c_str(), obj_verts, g_polys);
            if (!is_finite(obj_verts)) {
                std::cout << path << std::endl;
                std::cout << "  NaNs found!" << std::endl << std::endl;
                if (change_files) { std::filesystem::remove(path); }
                continue;
            }

            //Get actual intersections and crossings
            Edges dual_edges;
            Planes obj_planes;
            make_edges(g_polys, g_edges);
            dual_graph(g_polys, g_tris, dual_edges);
            v3ds_to_planes(obj_verts, g_polys, obj_planes);
            const int crossings = count_crossings(obj_verts, obj_planes);
            const int intersections = count_intersections(obj_verts, obj_planes);

            //Check degenerate scores for study samples
            if (name_split[0] == "study") {
              const double dp = dist_penalty(obj_verts);
              const double lp = length_penalty(obj_verts);
              const double ap = angle_penalty(obj_verts);
              const double pp = plane_penalty(obj_planes);
              if (dp >= 0.9999 || ap >= 0.9999 || lp >= 0.9999 || pp >= 0.9999) {
                std::cout << path << std::endl;
                std::cout << "  DEGENERATE! dp(" << dp << ") dp(" << ap << ") lp(" << lp << ") pp(" << pp << ")" << std::endl << std::endl;
                if (change_files) { std::filesystem::remove(path); }
                continue;
              }
            }

            //Compare
            if (crossings != f_crossings || intersections != f_intersections) {
                const std::string new_name = name_split[0] + "_c" + std::to_string(crossings) + "_i" + std::to_string(intersections) + "_" + name_split[3];
                const std::filesystem::path new_path = std::filesystem::path(path).replace_filename(new_name);
                std::cout << path << std::endl;
                std::cout << new_path << std::endl << std::endl;
                if (change_files) { std::filesystem::rename(path, new_path); }
                continue;
            }
        }
    }
}

void main_solver() {
    int iter = 0;
    while (true) {
        //Load the dual adjacency list
        // g_topology = iter % NUM_TOPOLOGIES;
        open_topology("data/topologies.txt", g_tris, g_topology);
        std::cout << "Topology[" << g_topology << "]" << std::endl;

        //Create directory for results
#if DUAL_PROBLEM != 0
        const std::string topology_folder = "results/topologies/Dual/topology_" + std::to_string(g_topology);
#else
        const std::string topology_folder = "results/topologies/topology_" + std::to_string(g_topology);
#endif
        if (!std::filesystem::exists(topology_folder)) {
            std::filesystem::create_directories(topology_folder);
        }

        //Find the dual graph to get the polygon and edge linkage
        dual_graph(g_tris, g_polys, g_edges);
        fix_face_ordering(g_polys, g_edges);
#if DUAL_PROBLEM != 0
        std::swap(g_tris, g_polys);
        make_edges(g_polys, g_edges);
        fix_face_ordering(g_polys, g_edges);
#endif

        //Run the optimizer
        VectorXd result;
        double score = my_optimizer(objective_sum, result, 16000, 0.5, 0.998, 32, false, DUAL_PROBLEM);

        //Get the actual values of crossings and intersection independent of score
        Planes planes;
        Verts3D v3ds;
#if DUAL_PROBLEM != 0
        y_to_v3ds(result, v3ds);
        v3ds_to_planes(v3ds, g_polys, planes);
#else
        x_to_planes(result, planes);
        planes_to_v3ds(g_tris, planes, v3ds);
#endif
        const int crossings = count_crossings(v3ds, planes);
        const int intersections = count_intersections(v3ds, planes);

        //Check we should save it
        std::cout << "Score    : " << score << std::endl;
#if DUAL_PROBLEM != 0
        const bool can_save = (intersections <= 8);
#else
        const bool can_save = (crossings == 0 || (crossings + intersections <= 10));
#endif
        const std::string save_str = topology_folder + "/shape";
        save_sample(save_str.c_str(), planes, v3ds, iter, can_save);
        iter += 1;
        break;
    }
}

void quality_solver() {
    int iter = 0;
    while (true) {
        iter += 1;
        g_topology = iter % NUM_TOPOLOGIES;
#if DUAL_PROBLEM != 0
        const std::string topology_folder = "results/topologies/Dual/topology_" + std::to_string(g_topology);
#else
        const std::string topology_folder = "results/topologies/topology_" + std::to_string(g_topology);
#endif
        const std::filesystem::path dir(topology_folder);
        std::vector<std::filesystem::path> paths;
        for (const auto& dir_entry : std::filesystem::directory_iterator(dir)) {
            //Get intersections and crossings reported from file
            const auto& path = dir_entry.path();
            if (path.extension() != ".obj") { continue; }
            std::vector<std::string> name_split = split(path.filename().string(), '_');
            if (name_split[0] != "shape") { continue; }
            const int f_crossings = std::stoi(name_split[1].substr(1));
            const int f_intersections = std::stoi(name_split[2].substr(1));
            if ((f_crossings == 0 && f_intersections <= 12) ||
                (f_crossings <= 1 && f_intersections <= 8) ||
                (f_crossings <= 2 && f_intersections <= 6) ||
                (f_crossings <= 4 && f_intersections <= 4)) {
                paths.push_back(path);
            }
        }

        if (paths.size() == 0) { continue; }
        const std::filesystem::path& fpath = paths[std::uniform_int_distribution<int>(0, (int)paths.size() - 1)(eng)];

        Verts3D obj_verts;
        Edges dual_edges;
        Planes obj_planes;
        VectorXd obj_x;
        std::vector<std::string> name_split = split(fpath.stem().string(), '_');
        const int f_iter = std::atoi(name_split[name_split.size() - 1].c_str());
        import_obj(fpath.string().c_str(), obj_verts, g_polys);
        make_edges(g_polys, g_edges);
        dual_graph(g_polys, g_tris, dual_edges);
        v3ds_to_planes(obj_verts, g_polys, obj_planes);
#if DUAL_PROBLEM != 0
        v3ds_to_y(obj_verts, obj_x);
#else
        planes_to_x(obj_planes, obj_x);
#endif

        //Print characteristics
        std::cout << "===================" << std::endl;
        std::cout << "Loaded:        " << fpath << std::endl;
        save_sample("study", obj_planes, obj_verts, f_iter, false);
        std::cout << "===================" << std::endl;

        //Run optimizer
        study_sample(objective_dual_q, obj_x, 360, 10000, 1e-2, 0.9, true);
#if DUAL_PROBLEM != 0
        y_to_v3ds(obj_x, obj_verts);
        v3ds_to_planes(obj_verts, g_polys, obj_planes);
#else
        x_to_planes(obj_x, obj_planes);
        planes_to_v3ds(g_tris, obj_planes, obj_verts);
#endif

        if (q_penalty(obj_verts) >= 0.9999) {
            std::cout << "DEGENERATE" << std::endl;
        } else {
            const std::string out_path = fpath.parent_path().string() + "/study";
            save_sample(out_path.c_str(), obj_planes, obj_verts, f_iter, true);
            //render_cutout(obj_verts, obj_planes, 3840);
        }
    }
}

void explore_shape(const char* load_fname) {
    //Import an example obj file
    g_topology = 6;
    Verts3D obj_verts;
    Edges dual_edges;
    Planes obj_planes;
    VectorXd obj_x;
    std::vector<std::string> name_split = split(std::filesystem::path(load_fname).stem().string(), '_');
    const int f_iter = std::atoi(name_split[name_split.size() - 1].c_str());
    import_obj(load_fname, obj_verts, g_polys);
    make_edges(g_polys, g_edges);
    dual_graph(g_polys, g_tris, dual_edges);
    v3ds_to_planes(obj_verts, g_polys, obj_planes);
#if DUAL_PROBLEM != 0
    v3ds_to_y(obj_verts, obj_x);
#else
    planes_to_x(obj_planes, obj_x);
#endif

    //Print characteristics
    std::cout << "===================" << std::endl;
    std::cout << "Loaded:        " << load_fname << std::endl;
    save_sample("study", obj_planes, obj_verts, f_iter, false);
    std::cout << "===================" << std::endl;

    //Run optimizer
    study_sample(objective_sum_q, obj_x, 360, 10000, 1e-2, 0.9, true);
#if DUAL_PROBLEM != 0
    y_to_v3ds(obj_x, obj_verts);
    v3ds_to_planes(obj_verts, g_polys, obj_planes);
#else
    x_to_planes(obj_x, obj_planes);
    planes_to_v3ds(g_tris, obj_planes, obj_verts);
#endif
    save_sample("study", obj_planes, obj_verts, f_iter, true);
#ifdef USE_CAIRO
    render_cutout(obj_verts, obj_planes, 3840);
#endif
}

int study_shape(const StudyOptions& options) {
    double (*objective_function)(const VectorXd&) = select_objective(options.objective);
    if (objective_function == nullptr) {
        std::cerr << "Unknown objective: " << options.objective << std::endl;
        return 2;
    }

    g_topology = options.topology;
    set_rand_seed(options.seed);

    const std::filesystem::path out_prefix_path(options.out_prefix);
    if (out_prefix_path.has_parent_path()) {
        std::filesystem::create_directories(out_prefix_path.parent_path());
    }
    const std::filesystem::path report_path(options.report_path);
    if (report_path.has_parent_path()) {
        std::filesystem::create_directories(report_path.parent_path());
    }

    Verts3D obj_verts;
    Edges dual_edges;
    Planes obj_planes;
    VectorXd obj_x;

    import_obj(options.obj_path.c_str(), obj_verts, g_polys);
    if (obj_verts.empty() || g_polys.empty()) {
        std::cerr << "Failed to load OBJ: " << options.obj_path << std::endl;
        return 2;
    }

    make_edges(g_polys, g_edges);
    dual_graph(g_polys, g_tris, dual_edges);
    v3ds_to_planes(obj_verts, g_polys, obj_planes);
#if DUAL_PROBLEM != 0
    v3ds_to_y(obj_verts, obj_x);
#else
    planes_to_x(obj_planes, obj_x);
#endif

    const int initial_crossings = count_crossings(obj_verts, obj_planes);
    const int initial_intersections = count_intersections(obj_verts, obj_planes);
    const double initial_score = objective_function(obj_x);

    std::cout << "===================" << std::endl;
    std::cout << "Mode      : study" << std::endl;
    std::cout << "Loaded    : " << options.obj_path << std::endl;
    std::cout << "Seed      : " << options.seed << std::endl;
    std::cout << "Topology  : " << g_topology << std::endl;
    std::cout << "Objective : " << options.objective << std::endl;
    std::cout << "Iters     : " << options.max_iters << std::endl;
    std::cout << "Clusters  : " << options.clusters << std::endl;
    std::cout << "Sigma     : " << options.sigma << std::endl;
    std::cout << "Beta      : " << options.beta << std::endl;
    std::cout << "Symmetry  : " << (options.use_symmetry ? "yes" : "no") << std::endl;
    std::cout << "Initial crossing/intersect: " << initial_crossings << "/" << initial_intersections << std::endl;
    std::cout << "Initial objective: " << initial_score << std::endl;
    std::cout << "===================" << std::endl;
    save_sample(options.out_prefix.c_str(), obj_planes, obj_verts, options.seed, false);

    study_sample(
        objective_function,
        obj_x,
        options.max_iters,
        options.clusters,
        options.sigma,
        options.beta,
        options.use_symmetry
    );

#if DUAL_PROBLEM != 0
    y_to_v3ds(obj_x, obj_verts);
    v3ds_to_planes(obj_verts, g_polys, obj_planes);
#else
    x_to_planes(obj_x, obj_planes);
    planes_to_v3ds(g_tris, obj_planes, obj_verts);
#endif

    const int final_crossings = count_crossings(obj_verts, obj_planes);
    const int final_intersections = count_intersections(obj_verts, obj_planes);
    const double final_score = objective_function(obj_x);
    save_sample(options.out_prefix.c_str(), obj_planes, obj_verts, options.seed, true);

    if (final_crossings == 0 && final_intersections == 0) {
        export_obj("runtime/candidates/FOUND_top4_candidate.obj", obj_verts, g_polys);
    }

    std::ofstream report(options.report_path, std::ios::app);
    report << "# C++ study run\n\n";
    report << "- OBJ: `" << options.obj_path << "`\n";
    report << "- Seed: `" << options.seed << "`\n";
    report << "- Topology: `" << options.topology << "`\n";
    report << "- Objective: `" << options.objective << "`\n";
    report << "- Iters: `" << options.max_iters << "`\n";
    report << "- Clusters: `" << options.clusters << "`\n";
    report << "- Sigma/Beta: `" << options.sigma << "/" << options.beta << "`\n";
    report << "- Symmetry: `" << (options.use_symmetry ? "yes" : "no") << "`\n";
    report << "- Initial C/I: `" << initial_crossings << "/" << initial_intersections << "`\n";
    report << "- Final C/I: `" << final_crossings << "/" << final_intersections << "`\n";
    report << "- Initial objective: `" << initial_score << "`\n";
    report << "- Final objective: `" << final_score << "`\n";
    report << "- Output prefix: `" << options.out_prefix << "`\n\n";

    std::cout << "===================" << std::endl;
    std::cout << "Study finished." << std::endl;
    std::cout << "Initial C/I: " << initial_crossings << "/" << initial_intersections << std::endl;
    std::cout << "Final C/I  : " << final_crossings << "/" << final_intersections << std::endl;
    if (final_crossings == 0 && final_intersections == 0) {
        std::cout << "FOUND candidate saved to runtime/candidates/FOUND_top4_candidate.obj" << std::endl;
    }
    std::cout << "===================" << std::endl;
    return 0;
}

double cross2d(const Vector2d& a, const Vector2d& b) {
    return a.x() * b.y() - a.y() * b.x();
}

double normalized_orient2d(const Vector2d& a, const Vector2d& b, const Vector2d& c) {
    const Vector2d ab = b - a;
    const Vector2d ac = c - a;
    const double denom = std::max(1e-8, ab.norm() * ac.norm());
    return cross2d(ab, ac) / denom;
}

double median_edge_length(const Verts3D& verts) {
    std::vector<double> lengths;
    lengths.reserve(g_edges.size());
    for (const Edge& edge : g_edges) {
        lengths.push_back((verts[edge.first] - verts[edge.second]).norm());
    }
    if (lengths.empty()) {
        return 1.0;
    }
    std::sort(lengths.begin(), lengths.end());
    return std::max(1e-6, lengths[lengths.size() / 2]);
}

double smooth_crossing_loss(const Verts3D& verts, const Planes& planes) {
    double loss = 0.0;
    static const double margin = 1e-4;
    thread_local Verts2D projected;
    for (size_t face_ix = 0; face_ix < g_polys.size(); ++face_ix) {
        const Face& face = g_polys[face_ix];
        make_2d_projection(verts, face, planes[face_ix], projected);

        const size_t n = face.size();
        for (size_t i = 0; i < n; ++i) {
            const size_t i_next = (i + 1) % n;
            for (size_t j = i + 1; j < n; ++j) {
                const size_t j_next = (j + 1) % n;
                if (i_next == j || j_next == i) {
                    continue;
                }
                std::array<int, 4> endpoints = {face[i], face[i_next], face[j], face[j_next]};
                std::sort(endpoints.begin(), endpoints.end());
                if (std::unique(endpoints.begin(), endpoints.end()) != endpoints.end()) {
                    continue;
                }

                const Vector2d& a = projected[i];
                const Vector2d& b = projected[i_next];
                const Vector2d& c = projected[j];
                const Vector2d& d = projected[j_next];
                const double o1 = normalized_orient2d(a, b, c);
                const double o2 = normalized_orient2d(a, b, d);
                const double o3 = normalized_orient2d(c, d, a);
                const double o4 = normalized_orient2d(c, d, b);
                const double side_a = o1 * o2;
                const double side_b = o3 * o4;

                // A crossing requires both endpoint pairs to lie on opposite sides.
                // Penalize only that joint condition so ordinary concave geometry is not flattened.
                const double active_a = std::max(0.0, margin - side_a);
                const double active_b = std::max(0.0, margin - side_b);
                loss += active_a * active_b;
            }
        }
    }
    return loss;
}

double planarity_loss(const Verts3D& verts, double scale) {
    double loss = 0.0;
    for (const Face& face : g_polys) {
        const Plane plane = get_plane(verts, face);
        for (int vertex_ix : face) {
            const double d = plane.signed_distance(verts[vertex_ix]) / scale;
            loss += d * d;
        }
    }
    return loss / std::max<size_t>(1, g_polys.size());
}

double edge_length_loss(const Verts3D& verts, const Verts3D& original) {
    double loss = 0.0;
    for (const Edge& edge : g_edges) {
        const double current = (verts[edge.first] - verts[edge.second]).norm();
        const double base = std::max(1e-6, (original[edge.first] - original[edge.second]).norm());
        const double relative = current / base - 1.0;
        loss += relative * relative;
    }
    return loss / std::max<size_t>(1, g_edges.size());
}

double anchor_loss(const Verts3D& verts, const Verts3D& original, const std::vector<int>& movable, double scale) {
    double loss = 0.0;
    for (int vertex_ix : movable) {
        const double d = (verts[vertex_ix] - original[vertex_ix]).norm() / scale;
        loss += d * d;
    }
    return loss / std::max<size_t>(1, movable.size());
}

constexpr int PRECISE_DEFECT_THRESHOLD = 12;

void precise_counts(
    const Verts3D& verts,
    const Planes& planes,
    int& crossings,
    int& intersections,
    double eps = 1e-12
) {
    crossings = count_self_crossings_precise(verts, planes, eps);
    intersections = count_edge_face_intersections_precise(verts, planes, eps);
}

bool promote_precise_counts_if_close(
    const Verts3D& verts,
    const Planes& planes,
    int& crossings,
    int& intersections
) {
    if (crossings + intersections > PRECISE_DEFECT_THRESHOLD) {
        return false;
    }
    precise_counts(verts, planes, crossings, intersections);
    return true;
}

double local_repair_objective(
    const Verts3D& verts,
    const Verts3D& original,
    const std::vector<int>& movable,
    int& crossings,
    int& intersections
) {
    if (!is_finite(verts)) {
        return std::numeric_limits<double>::infinity();
    }

    Planes planes;
    v3ds_to_planes(verts, g_polys, planes);
    crossings = count_self_crossings_strict(verts, planes);
    intersections = count_edge_face_intersections_strict(verts, planes);
    promote_precise_counts_if_close(verts, planes, crossings, intersections);
    const double scale = median_edge_length(original);

    const double strict = double(crossings * 10000 + intersections * 10000);
    const double smooth = smooth_crossing_loss(verts, planes) * 2500.0;
    const double planar = planarity_loss(verts, scale) * 5000.0;
    const double lengths = edge_length_loss(verts, original) * 25.0;
    const double anchor = anchor_loss(verts, original, movable, scale) * 2.0;
    return strict + smooth + planar + lengths + anchor;
}

struct CrossingTarget {
    int face_ix;
    int a;
    int b;
    int c;
    int d;
    const char* label;
};

struct HuntMetrics {
    int target_crossings = 0;
    int crossings = -1;
    int intersections = -1;
    double target_loss = 0.0;
    double crossing_loss = 0.0;
    double regularization_loss = 0.0;
    double local_loss = 0.0;
    double total_loss = 0.0;
    bool verified = false;
};

const std::array<CrossingTarget, 2>& known_targets() {
    static const std::array<CrossingTarget, 2> targets = {{
        {4, 33, 32, 27, 4, "OBJ face 5: edge 34-33 vs 28-5"},
        {5, 33, 32, 12, 23, "OBJ face 6: edge 34-33 vs 13-24"}
    }};
    return targets;
}

int face_vertex_index(const Face& face, int vertex_ix) {
    for (size_t i = 0; i < face.size(); ++i) {
        if (face[i] == vertex_ix) {
            return static_cast<int>(i);
        }
    }
    return -1;
}

bool strict_segment_cross_2d(const Vector2d& a, const Vector2d& b, const Vector2d& c, const Vector2d& d) {
    static const double eps = 1e-7;
    const double o1 = cross2d(b - a, c - a);
    const double o2 = cross2d(b - a, d - a);
    const double o3 = cross2d(d - c, a - c);
    const double o4 = cross2d(d - c, b - c);
    return ((o1 > eps && o2 < -eps) || (o1 < -eps && o2 > eps)) &&
           ((o3 > eps && o4 < -eps) || (o3 < -eps && o4 > eps));
}

double crossing_barrier_2d(const Vector2d& a, const Vector2d& b, const Vector2d& c, const Vector2d& d) {
    static const double margin = 2e-4;
    const double o1 = normalized_orient2d(a, b, c);
    const double o2 = normalized_orient2d(a, b, d);
    const double o3 = normalized_orient2d(c, d, a);
    const double o4 = normalized_orient2d(c, d, b);
    const double side_a = o1 * o2;
    const double side_b = o3 * o4;
    const double active_a = std::max(0.0, margin - side_a);
    const double active_b = std::max(0.0, margin - side_b);
    return active_a * active_b;
}

std::vector<int> collect_affected_faces(const std::vector<int>& movable) {
    std::unordered_set<int> movable_set(movable.begin(), movable.end());
    std::vector<int> faces;
    for (size_t face_ix = 0; face_ix < g_polys.size(); ++face_ix) {
        for (int vertex_ix : g_polys[face_ix]) {
            if (movable_set.count(vertex_ix) != 0) {
                faces.push_back(static_cast<int>(face_ix));
                break;
            }
        }
    }
    return faces;
}

Edges collect_affected_edges(const std::vector<int>& movable) {
    std::unordered_set<int> movable_set(movable.begin(), movable.end());
    Edges edges;
    for (const Edge& edge : g_edges) {
        if (movable_set.count(edge.first) != 0 || movable_set.count(edge.second) != 0) {
            edges.push_back(edge);
        }
    }
    return edges;
}

double local_planarity_loss(const Verts3D& verts, const std::vector<int>& affected_faces, double scale) {
    double loss = 0.0;
    int samples = 0;
    for (int face_ix : affected_faces) {
        const Face& face = g_polys[face_ix];
        const Plane plane = get_plane(verts, face);
        for (int vertex_ix : face) {
            const double d = plane.signed_distance(verts[vertex_ix]) / scale;
            loss += d * d;
            samples += 1;
        }
    }
    return loss / std::max(1, samples);
}

double local_edge_length_loss(const Verts3D& verts, const Verts3D& original, const Edges& affected_edges) {
    double loss = 0.0;
    for (const Edge& edge : affected_edges) {
        const double current = (verts[edge.first] - verts[edge.second]).norm();
        const double base = std::max(1e-6, (original[edge.first] - original[edge.second]).norm());
        const double relative = current / base - 1.0;
        loss += relative * relative;
    }
    return loss / std::max<size_t>(1, affected_edges.size());
}

HuntMetrics evaluate_hunt_candidate(
    const Verts3D& verts,
    const Verts3D& original,
    const std::vector<int>& movable,
    const std::vector<int>& affected_faces,
    const Edges& affected_edges,
    double scale,
    bool verify_full
) {
    HuntMetrics metrics;
    if (!is_finite(verts)) {
        metrics.target_crossings = std::numeric_limits<int>::max() / 4;
        metrics.crossings = std::numeric_limits<int>::max() / 4;
        metrics.intersections = std::numeric_limits<int>::max() / 4;
        metrics.local_loss = std::numeric_limits<double>::infinity();
        metrics.total_loss = std::numeric_limits<double>::infinity();
        return metrics;
    }

    Verts2D projected;
    for (const CrossingTarget& target : known_targets()) {
        const Face& face = g_polys[target.face_ix];
        const Plane target_plane = get_plane(verts, face);
        make_2d_projection(verts, face, target_plane, projected);
        const int ia = face_vertex_index(face, target.a);
        const int ib = face_vertex_index(face, target.b);
        const int ic = face_vertex_index(face, target.c);
        const int id = face_vertex_index(face, target.d);
        if (ia < 0 || ib < 0 || ic < 0 || id < 0) {
            metrics.target_crossings += 1;
            metrics.target_loss += 10.0;
            continue;
        }
        const Vector2d& a = projected[ia];
        const Vector2d& b = projected[ib];
        const Vector2d& c = projected[ic];
        const Vector2d& d = projected[id];
        if (strict_segment_cross_2d(a, b, c, d)) {
            metrics.target_crossings += 1;
        }
        metrics.target_loss += crossing_barrier_2d(a, b, c, d);
    }

    const double planar = local_planarity_loss(verts, affected_faces, scale);
    const double lengths = local_edge_length_loss(verts, original, affected_edges);
    const double anchor = anchor_loss(verts, original, movable, scale);
    metrics.local_loss =
        metrics.target_crossings * 1000000.0 +
        metrics.target_loss * 250000.0 +
        planar * 2000.0 +
        lengths * 5.0 +
        anchor * 0.25;
    metrics.total_loss = metrics.local_loss;

    if (verify_full || metrics.target_crossings == 0) {
        Planes planes;
        v3ds_to_planes(verts, g_polys, planes);
        metrics.crossings = count_self_crossings_strict(verts, planes);
        metrics.intersections = count_edge_face_intersections_strict(verts, planes);
        promote_precise_counts_if_close(
            verts, planes, metrics.crossings, metrics.intersections);
        metrics.verified = true;
        metrics.total_loss += double(metrics.crossings * 10000 + metrics.intersections * 10000);
    }
    return metrics;
}

bool better_verified(const HuntMetrics& left, const HuntMetrics& right) {
    if (!left.verified) {
        return false;
    }
    if (!right.verified) {
        return true;
    }
    const auto rank = [](const HuntMetrics& value) {
        const int defects = value.crossings + value.intersections;
        return std::make_tuple(
            defects,
            std::max(value.crossings, value.intersections),
            value.target_crossings,
            value.intersections,
            value.total_loss);
    };
    return rank(left) < rank(right);
}

double search_loss(const HuntMetrics& metrics) {
    return metrics.local_loss;
}

bool better_fast(const HuntMetrics& left, const HuntMetrics& right) {
    if (!std::isfinite(search_loss(left))) {
        return false;
    }
    if (!std::isfinite(search_loss(right))) {
        return true;
    }
    if (left.target_crossings != right.target_crossings) {
        return left.target_crossings < right.target_crossings;
    }
    const int left_defects = left.crossings + left.intersections;
    const int right_defects = right.crossings + right.intersections;
    if (left_defects != right_defects) {
        return left_defects < right_defects;
    }
    const int left_peak = std::max(left.crossings, left.intersections);
    const int right_peak = std::max(right.crossings, right.intersections);
    if (left_peak != right_peak) {
        return left_peak < right_peak;
    }
    return search_loss(left) < search_loss(right);
}

bool better_search_state(const HuntMetrics& left, const HuntMetrics& right) {
    if (!std::isfinite(search_loss(left))) {
        return false;
    }
    if (!std::isfinite(search_loss(right))) {
        return true;
    }
    const int left_defects = left.crossings + left.intersections;
    const int right_defects = right.crossings + right.intersections;
    if (left_defects != right_defects) {
        return left_defects < right_defects;
    }
    if (left.target_crossings != right.target_crossings) {
        return left.target_crossings < right.target_crossings;
    }
    return search_loss(left) < search_loss(right);
}

Vector3d random_vec3(std::normal_distribution<double>& normal) {
    return Vector3d(normal(eng), normal(eng), normal(eng));
}

Vector3d random_vec3(std::normal_distribution<double>& normal, RNG& rng) {
    return Vector3d(normal(rng), normal(rng), normal(rng));
}

void apply_hunt_move(
    Verts3D& candidate,
    const std::vector<int>& movable,
    double step,
    bool large_jump,
    std::normal_distribution<double>& normal,
    std::uniform_real_distribution<double>& uniform
) {
    const double s = step * (large_jump ? 4.0 : 1.0);
    const int move_type = std::uniform_int_distribution<int>(0, 6)(eng);
    auto add = [&](int ix, const Vector3d& delta) {
        candidate[ix] += delta;
    };

    if (move_type == 0) {
        const int vertex_ix = movable[std::uniform_int_distribution<int>(0, static_cast<int>(movable.size()) - 1)(eng)];
        add(vertex_ix, random_vec3(normal) * s);
    } else if (move_type == 1) {
        const Vector3d delta = random_vec3(normal) * s;
        add(32, delta);
        add(33, delta);
    } else if (move_type == 2) {
        const Vector3d delta = random_vec3(normal) * s;
        add(32, delta);
        add(33, -delta);
    } else if (move_type == 3) {
        const Vector3d delta = random_vec3(normal) * s;
        add(27, delta);
        add(4, delta);
    } else if (move_type == 4) {
        const Vector3d delta = random_vec3(normal) * s;
        add(12, delta);
        add(23, delta);
    } else if (move_type == 5) {
        for (int vertex_ix : movable) {
            add(vertex_ix, random_vec3(normal) * (s * 0.35));
        }
    } else {
        const Vector3d common = random_vec3(normal) * (s * 0.7);
        add(32, common + random_vec3(normal) * (s * 0.25));
        add(33, common + random_vec3(normal) * (s * 0.25));
        add(27, -common * 0.5 + random_vec3(normal) * (s * 0.25));
        add(4, -common * 0.5 + random_vec3(normal) * (s * 0.25));
        add(12, -common * 0.5 + random_vec3(normal) * (s * 0.25));
        add(23, -common * 0.5 + random_vec3(normal) * (s * 0.25));
    }

    if (large_jump && uniform(eng) < 0.35) {
        const Vector3d twist = random_vec3(normal) * (s * 0.75);
        add(32, twist);
        add(33, -twist);
    }
}

void jitter_seed_shape(
    Verts3D& verts,
    const std::vector<int>& movable,
    double step,
    std::normal_distribution<double>& normal
) {
    for (int vertex_ix : movable) {
        verts[vertex_ix] += random_vec3(normal) * step;
    }
}

double plane_anchor_loss(const VectorXd& x, const VectorXd& original_x) {
    if (x.size() == 0) {
        return 0.0;
    }
    return (x - original_x).squaredNorm() / static_cast<double>(x.size());
}

struct PlaneEvaluationScratch {
    Planes planes;
    Planes canonical_planes;
    Verts3D verts;
};

HuntMetrics invalid_hunt_metrics() {
    HuntMetrics metrics;
    metrics.target_crossings = std::numeric_limits<int>::max() / 4;
    metrics.crossings = std::numeric_limits<int>::max() / 4;
    metrics.intersections = std::numeric_limits<int>::max() / 4;
    metrics.target_loss = std::numeric_limits<double>::infinity();
    metrics.crossing_loss = std::numeric_limits<double>::infinity();
    metrics.regularization_loss = std::numeric_limits<double>::infinity();
    metrics.local_loss = std::numeric_limits<double>::infinity();
    metrics.total_loss = std::numeric_limits<double>::infinity();
    return metrics;
}

bool valid_plane_state(const VectorXd& x) {
    if (x.size() != static_cast<int>(g_polys.size() * 3) || !x.allFinite()) {
        return false;
    }
    for (int i = 0; i < x.size(); i += 3) {
        if (Eigen::Map<const Vector3d>(x.data() + i).squaredNorm() < 1e-20) {
            return false;
        }
    }
    return true;
}

HuntMetrics evaluate_hunt_planes_with_scratch(
    const VectorXd& x,
    const VectorXd& original_x,
    const Verts3D& original_verts,
    const Edges& affected_edges,
    bool canonical,
    PlaneEvaluationScratch& scratch
) {
    if (!valid_plane_state(x)) {
        return invalid_hunt_metrics();
    }

    scratch.planes.reserve(g_polys.size());
    scratch.canonical_planes.reserve(g_polys.size());
    scratch.verts.reserve(g_tris.size());
    x_to_planes(x, scratch.planes);
    planes_to_v3ds(g_tris, scratch.planes, scratch.verts);
    if (!is_finite(scratch.verts)) {
        return invalid_hunt_metrics();
    }

    const Planes* metric_planes = &scratch.planes;
    if (canonical) {
        v3ds_to_planes(scratch.verts, g_polys, scratch.canonical_planes);
        metric_planes = &scratch.canonical_planes;
    }

    HuntMetrics metrics;
    thread_local Verts2D projected;
    for (const CrossingTarget& target : known_targets()) {
        const Face& face = g_polys[target.face_ix];
        make_2d_projection(
            scratch.verts, face, (*metric_planes)[target.face_ix], projected);
        const int ia = face_vertex_index(face, target.a);
        const int ib = face_vertex_index(face, target.b);
        const int ic = face_vertex_index(face, target.c);
        const int id = face_vertex_index(face, target.d);
        if (ia < 0 || ib < 0 || ic < 0 || id < 0) {
            metrics.target_crossings += 1;
            metrics.target_loss += 10.0;
            continue;
        }
        const Vector2d& a = projected[ia];
        const Vector2d& b = projected[ib];
        const Vector2d& c = projected[ic];
        const Vector2d& d = projected[id];
        if (strict_segment_cross_2d(a, b, c, d)) {
            metrics.target_crossings += 1;
        }
        metrics.target_loss += crossing_barrier_2d(a, b, c, d);
    }

    metrics.crossings = count_self_crossings_strict(
        scratch.verts, *metric_planes, 1e-8, &metrics.crossing_loss);
    metrics.intersections = count_edge_face_intersections_strict(scratch.verts, *metric_planes);
    promote_precise_counts_if_close(
        scratch.verts, *metric_planes, metrics.crossings, metrics.intersections);

    const double lengths = local_edge_length_loss(scratch.verts, original_verts, affected_edges);
    const double plane_anchor = plane_anchor_loss(x, original_x);
    metrics.regularization_loss =
        0.0005 * std::min(lengths, 100.0) +
        0.00005 * std::min(plane_anchor, 1000.0);

    // Keep every term near unit scale so the GUI temperature has a real
    // annealing effect. One temporary defect costs 0.05 instead of 10000.
    const int defects = metrics.crossings + metrics.intersections;
    metrics.local_loss =
        0.05 * static_cast<double>(defects) +
        0.01 * static_cast<double>(metrics.crossings) +
        0.0125 * static_cast<double>(metrics.target_crossings) +
        0.20 * std::min(metrics.target_loss, 10.0) +
        0.02 * std::min(metrics.crossing_loss, 10.0) +
        metrics.regularization_loss;
    metrics.total_loss = metrics.local_loss;
    metrics.verified = canonical;
    return metrics;
}

HuntMetrics evaluate_hunt_planes(
    const VectorXd& x,
    const VectorXd& original_x,
    const Verts3D& original_verts,
    const Edges& affected_edges,
    bool canonical
) {
    PlaneEvaluationScratch scratch;
    return evaluate_hunt_planes_with_scratch(
        x, original_x, original_verts, affected_edges, canonical, scratch);
}

void apply_plane_hunt_move(
    VectorXd& x,
    const std::vector<int>& movable_faces,
    double step,
    bool large_jump,
    std::normal_distribution<double>& normal
) {
    const double s = step * (large_jump ? 4.0 : 1.0);
    const int move_type = std::uniform_int_distribution<int>(0, 5)(eng);
    auto add_plane = [&](int face_ix, const Vector3d& delta) {
        Eigen::Map<Vector3d>(x.data() + face_ix * 3) += delta;
    };

    if (move_type == 0) {
        const int face_ix = movable_faces[std::uniform_int_distribution<int>(0, static_cast<int>(movable_faces.size()) - 1)(eng)];
        add_plane(face_ix, random_vec3(normal) * s);
    } else if (move_type == 1) {
        const Vector3d delta = random_vec3(normal) * s;
        add_plane(4, delta);
        add_plane(5, delta);
    } else if (move_type == 2) {
        const Vector3d delta = random_vec3(normal) * s;
        add_plane(4, delta);
        add_plane(5, -delta);
    } else if (move_type == 3) {
        add_plane(4, random_vec3(normal) * (s * 1.4));
    } else if (move_type == 4) {
        add_plane(5, random_vec3(normal) * (s * 1.4));
    } else {
        for (int face_ix : movable_faces) {
            add_plane(face_ix, random_vec3(normal) * (s * 0.20));
        }
    }
}

void apply_plane_hunt_move(
    VectorXd& x,
    const std::vector<int>& movable_faces,
    double step,
    bool large_jump,
    std::normal_distribution<double>& normal,
    RNG& rng
) {
    const double s = step * (large_jump ? 4.0 : 1.0);
    const int move_type = std::uniform_int_distribution<int>(0, 5)(rng);
    auto add_plane = [&](int face_ix, const Vector3d& delta) {
        Eigen::Map<Vector3d>(x.data() + face_ix * 3) += delta;
    };

    if (move_type == 0) {
        const int face_ix = movable_faces[std::uniform_int_distribution<int>(0, static_cast<int>(movable_faces.size()) - 1)(rng)];
        add_plane(face_ix, random_vec3(normal, rng) * s);
    } else if (move_type == 1) {
        const Vector3d delta = random_vec3(normal, rng) * s;
        add_plane(4, delta);
        add_plane(5, delta);
    } else if (move_type == 2) {
        const Vector3d delta = random_vec3(normal, rng) * s;
        add_plane(4, delta);
        add_plane(5, -delta);
    } else if (move_type == 3) {
        add_plane(4, random_vec3(normal, rng) * (s * 1.4));
    } else if (move_type == 4) {
        add_plane(5, random_vec3(normal, rng) * (s * 1.4));
    } else {
        for (int face_ix : movable_faces) {
            add_plane(face_ix, random_vec3(normal, rng) * (s * 0.20));
        }
    }
}

void jitter_seed_planes(
    VectorXd& x,
    const std::vector<int>& movable_faces,
    double step,
    std::normal_distribution<double>& normal
) {
    for (int face_ix : movable_faces) {
        Eigen::Map<Vector3d>(x.data() + face_ix * 3) += random_vec3(normal) * (step * 0.35);
    }
}

void jitter_seed_planes(
    VectorXd& x,
    const std::vector<int>& movable_faces,
    double step,
    std::normal_distribution<double>& normal,
    RNG& rng
) {
    for (int face_ix : movable_faces) {
        Eigen::Map<Vector3d>(x.data() + face_ix * 3) += random_vec3(normal, rng) * (step * 0.35);
    }
}

int hunt_local_shape(const LocalRepairOptions& options) {
    g_topology = options.topology;
    set_rand_seed(options.seed);

    const std::filesystem::path out_prefix_path(options.out_prefix);
    if (out_prefix_path.has_parent_path()) {
        std::filesystem::create_directories(out_prefix_path.parent_path());
    }
    const std::filesystem::path report_path(options.report_path);
    if (report_path.has_parent_path()) {
        std::filesystem::create_directories(report_path.parent_path());
    }

    Verts3D original;
    import_obj(options.obj_path.c_str(), original, g_polys);
    if (original.empty() || g_polys.empty()) {
        std::cerr << "Failed to load OBJ: " << options.obj_path << std::endl;
        return 2;
    }
    make_edges(g_polys, g_edges);
    Edges dual_edges;
    dual_graph(g_polys, g_tris, dual_edges);

    Planes original_planes;
    VectorXd original_x;
    if (!options.start_planes_path.empty()) {
        if (!load_plane_state(options.start_planes_path, original_x) ||
            !valid_plane_state(original_x)) {
            std::cerr << "Failed to load plane state: " << options.start_planes_path << std::endl;
            return 2;
        }
        x_to_planes(original_x, original_planes);
        planes_to_v3ds(g_tris, original_planes, original);
        if (!is_finite(original)) {
            std::cerr << "Plane state reconstructs non-finite vertices: "
                      << options.start_planes_path << std::endl;
            return 2;
        }
    } else {
        v3ds_to_planes(original, g_polys, original_planes);
        planes_to_x(original_planes, original_x);
    }

    std::vector<int> movable_faces = collect_affected_faces(options.movable_vertices);
    movable_faces.push_back(4);
    movable_faces.push_back(5);
    std::sort(movable_faces.begin(), movable_faces.end());
    movable_faces.erase(std::unique(movable_faces.begin(), movable_faces.end()), movable_faces.end());
    const Edges affected_edges = collect_affected_edges(options.movable_vertices);

    VectorXd current = original_x;
    VectorXd best_verified_x = original_x;
    VectorXd best_fast_x = original_x;
    HuntMetrics current_metrics = evaluate_hunt_planes(
        current, original_x, original, affected_edges, true);
    HuntMetrics best_verified = evaluate_hunt_planes(
        original_x, original_x, original, affected_edges, true);
    HuntMetrics best_fast = current_metrics;

    double step = options.step;
    double temperature = options.temperature;
    const double min_step = std::max(1e-8, options.step * options.min_step_ratio);
    int stagnant = 0;
    int restart_count = 0;

    std::normal_distribution<double> normal(0.0, 1.0);
    std::uniform_real_distribution<double> uniform(0.0, 1.0);

    std::ofstream report(options.report_path, std::ios::app);
    report << "# C++ focused plane hunt run\n\n";
    report << "- OBJ: `" << options.obj_path << "`\n";
    if (!options.start_planes_path.empty()) {
        report << "- Start planes: `" << options.start_planes_path << "`\n";
    }
    report << "- Seed: `" << options.seed << "`\n";
    report << "- Topology: `" << options.topology << "`\n";
    report << "- Iterations: `" << options.iterations << "`\n";
    report << "- Initial step: `" << options.step << "`\n";
    report << "- Min step: `" << min_step << "`\n";
    report << "- Restarts: `" << options.restarts << "`\n";
    report << "- Search space: `12 face planes; local affected plane subset`\n";
    report << "- Targets: `34-33 x 28-5`, `34-33 x 13-24`\n";
    report << "- Initial verified C/I: `" << best_verified.crossings << "/" << best_verified.intersections << "`\n\n";

    std::cout << "===================" << std::endl;
    std::cout << "Mode      : hunt-local" << std::endl;
    std::cout << "Space     : face planes" << std::endl;
    std::cout << "Loaded    : " << options.obj_path << std::endl;
    if (!options.start_planes_path.empty()) {
        std::cout << "Planes    : " << options.start_planes_path << std::endl;
    }
    std::cout << "Seed      : " << options.seed << std::endl;
    std::cout << "Topology  : " << g_topology << std::endl;
    std::cout << "Iters     : " << options.iterations << std::endl;
    std::cout << "Step      : " << options.step << " floor " << min_step << std::endl;
    std::cout << "Targets   : edge 34-33 against 28-5 and 13-24" << std::endl;
    std::cout << "Initial verified C/I: " << best_verified.crossings << "/" << best_verified.intersections << std::endl;
    std::cout << "===================" << std::endl;

    for (int iter = 1; iter <= options.iterations; ++iter) {
        VectorXd candidate = current;
        const bool large_jump = uniform(eng) < options.jump_chance;
        apply_plane_hunt_move(candidate, movable_faces, step, large_jump, normal);

        HuntMetrics candidate_metrics = evaluate_hunt_planes(
            candidate, original_x, original, affected_edges, true);

        const bool improves_current = search_loss(candidate_metrics) < search_loss(current_metrics);
        const double accept_probability = std::exp(
            (search_loss(current_metrics) - search_loss(candidate_metrics)) / std::max(1e-6, temperature));
        if (improves_current || uniform(eng) < accept_probability) {
            current = candidate;
            current_metrics = candidate_metrics;
        }

        bool improved = false;
        if (better_fast(candidate_metrics, best_fast)) {
            best_fast = candidate_metrics;
            best_fast_x = candidate;
            improved = true;
            std::cout << "Fast best iter " << iter
                      << ": target " << best_fast.target_crossings
                      << ", loss " << search_loss(best_fast)
                      << ", step " << step << std::endl;
            if (best_fast.verified) {
                std::cout << "  Fast best verified C/I "
                          << best_fast.crossings << "/" << best_fast.intersections << std::endl;
            }
        }

        if (candidate_metrics.verified && better_verified(candidate_metrics, best_verified)) {
            best_verified = candidate_metrics;
            best_verified_x = candidate;
            improved = true;
            std::cout << "Verified best iter " << iter << ": C/I "
                      << best_verified.crossings << "/" << best_verified.intersections
                      << ", target " << best_verified.target_crossings
                      << ", loss " << best_verified.total_loss
                      << ", step " << step << std::endl;
            report << "- Verified best iter `" << iter << "`: C/I `"
                   << best_verified.crossings << "/" << best_verified.intersections
                   << "`, target `" << best_verified.target_crossings
                   << "`, loss `" << best_verified.total_loss << "`\n";
        }

        stagnant = improved ? 0 : stagnant + 1;
        step = std::max(min_step, step * options.beta);
        temperature = std::max(options.temperature * 0.05, temperature * options.beta);

        if (best_verified.verified && best_verified.crossings == 0 && best_verified.intersections == 0) {
            break;
        }

        if (stagnant >= options.stagnation && restart_count < options.restarts) {
            restart_count += 1;
            const bool use_best = uniform(eng) < 0.70;
            current = use_best ? best_verified_x : original_x;
            const double restart_step = options.step * (0.35 + uniform(eng) * 1.25);
            jitter_seed_planes(current, movable_faces, restart_step, normal);
            current_metrics = evaluate_hunt_planes(
                current, original_x, original, affected_edges, true);
            step = std::max(min_step, restart_step);
            temperature = options.temperature * (0.75 + uniform(eng));
            stagnant = 0;
            std::cout << "Restart " << restart_count
                      << ": target " << current_metrics.target_crossings
                      << ", step " << step << std::endl;
        }

        if (iter % options.report_every == 0) {
            std::cout << "Iter " << iter
                      << ": current target " << current_metrics.target_crossings
                      << ", best verified " << best_verified.crossings << "/" << best_verified.intersections
                      << ", best target " << best_fast.target_crossings
                      << ", restarts " << restart_count
                      << ", step " << step << std::endl;
        }
    }

    const std::string verified_prefix = options.out_prefix + "_verified";
    const std::string fast_prefix = options.out_prefix + "_fast";

    Planes fast_planes;
    Verts3D best_fast_shape;
    x_to_planes(best_fast_x, fast_planes);
    planes_to_v3ds(g_tris, fast_planes, best_fast_shape);
    save_sample(fast_prefix.c_str(), fast_planes, best_fast_shape, options.seed, true);

    Planes best_planes;
    Verts3D best_verified_shape;
    x_to_planes(best_verified_x, best_planes);
    planes_to_v3ds(g_tris, best_planes, best_verified_shape);
    save_sample(verified_prefix.c_str(), best_planes, best_verified_shape, options.seed, true);
    if (best_verified.crossings == 0 && best_verified.intersections == 0) {
        export_obj("runtime/candidates/FOUND_top4_candidate.obj", best_verified_shape, g_polys);
    }

    report << "\n- Final verified C/I: `" << best_verified.crossings << "/" << best_verified.intersections << "`\n";
    report << "- Final target crossings: `" << best_verified.target_crossings << "`\n";
    report << "- Fast best target crossings: `" << best_fast.target_crossings << "`\n";
    if (best_fast.verified) {
        report << "- Fast best verified C/I: `" << best_fast.crossings << "/" << best_fast.intersections << "`\n";
    }
    report << "- Restarts used: `" << restart_count << "`\n";
    report << "- Output prefix: `" << options.out_prefix << "`\n\n";

    std::cout << "===================" << std::endl;
    std::cout << "Focused hunt finished." << std::endl;
    std::cout << "Best verified C/I: " << best_verified.crossings << "/" << best_verified.intersections << std::endl;
    std::cout << "Best target crossings: " << best_fast.target_crossings << std::endl;
    if (best_fast.verified) {
        std::cout << "Fast candidate verified C/I: " << best_fast.crossings << "/" << best_fast.intersections << std::endl;
    }
    if (best_verified.crossings == 0 && best_verified.intersections == 0) {
        std::cout << "FOUND candidate saved to runtime/candidates/FOUND_top4_candidate.obj" << std::endl;
    }
    std::cout << "===================" << std::endl;
    return 0;
}

struct BatchTrialResult {
    int trial = 0;
    int seed = 0;
    double step = 0.0;
    double beta = 0.0;
    double temperature = 0.0;
    double jump_chance = 0.0;
    int face_count = 0;
    int restarts = 0;
    HuntMetrics best_verified;
    HuntMetrics best_fast;
    VectorXd best_verified_x;
    VectorXd best_fast_x;
};

std::vector<int> unique_faces(std::vector<int> faces) {
    std::sort(faces.begin(), faces.end());
    faces.erase(std::unique(faces.begin(), faces.end()), faces.end());
    return faces;
}

std::vector<int> all_face_indices() {
    std::vector<int> faces;
    faces.reserve(g_polys.size());
    for (int i = 0; i < static_cast<int>(g_polys.size()); ++i) {
        faces.push_back(i);
    }
    return faces;
}

double random_log_scale(RNG& rng, double lo_exp, double hi_exp) {
    std::uniform_real_distribution<double> uniform(0.0, 1.0);
    return std::pow(10.0, lo_exp + (hi_exp - lo_exp) * uniform(rng));
}

void save_plane_candidate(const std::string& prefix, const VectorXd& x, int seed) {
    Planes planes;
    Verts3D verts;
    x_to_planes(x, planes);
    planes_to_v3ds(g_tris, planes, verts);
    save_sample(prefix.c_str(), planes, verts, seed, true);
}

bool save_plane_state(const std::filesystem::path& path, const VectorXd& x) {
    std::ofstream out(path);
    if (!out) {
        return false;
    }
    out << std::setprecision(17);
    out << x.size() << "\n";
    for (int i = 0; i < x.size(); ++i) {
        out << static_cast<double>(x[i]) << "\n";
    }
    return true;
}

bool load_plane_state(const std::filesystem::path& path, VectorXd& x) {
    std::ifstream in(path);
    if (!in) {
        return false;
    }
    int size = 0;
    in >> size;
    if (size <= 0 || size % 3 != 0) {
        return false;
    }
    x.resize(size);
    for (int i = 0; i < size; ++i) {
        double value = 0.0;
        in >> value;
        if (!in) {
            return false;
        }
        x[i] = static_cast<double>(value);
    }
    return true;
}

bool roundtrip_plane_state(const VectorXd& source, VectorXd& result) {
    std::ostringstream out;
    out << std::setprecision(17);
    out << source.size() << "\n";
    for (int i = 0; i < source.size(); ++i) {
        out << static_cast<double>(source[i]) << "\n";
    }

    std::istringstream in(out.str());
    int size = 0;
    in >> size;
    if (size <= 0 || size % 3 != 0) {
        return false;
    }
    result.resize(size);
    for (int i = 0; i < size; ++i) {
        double value = 0.0;
        in >> value;
        if (!in) {
            return false;
        }
        result[i] = static_cast<double>(value);
    }
    return true;
}

void save_plane_candidate_quiet(const std::string& prefix, const VectorXd& x, int seed, const HuntMetrics& metrics) {
    VectorXd stored_x;
    if (!roundtrip_plane_state(x, stored_x)) {
        stored_x = x;
    }
    Planes planes;
    Verts3D verts;
    x_to_planes(stored_x, planes);
    planes_to_v3ds(g_tris, planes, verts);

    HuntMetrics stored_metrics = metrics;
    Planes verification_planes;
    v3ds_to_planes(verts, g_polys, verification_planes);
    precise_counts(
        verts, verification_planes,
        stored_metrics.crossings, stored_metrics.intersections);
    stored_metrics.verified = true;

    std::ostringstream name;
    name << prefix
         << "_strict_c" << stored_metrics.crossings
         << "_i" << stored_metrics.intersections
         << "_" << seed
         << ".obj";
    const std::filesystem::path path(name.str());
    if (path.has_parent_path()) {
        std::filesystem::create_directories(path.parent_path());
    }
    export_obj(path.string().c_str(), verts, g_polys);
    std::filesystem::path state_path = path;
    state_path.replace_extension(".planes");
    save_plane_state(state_path, stored_x);
}

void export_plane_candidate(const char* path, const VectorXd& x) {
    Planes planes;
    Verts3D verts;
    x_to_planes(x, planes);
    planes_to_v3ds(g_tris, planes, verts);
    export_obj(path, verts, g_polys);
}

bool export_and_validate_found_candidate(
    const std::filesystem::path& path,
    const VectorXd& x,
    int& crossings,
    int& intersections
) {
    std::filesystem::path temp_path = path;
    temp_path += ".tmp";
    export_plane_candidate(temp_path.string().c_str(), x);
    const auto remove_temp = [&]() {
        std::error_code remove_error;
        std::filesystem::remove(temp_path, remove_error);
    };

    Verts3D saved_verts;
    Faces saved_faces;
    import_obj(temp_path.string().c_str(), saved_verts, saved_faces);
    if (saved_verts.size() != g_tris.size() || saved_faces != g_polys || !is_finite(saved_verts)) {
        crossings = std::numeric_limits<int>::max() / 4;
        intersections = std::numeric_limits<int>::max() / 4;
        remove_temp();
        return false;
    }

    Planes saved_planes;
    v3ds_to_planes(saved_verts, saved_faces, saved_planes);
    const double scale = median_edge_length(saved_verts);
    if (!std::isfinite(scale) || scale <= 0.0) {
        remove_temp();
        return false;
    }
    for (const Edge& edge : g_edges) {
        const double length = (saved_verts[edge.first] - saved_verts[edge.second]).norm();
        if (!std::isfinite(length) || length <= scale * 1e-8) {
            remove_temp();
            return false;
        }
    }
    for (size_t face_ix = 0; face_ix < saved_faces.size(); ++face_ix) {
        for (int vertex_ix : saved_faces[face_ix]) {
            if (std::abs(saved_planes[face_ix].signed_distance(saved_verts[vertex_ix])) > scale * 1e-4) {
                remove_temp();
                return false;
            }
        }
    }
    precise_counts(saved_verts, saved_planes, crossings, intersections);
    for (double eps : {1e-7, 1e-9, 1e-11, 1e-13}) {
        if (count_self_crossings_precise(saved_verts, saved_planes, eps) != 0 ||
            count_edge_face_intersections_precise(saved_verts, saved_planes, eps) != 0) {
            remove_temp();
            return false;
        }
    }

    std::error_code copy_error;
    std::filesystem::copy_file(
        temp_path,
        path,
        std::filesystem::copy_options::overwrite_existing,
        copy_error);
    remove_temp();
    if (copy_error) {
        return false;
    }

    std::filesystem::path state_path = path;
    state_path.replace_extension(".planes");
    return save_plane_state(state_path, x);
}

BatchTrialResult run_plane_batch_trial(
    const LocalRepairOptions& trial_options,
    const Verts3D& original,
    const VectorXd& original_x,
    const VectorXd& start_x,
    const Edges& affected_edges,
    const std::vector<int>& movable_faces,
    int trial_ix,
    const std::atomic<bool>& stop_requested
) {
    BatchTrialResult result;
    result.trial = trial_ix;
    result.seed = trial_options.seed;
    result.step = trial_options.step;
    result.beta = trial_options.beta;
    result.temperature = trial_options.temperature;
    result.jump_chance = trial_options.jump_chance;
    result.face_count = static_cast<int>(movable_faces.size());

    RNG trial_rng(trial_options.seed);
    std::normal_distribution<double> normal(0.0, 1.0);
    std::uniform_real_distribution<double> uniform(0.0, 1.0);
    PlaneEvaluationScratch scratch;

    VectorXd current = start_x;
    if (trial_ix > 1) {
        jitter_seed_planes(current, movable_faces, trial_options.step * (0.25 + uniform(trial_rng)), normal, trial_rng);
    }

    HuntMetrics current_metrics = evaluate_hunt_planes_with_scratch(
        current, original_x, original, affected_edges, false, scratch);
    HuntMetrics current_canonical = evaluate_hunt_planes_with_scratch(
        current, original_x, original, affected_edges, true, scratch);
    HuntMetrics best_search = current_metrics;
    VectorXd best_search_x = current;
    result.best_verified = current_canonical;
    result.best_fast = current_canonical;
    result.best_verified_x = current;
    result.best_fast_x = current;

    double step = trial_options.step;
    double temperature = trial_options.temperature;
    const double min_step = std::max(1e-8, trial_options.step * trial_options.min_step_ratio);
    int stagnant = 0;

    for (int iter = 1; iter <= trial_options.iterations; ++iter) {
        if ((iter & 255) == 0 && stop_requested.load(std::memory_order_relaxed)) {
            break;
        }
        VectorXd candidate = current;
        const bool large_jump = uniform(trial_rng) < trial_options.jump_chance;
        apply_plane_hunt_move(candidate, movable_faces, step, large_jump, normal, trial_rng);

        HuntMetrics candidate_metrics = evaluate_hunt_planes_with_scratch(
            candidate, original_x, original, affected_edges, false, scratch);

        const bool improves_current = search_loss(candidate_metrics) < search_loss(current_metrics);
        const double exponent = std::clamp(
            (search_loss(current_metrics) - search_loss(candidate_metrics)) /
                std::max(1e-6, temperature),
            -80.0,
            0.0);
        const double accept_probability = std::exp(exponent);
        if (improves_current || uniform(trial_rng) < accept_probability) {
            current = candidate;
            current_metrics = candidate_metrics;
        }

        bool improved = false;
        const bool search_improved = better_search_state(candidate_metrics, best_search);
        const bool target_improved = better_fast(candidate_metrics, result.best_fast);
        const bool possible_solution =
            candidate_metrics.crossings == 0 && candidate_metrics.intersections == 0;
        if (search_improved) {
            best_search = candidate_metrics;
            best_search_x = candidate;
            improved = true;
        }
        if (search_improved || target_improved || possible_solution) {
            HuntMetrics canonical_metrics = evaluate_hunt_planes_with_scratch(
                candidate, original_x, original, affected_edges, true, scratch);
            if (better_fast(canonical_metrics, result.best_fast)) {
                result.best_fast = canonical_metrics;
                result.best_fast_x = candidate;
                improved = true;
            }
            if (better_verified(canonical_metrics, result.best_verified)) {
                result.best_verified = canonical_metrics;
                result.best_verified_x = candidate;
                improved = true;
            }
        }

        stagnant = improved ? 0 : stagnant + 1;
        step = std::max(min_step, step * trial_options.beta);
        temperature = std::max(trial_options.temperature * 0.05, temperature * trial_options.beta);

        if (result.best_verified.verified &&
            result.best_verified.crossings == 0 &&
            result.best_verified.intersections == 0) {
            break;
        }

        if (stagnant >= trial_options.stagnation && result.restarts < trial_options.restarts) {
            result.restarts += 1;
            const double restart_source = uniform(trial_rng);
            if (restart_source < 0.55) {
                current = best_search_x;
            } else if (restart_source < 0.85) {
                current = result.best_verified_x;
            } else {
                current = start_x;
            }
            const double restart_step = trial_options.step * (0.25 + uniform(trial_rng) * 1.75);
            jitter_seed_planes(current, movable_faces, restart_step, normal, trial_rng);
            current_metrics = evaluate_hunt_planes_with_scratch(
                current, original_x, original, affected_edges, false, scratch);
            step = std::max(min_step, restart_step);
            temperature = trial_options.temperature * (0.75 + uniform(trial_rng));
            stagnant = 0;
        }
    }

    HuntMetrics final_canonical = evaluate_hunt_planes_with_scratch(
        best_search_x, original_x, original, affected_edges, true, scratch);
    if (better_verified(final_canonical, result.best_verified)) {
        result.best_verified = final_canonical;
        result.best_verified_x = best_search_x;
    }
    if (better_fast(final_canonical, result.best_fast)) {
        result.best_fast = final_canonical;
        result.best_fast_x = best_search_x;
    }

    return result;
}

int batch_hunt_shape(const LocalRepairOptions& options) {
    g_topology = options.topology;
    set_rand_seed(options.seed);

    const std::filesystem::path out_prefix_path(options.out_prefix);
    if (out_prefix_path.has_parent_path()) {
        std::filesystem::create_directories(out_prefix_path.parent_path());
    }
    const std::filesystem::path report_path(options.report_path);
    if (report_path.has_parent_path()) {
        std::filesystem::create_directories(report_path.parent_path());
    }
    if (!options.stop_file_path.empty()) {
        std::error_code remove_error;
        std::filesystem::remove(options.stop_file_path, remove_error);
    }

    Verts3D original;
    import_obj(options.obj_path.c_str(), original, g_polys);
    if (original.empty() || g_polys.empty()) {
        std::cerr << "Failed to load OBJ: " << options.obj_path << std::endl;
        return 2;
    }
    make_edges(g_polys, g_edges);
    Edges dual_edges;
    dual_graph(g_polys, g_tris, dual_edges);

    Planes original_planes;
    VectorXd original_x;
    if (!options.start_planes_path.empty()) {
        if (!load_plane_state(options.start_planes_path, original_x) ||
            !valid_plane_state(original_x)) {
            std::cerr << "Failed to load valid plane state: "
                      << options.start_planes_path << std::endl;
            return 2;
        }
        x_to_planes(original_x, original_planes);
        planes_to_v3ds(g_tris, original_planes, original);
        if (!is_finite(original)) {
            std::cerr << "Plane state reconstructs non-finite vertices: "
                      << options.start_planes_path << std::endl;
            return 2;
        }
    } else {
        v3ds_to_planes(original, g_polys, original_planes);
        planes_to_x(original_planes, original_x);
    }

    std::vector<int> local_faces = collect_affected_faces(options.movable_vertices);
    local_faces.push_back(4);
    local_faces.push_back(5);
    local_faces = unique_faces(local_faces);
    std::vector<int> target_faces = unique_faces({4, 5});
    std::vector<int> global_faces = all_face_indices();
    const Edges affected_edges = collect_affected_edges(options.movable_vertices);

    HuntMetrics global_best = evaluate_hunt_planes(
        original_x, original_x, original, affected_edges, true);
    HuntMetrics global_fast = global_best;
    VectorXd global_best_x = original_x;
    VectorXd global_fast_x = original_x;

    std::ofstream report(options.report_path, std::ios::app);
    report << "# C++ batch hunt run\n\n";
    report << "- OBJ: `" << options.obj_path << "`\n";
    if (!options.start_planes_path.empty()) {
        report << "- Start planes: `" << options.start_planes_path << "`\n";
    }
    report << "- Base seed: `" << options.seed << "`\n";
    report << "- Topology: `" << options.topology << "`\n";
    report << "- Trials: `" << options.trials << "`\n";
    report << "- Iterations per trial: `" << options.iterations << "`\n";
    report << "- Time limit seconds: `" << options.time_limit_seconds << "`\n";
    report << "- Initial C/I: `" << global_best.crossings << "/" << global_best.intersections << "`\n\n";

    std::cout << "===================" << std::endl;
    std::cout << "Mode      : batch-hunt" << std::endl;
    std::cout << "Space     : face planes, mixed local/global faces" << std::endl;
    std::cout << "Loaded    : " << options.obj_path << std::endl;
    if (!options.start_planes_path.empty()) {
        std::cout << "Planes    : " << options.start_planes_path << std::endl;
    }
    std::cout << "Base seed : " << options.seed << std::endl;
    std::cout << "Topology  : " << g_topology << std::endl;
    std::cout << "Trials    : " << options.trials << std::endl;
    std::cout << "Iters/trial: " << options.iterations << std::endl;
    if (options.time_limit_seconds > 0) {
        std::cout << "Time limit: " << options.time_limit_seconds << " sec" << std::endl;
    }
    std::cout << "Initial C/I: " << global_best.crossings << "/" << global_best.intersections << std::endl;
    std::cout << "===================" << std::endl;

    const unsigned int hardware_threads = std::thread::hardware_concurrency();
    const int worker_count = std::clamp(
        options.threads > 0 ? options.threads : static_cast<int>(hardware_threads == 0 ? 1 : hardware_threads),
        1,
        std::max(1, options.trials)
    );

    report << "- Threads: `" << worker_count << "`\n\n";
    std::cout << "Threads   : " << worker_count << std::endl;
    std::cout << "===================" << std::endl;

    auto make_trial_options = [&](int trial, std::vector<int>& movable_faces, const char*& strategy_name) {
        const int strategy = trial % 8;
        if (strategy < 3) {
            movable_faces = target_faces;
            strategy_name = "target";
        } else if (strategy < 6) {
            movable_faces = local_faces;
            strategy_name = "local";
        } else {
            movable_faces = global_faces;
            strategy_name = "global";
        }

        RNG param_rng(static_cast<RNG::result_type>(
            static_cast<unsigned int>(options.seed) +
            static_cast<unsigned int>(trial * 104729 + strategy * 7919)
        ));
        std::uniform_real_distribution<double> uniform(0.0, 1.0);
        std::uniform_int_distribution<int> seed_noise(0, 999999);

        LocalRepairOptions trial_options = options;
        const long long seed_value =
            static_cast<long long>(options.seed) +
            1000003LL * trial +
            7919LL * strategy +
            seed_noise(param_rng);
        trial_options.seed = static_cast<int>(
            1 + static_cast<unsigned long long>(seed_value) % 2147483646ULL);
        trial_options.step = std::max(0.001, options.step * random_log_scale(param_rng, -1.2, 0.9));
        if (strategy >= 6) {
            trial_options.step *= 0.35;
        } else if (strategy < 3) {
            trial_options.step *= 1.4;
        }
        trial_options.beta = std::clamp(
            options.beta + (uniform(param_rng) - 0.5) * 0.0015,
            0.995,
            0.99999);
        trial_options.temperature = std::max(
            1e-5,
            options.temperature * random_log_scale(param_rng, -0.4, 0.7));
        trial_options.jump_chance = 0.02 + uniform(param_rng) * 0.22;
        trial_options.stagnation = std::max(128, options.stagnation / 4 + static_cast<int>(uniform(param_rng) * options.stagnation));
        trial_options.restarts = std::max(1, options.restarts / 4);
        return trial_options;
    };

    std::atomic<int> next_trial{1};
    std::atomic<int> completed_trials{0};
    std::atomic<int> active_workers{worker_count};
    std::atomic<bool> found{false};
    std::atomic<bool> stop_requested{false};
    std::mutex result_mutex;
    std::vector<std::thread> workers;
    workers.reserve(worker_count);
    const auto started_at = std::chrono::steady_clock::now();
    const auto deadline = options.time_limit_seconds > 0
        ? started_at + std::chrono::seconds(options.time_limit_seconds)
        : std::chrono::steady_clock::time_point::max();

    auto worker = [&]() {
        while (!found.load(std::memory_order_relaxed) &&
               !stop_requested.load(std::memory_order_relaxed)) {
            const int trial = next_trial.fetch_add(1, std::memory_order_relaxed);
            if (trial > options.trials) {
                break;
            }

            std::vector<int> movable_faces;
            const char* strategy_name = "local";
            LocalRepairOptions trial_options = make_trial_options(trial, movable_faces, strategy_name);

            VectorXd trial_start_x;
            {
                std::lock_guard<std::mutex> lock(result_mutex);
                const int source = trial % 10;
                if (trial > worker_count && source < 6) {
                    trial_start_x = global_best_x;
                } else if (trial > worker_count && source < 9) {
                    trial_start_x = global_fast_x;
                } else {
                    trial_start_x = original_x;
                }
            }

            BatchTrialResult result = run_plane_batch_trial(
                trial_options,
                original,
                original_x,
                trial_start_x,
                affected_edges,
                movable_faces,
                trial,
                stop_requested);

            PlaneEvaluationScratch canonical_scratch;
            result.best_verified = evaluate_hunt_planes_with_scratch(
                result.best_verified_x,
                original_x,
                original,
                affected_edges,
                true,
                canonical_scratch);
            result.best_fast = evaluate_hunt_planes_with_scratch(
                result.best_fast_x,
                original_x,
                original,
                affected_edges,
                true,
                canonical_scratch);

            bool saved = false;
            std::lock_guard<std::mutex> lock(result_mutex);
            if (better_verified(result.best_verified, global_best)) {
                global_best = result.best_verified;
                global_best_x = result.best_verified_x;
                const std::string prefix = options.out_prefix + "_best_t" + std::to_string(trial);
                save_plane_candidate_quiet(prefix, global_best_x, result.seed, global_best);
                saved = true;
                std::cout << "GLOBAL BEST trial " << trial << ": C/I "
                          << global_best.crossings << "/" << global_best.intersections
                          << ", target " << global_best.target_crossings
                          << ", strategy " << strategy_name
                          << ", seed " << result.seed << std::endl;
                report << "- GLOBAL BEST trial `" << trial << "`: C/I `"
                       << global_best.crossings << "/" << global_best.intersections
                       << "`, target `" << global_best.target_crossings
                       << "`, strategy `" << strategy_name
                       << "`, seed `" << result.seed << "`\n";
            }
            if (result.best_fast.verified && better_fast(result.best_fast, global_fast)) {
                global_fast = result.best_fast;
                global_fast_x = result.best_fast_x;
                const std::string prefix = options.out_prefix + "_target_t" + std::to_string(trial);
                save_plane_candidate_quiet(prefix, global_fast_x, result.seed, global_fast);
                saved = true;
                std::cout << "TARGET BEST trial " << trial << ": target "
                          << global_fast.target_crossings << ", verified C/I "
                          << global_fast.crossings << "/" << global_fast.intersections
                          << ", strategy " << strategy_name
                          << ", seed " << result.seed << std::endl;
            }

            const bool show_trial =
                saved || (options.report_every > 0 && trial % options.report_every == 0);
            if (show_trial) {
                std::cout << "Trial " << trial << "/" << options.trials
                          << " [" << strategy_name << "] seed " << result.seed
                          << " -> best " << result.best_verified.crossings << "/" << result.best_verified.intersections
                          << ", target " << result.best_verified.target_crossings
                          << " | global " << global_best.crossings << "/" << global_best.intersections
                          << (saved ? " saved" : "")
                          << std::endl;
            }

            if (global_best.crossings == 0 && global_best.intersections == 0) {
                int saved_crossings = 0;
                int saved_intersections = 0;
                if (export_and_validate_found_candidate(
                        "runtime/candidates/FOUND_top4_candidate.obj",
                        global_best_x,
                        saved_crossings,
                        saved_intersections)) {
                    std::cout << "FOUND candidate saved and double-double checked at eps 1e-7/1e-9/1e-11/1e-13: "
                              << "runtime/candidates/FOUND_top4_candidate.obj" << std::endl;
                    report << "\nFOUND candidate saved and double-double checked at eps `1e-7/1e-9/1e-11/1e-13`: "
                           << "`runtime/candidates/FOUND_top4_candidate.obj`\n";
                    found.store(true, std::memory_order_relaxed);
                    stop_requested.store(true, std::memory_order_relaxed);
                } else {
                    std::cout << "Rejected fragile 0/0 after OBJ recheck: C/I "
                              << saved_crossings << "/" << saved_intersections << std::endl;
                }
            }

            const int completed = completed_trials.fetch_add(1, std::memory_order_relaxed) + 1;
            if (options.report_every > 0 && completed % options.report_every == 0) {
                const std::string prefix = options.out_prefix + "_checkpoint";
                save_plane_candidate_quiet(prefix, global_best_x, options.seed + completed, global_best);
                report << "- Checkpoint completed `" << completed << "`: global C/I `"
                       << global_best.crossings << "/" << global_best.intersections << "`\n";
            }
        }
        active_workers.fetch_sub(1, std::memory_order_relaxed);
    };

    for (int i = 0; i < worker_count; ++i) {
        workers.emplace_back(worker);
    }
    while (active_workers.load(std::memory_order_relaxed) > 0) {
        std::this_thread::sleep_for(std::chrono::seconds(1));
        const auto now = std::chrono::steady_clock::now();
        if (now >= deadline) {
            stop_requested.store(true, std::memory_order_relaxed);
        }
        if (!options.stop_file_path.empty() && std::filesystem::exists(options.stop_file_path)) {
            stop_requested.store(true, std::memory_order_relaxed);
        }
        const double elapsed = std::max(
            0.001,
            std::chrono::duration<double>(now - started_at).count());
        const int completed = completed_trials.load(std::memory_order_relaxed);
        const int issued = std::min(
            options.trials,
            next_trial.load(std::memory_order_relaxed) - 1);
        std::lock_guard<std::mutex> lock(result_mutex);
        std::cout << "Progress  : " << completed << "/" << options.trials
                  << " complete, " << issued - completed << " active, "
                  << std::fixed << std::setprecision(2)
                  << (static_cast<double>(completed) / elapsed) << " trials/s, best "
                  << global_best.crossings << "/" << global_best.intersections
                  << std::defaultfloat << std::setprecision(6)
                  << (stop_requested.load(std::memory_order_relaxed) && !found.load(std::memory_order_relaxed)
                      ? " (stopping)"
                      : "")
                  << std::endl;
    }
    for (std::thread& thread : workers) {
        thread.join();
    }
    if (!options.stop_file_path.empty()) {
        std::error_code remove_error;
        std::filesystem::remove(options.stop_file_path, remove_error);
    }

    PlaneEvaluationScratch final_scratch;
    global_best = evaluate_hunt_planes_with_scratch(
        global_best_x, original_x, original, affected_edges, true, final_scratch);
    global_fast = evaluate_hunt_planes_with_scratch(
        global_fast_x, original_x, original, affected_edges, true, final_scratch);

    save_plane_candidate_quiet(options.out_prefix + "_final_target", global_fast_x, options.seed, global_fast);
    save_plane_candidate_quiet(options.out_prefix + "_final_best", global_best_x, options.seed, global_best);

    report << "\n- Final global C/I: `" << global_best.crossings << "/" << global_best.intersections << "`\n";
    report << "- Final target-best C/I: `" << global_fast.crossings << "/" << global_fast.intersections << "`\n";
    report << "- Final target-best target crossings: `" << global_fast.target_crossings << "`\n";
    report << "- Completed trials: `" << completed_trials.load(std::memory_order_relaxed) << "`\n";
    report << "- Stopped early (time limit or GUI): `"
           << (stop_requested.load(std::memory_order_relaxed) && !found.load(std::memory_order_relaxed)
               ? "yes"
               : "no")
           << "`\n";
    report << "- Output prefix: `" << options.out_prefix << "`\n\n";

    std::cout << "===================" << std::endl;
    std::cout << "Batch hunt finished." << std::endl;
    std::cout << "Best verified C/I: " << global_best.crossings << "/" << global_best.intersections << std::endl;
    std::cout << "Best target candidate C/I: " << global_fast.crossings << "/" << global_fast.intersections
              << ", target " << global_fast.target_crossings << std::endl;
    std::cout << "===================" << std::endl;
    return 0;
}

struct GlobalMetrics {
    int crossings = std::numeric_limits<int>::max() / 4;
    int intersections = std::numeric_limits<int>::max() / 4;
    double crossing_loss = std::numeric_limits<double>::infinity();
    double geometry_penalty = std::numeric_limits<double>::infinity();
    double degeneracy_penalty = std::numeric_limits<double>::infinity();
    double energy = std::numeric_limits<double>::infinity();
    // Geometry descriptors are kept separately from the objective.  They are
    // used by the quality-diversity archive; changing archive bins therefore
    // cannot silently change which candidate is considered globally best.
    double min_plane_determinant = 0.0;
    double relative_min_edge = 0.0;
    double min_turn_sine = 0.0;
    double max_vertex_norm = std::numeric_limits<double>::infinity();
    bool canonical = false;
    bool precise = false;
};

double g_global_degeneracy_weight = 0.01;
std::string g_search_device_id = "CPU";

std::uint64_t topology_fingerprint() {
    std::uint64_t hash = 1469598103934665603ULL;
    auto add = [&](int value) {
        hash ^= static_cast<std::uint64_t>(static_cast<std::uint32_t>(value));
        hash *= 1099511628211ULL;
    };
    for (const Face& triangle : g_tris) {
        add(static_cast<int>(triangle.size()));
        for (int value : triangle) {
            add(value);
        }
    }
    for (const Face& polygon : g_polys) {
        add(static_cast<int>(polygon.size()));
        for (int value : polygon) {
            add(value);
        }
    }
    for (const Edge& edge : g_edges) {
        add(edge.first);
        add(edge.second);
    }
    return hash;
}

struct GlobalTopologyState {
    int topology = 0;
    std::uint64_t visits = 0;
    std::uint64_t trials = 0;
    std::uint64_t iterations = 0;
    std::uint64_t run_visits = 0;
    std::uint64_t run_trials = 0;
    std::uint64_t run_iterations = 0;
    std::uint64_t checkpoint_sequence = 0;
    std::uint64_t last_checkpointed_run_visits = 0;
    std::chrono::steady_clock::time_point last_checkpoint_at{};
    bool pending_checkpoint_improvement = false;
    bool has_state = false;
    GlobalMetrics best;
    VectorXd best_x;

    struct DiverseElite {
        std::uint64_t descriptor = 0;
        GlobalMetrics metrics;
        VectorXd x;
        std::string source_run_id;
        std::uint64_t source_seed = 0;
        std::uint64_t source_sequence = 0;
        std::uint64_t selections = 0;
    };
    std::map<std::uint64_t, DiverseElite> archive;
    std::unordered_set<std::uint64_t> archive_dirty;

    struct DiagonalCemState {
        bool initialized = false;
        std::array<double, GLOBAL_PLANE_VALUE_COUNT> mean{};
        std::array<double, GLOBAL_PLANE_VALUE_COUNT> variance{};
        std::uint64_t updates = 0;
    } cem;
    std::unordered_set<std::uint64_t> cem_seen_state_hashes;

    // The scheduler state is intentionally run-local.  Durable geometry and
    // counters are mergeable; a bandit can relearn its allocation cheaply and
    // must not make two independent computers contend for mutable state.
    std::uint64_t scheduler_pulls = 0;
    std::uint64_t scheduler_last_pull = 0;
    std::uint64_t scheduler_last_improvement = 0;
    double scheduler_reward_ema = 0.0;

    std::uint64_t archive_sequence = 0;

    struct RoundTelemetry {
        bool completed = false;
        bool improved = false;
        bool backend_error = false;
        bool spsa_attempted = false;
        bool spsa_accepted = false;
        std::uint64_t evaluated_states = 0;
        std::uint64_t new_trials = 0;
        std::uint64_t verified_candidates = 0;
        std::uint64_t archive_improvements = 0;
        std::uint64_t replica_swaps_attempted = 0;
        std::uint64_t replica_swaps_accepted = 0;
        double kernel_milliseconds = 0.0;
        double transfer_milliseconds = 0.0;
        double wall_milliseconds = 0.0;
        std::array<std::uint64_t, 8> strategy_evaluated{};
        std::array<std::uint64_t, 8> strategy_verified{};
        std::array<std::uint64_t, 8> strategy_archive_improvements{};
        std::array<std::uint64_t, 8> strategy_global_improvements{};
    } last_round;
};

struct GlobalTrialResult {
    GlobalMetrics best;
    VectorXd best_x;
    int iterations = 0;
};

struct GlobalVerifiedCandidate {
    GlobalMetrics metrics;
    VectorXd x;
    std::uint32_t strategy = 0;
    std::uint64_t chain_id = 0;
};

int global_defects(const GlobalMetrics& metrics) {
    return metrics.crossings + metrics.intersections;
}

bool better_global_metrics(const GlobalMetrics& left, const GlobalMetrics& right) {
    if (!std::isfinite(left.energy)) {
        return false;
    }
    if (!std::isfinite(right.energy)) {
        return true;
    }
    return std::make_tuple(
               global_defects(left),
               std::max(left.crossings, left.intersections),
               left.crossings,
               left.energy) <
           std::make_tuple(
               global_defects(right),
               std::max(right.crossings, right.intersections),
               right.crossings,
               right.energy);
}

int logarithmic_archive_bin(
    double value,
    double minimum_log10,
    double maximum_log10,
    int bin_count
) {
    if (!std::isfinite(value) || value <= 0.0 || bin_count <= 1) {
        return 0;
    }
    const double scaled = (std::log10(value) - minimum_log10) /
        (maximum_log10 - minimum_log10);
    return std::clamp(
        static_cast<int>(std::floor(scaled * static_cast<double>(bin_count))),
        0,
        bin_count - 1);
}

std::uint64_t global_archive_descriptor(const GlobalMetrics& metrics) {
    if (!std::isfinite(metrics.energy)) {
        return std::numeric_limits<std::uint64_t>::max();
    }
    const std::uint64_t determinant = static_cast<std::uint64_t>(
        logarithmic_archive_bin(metrics.min_plane_determinant, -10.0, 0.0, 8));
    const std::uint64_t edge = static_cast<std::uint64_t>(
        logarithmic_archive_bin(metrics.relative_min_edge, -8.0, 0.0, 8));
    const std::uint64_t turn = static_cast<std::uint64_t>(
        logarithmic_archive_bin(metrics.min_turn_sine, -8.0, 0.0, 8));
    const std::uint64_t extent = static_cast<std::uint64_t>(
        logarithmic_archive_bin(1.0 + metrics.max_vertex_norm, 0.0, 6.0, 8));
    // C/I and energy are quality, not behavior descriptors.  Consequently a
    // cell can only improve as more computers contribute data.
    return determinant | (edge << 3U) | (turn << 6U) | (extent << 9U);
}

bool add_global_archive_elite(
    GlobalTopologyState& state,
    const GlobalMetrics& metrics,
    const VectorXd& x,
    const std::string& source_run_id,
    std::uint64_t source_seed,
    std::uint64_t source_sequence,
    bool mark_dirty
) {
    constexpr std::size_t kMaximumArchiveCells = 4096;
    if (!std::isfinite(metrics.energy) || x.size() != GLOBAL_PLANE_VALUE_COUNT ||
        !x.allFinite()) {
        return false;
    }
    const std::uint64_t descriptor = global_archive_descriptor(metrics);
    if (descriptor == std::numeric_limits<std::uint64_t>::max()) {
        return false;
    }
    auto existing = state.archive.find(descriptor);
    if (existing != state.archive.end()) {
        if (!better_global_metrics(metrics, existing->second.metrics)) {
            return false;
        }
        const std::uint64_t selections = existing->second.selections;
        existing->second = GlobalTopologyState::DiverseElite{
            descriptor,
            metrics,
            x,
            source_run_id,
            source_seed,
            source_sequence,
            selections};
        if (mark_dirty) {
            state.archive_dirty.insert(descriptor);
        }
        return true;
    }

    state.archive.emplace(
        descriptor,
        GlobalTopologyState::DiverseElite{
            descriptor,
            metrics,
            x,
            source_run_id,
            source_seed,
            source_sequence,
            0});
    if (mark_dirty) {
        state.archive_dirty.insert(descriptor);
    }
    if (state.archive.size() <= kMaximumArchiveCells) {
        return true;
    }

    auto worst = state.archive.begin();
    for (auto it = std::next(state.archive.begin()); it != state.archive.end(); ++it) {
        if (better_global_metrics(worst->second.metrics, it->second.metrics)) {
            worst = it;
        }
    }
    const bool retained = worst->first != descriptor;
    state.archive_dirty.erase(worst->first);
    state.archive.erase(worst);
    return retained;
}

double topology_scheduler_severity(const GlobalMetrics& metrics) {
    return static_cast<double>(global_defects(metrics)) +
        0.10 * static_cast<double>(std::max(metrics.crossings, metrics.intersections)) +
        0.01 * static_cast<double>(metrics.crossings);
}

double topology_scheduler_worstness(
    const GlobalTopologyState& state,
    double best_severity,
    double worst_severity
) {
    if (!state.has_state) {
        return 1.0;
    }
    const double span = worst_severity - best_severity;
    if (span <= 1.0e-9) {
        return 0.0;
    }
    return std::clamp(
        (topology_scheduler_severity(state.best) - best_severity) / span,
        0.0,
        1.0);
}

double topology_bandit_score(
    const GlobalTopologyState& state,
    std::uint64_t total_pulls,
    std::uint64_t current_round,
    double best_severity,
    double worst_severity,
    int best_defects,
    bool prioritize_worst
) {
    const double exploration = 0.60 * std::sqrt(
        std::log(static_cast<double>(total_pulls) + 2.0) /
        (static_cast<double>(state.scheduler_pulls) + 1.0));
    double quality_prior = prioritize_worst ? 0.70 : 0.35;
    if (state.has_state) {
        if (prioritize_worst) {
            quality_prior = 0.65 * topology_scheduler_worstness(
                state,
                best_severity,
                worst_severity);
        } else {
            quality_prior = 0.30 /
                (1.0 + std::max(0, global_defects(state.best) - best_defects));
        }
    }
    const double rounds_since_pull = state.scheduler_last_pull == 0
        ? static_cast<double>(current_round + 1)
        : static_cast<double>(current_round - state.scheduler_last_pull);
    const double staleness = std::min(0.25, rounds_since_pull * 0.01);
    return state.scheduler_reward_ema + exploration + quality_prior + staleness;
}

void update_topology_bandit_reward(
    GlobalTopologyState& state,
    bool had_before,
    const GlobalMetrics& before,
    std::uint64_t archive_improvements,
    std::uint64_t current_round
) {
    double reward = 0.0;
    if (!had_before && state.has_state) {
        reward = 1.0;
    } else if (had_before && state.has_state) {
        const int defect_gain = global_defects(before) - global_defects(state.best);
        reward += static_cast<double>(std::clamp(defect_gain, 0, 4));
        if (defect_gain == 0) {
            const int balance_gain =
                std::max(before.crossings, before.intersections) -
                std::max(state.best.crossings, state.best.intersections);
            reward += 0.50 * static_cast<double>(std::max(0, balance_gain));
            if (balance_gain == 0) {
                reward += 0.25 * static_cast<double>(
                    std::max(0, before.crossings - state.best.crossings));
                if (std::isfinite(before.energy) && before.energy > 0.0 &&
                    state.best.energy < before.energy) {
                    reward += std::min(
                        0.20,
                        (before.energy - state.best.energy) / before.energy);
                }
            }
        }
    }
    reward += std::min(0.10, static_cast<double>(archive_improvements) * 0.01);
    if (state.scheduler_pulls == 0) {
        state.scheduler_reward_ema = reward;
    } else {
        state.scheduler_reward_ema =
            0.90 * state.scheduler_reward_ema + 0.10 * reward;
    }
    state.scheduler_pulls += 1;
    state.scheduler_last_pull = current_round;
    if (reward >= 0.20) {
        state.scheduler_last_improvement = current_round;
    }
}

bool load_global_topology_context(int topology) {
    g_topology = topology;
    open_topology("data/topologies.txt", g_tris, topology);
    if (g_tris.size() != 44) {
        return false;
    }
    dual_graph(g_tris, g_polys, g_edges);
    fix_face_ordering(g_polys, g_edges);
    if (g_polys.size() != 12 || g_edges.size() != 66) {
        return false;
    }
    return std::all_of(g_polys.begin(), g_polys.end(), [](const Face& face) {
        return face.size() == 11;
    });
}

bool normalize_global_plane_state(VectorXd& x) {
    if (!x.allFinite() || x.size() == 0) {
        return false;
    }
    const double rms = x.norm() / std::sqrt(static_cast<double>(x.size() / 3));
    if (!std::isfinite(rms) || rms < 1e-8) {
        return false;
    }
    x /= rms;
    return true;
}

void randomize_global_plane(VectorXd& x, int face_ix, RNG& rng) {
    std::normal_distribution<double> normal(0.0, 1.0);
    Vector3d direction(normal(rng), normal(rng), normal(rng));
    if (direction.squaredNorm() < 1e-12) {
        direction = Vector3d::UnitX();
    }
    direction.normalize();
    const double log_distance = std::clamp(normal(rng) * 0.55, -1.8, 1.8);
    Eigen::Map<Vector3d>(x.data() + face_ix * 3) = direction * std::exp(log_distance);
}

VectorXd make_random_global_state(RNG& rng) {
    VectorXd x(static_cast<int>(g_polys.size() * 3));
    for (int face_ix = 0; face_ix < static_cast<int>(g_polys.size()); ++face_ix) {
        randomize_global_plane(x, face_ix, rng);
    }
    normalize_global_plane_state(x);
    return x;
}

void apply_global_plane_move(
    VectorXd& x,
    double step,
    bool large_jump,
    RNG& rng,
    std::normal_distribution<double>& normal
) {
    const int face_count = static_cast<int>(g_polys.size());
    const double s = step * (large_jump ? 3.5 : 1.0);
    const int move_type = std::uniform_int_distribution<int>(0, 6)(rng);
    const int a = std::uniform_int_distribution<int>(0, face_count - 1)(rng);
    const int b = std::uniform_int_distribution<int>(0, face_count - 1)(rng);
    auto block = [&](int face_ix) {
        return Eigen::Map<Vector3d>(x.data() + face_ix * 3);
    };
    const auto noise = [&]() {
        return Vector3d(normal(rng), normal(rng), normal(rng));
    };

    if (move_type == 0) {
        block(a) += noise() * s;
    } else if (move_type == 1) {
        const Vector3d delta = noise() * s;
        block(a) += delta;
        block(b) -= delta;
    } else if (move_type == 2) {
        block(a) += noise() * s;
        block(b) += noise() * s;
    } else if (move_type == 3) {
        for (int face_ix = 0; face_ix < face_count; ++face_ix) {
            block(face_ix) += noise() * (s * 0.18);
        }
    } else if (move_type == 4) {
        const Vector3d temp = block(a);
        block(a) = block(b);
        block(b) = temp;
    } else if (move_type == 5) {
        randomize_global_plane(x, a, rng);
        if (large_jump) {
            randomize_global_plane(x, b, rng);
        }
    } else {
        const double scale = std::exp(std::clamp(normal(rng) * s * 0.35, -1.2, 1.2));
        block(a) *= scale;
    }
    normalize_global_plane_state(x);
}

GlobalMetrics evaluate_global_state(
    const VectorXd& x,
    bool canonical,
    PlaneEvaluationScratch& scratch
) {
    GlobalMetrics metrics;
    if (!valid_plane_state(x)) {
        return metrics;
    }
    scratch.planes.reserve(g_polys.size());
    scratch.canonical_planes.reserve(g_polys.size());
    scratch.verts.reserve(g_tris.size());
    x_to_planes(x, scratch.planes);
    planes_to_v3ds(g_tris, scratch.planes, scratch.verts);
    if (!is_finite(scratch.verts)) {
        return metrics;
    }

    double min_edge_sq = std::numeric_limits<double>::infinity();
    double max_edge_sq = 0.0;
    double max_vertex_norm = 0.0;
    for (const Vector3d& vertex : scratch.verts) {
        max_vertex_norm = std::max(max_vertex_norm, vertex.norm());
    }
    for (const Edge& edge : g_edges) {
        const double length_sq =
            (scratch.verts[edge.first] - scratch.verts[edge.second]).squaredNorm();
        min_edge_sq = std::min(min_edge_sq, length_sq);
        max_edge_sq = std::max(max_edge_sq, length_sq);
    }
    if (!std::isfinite(min_edge_sq) || !std::isfinite(max_edge_sq) ||
        !std::isfinite(max_vertex_norm) || min_edge_sq < 1e-14 ||
        max_edge_sq > 1e18 || max_vertex_norm > 1e9) {
        return metrics;
    }

    const Planes* metric_planes = &scratch.planes;
    if (canonical) {
        v3ds_to_planes(scratch.verts, g_polys, scratch.canonical_planes);
        metric_planes = &scratch.canonical_planes;
    }
    metrics.crossings = count_self_crossings_strict(
        scratch.verts, *metric_planes, 1e-8, &metrics.crossing_loss);
    metrics.intersections = count_edge_face_intersections_strict(
        scratch.verts, *metric_planes, 1e-8);
    metrics.precise = canonical && promote_precise_counts_if_close(
        scratch.verts, *metric_planes, metrics.crossings, metrics.intersections);

    const double condition = std::sqrt(max_edge_sq / min_edge_sq);
    const double relative_min_edge = std::sqrt(min_edge_sq / max_edge_sq);
    double min_plane_determinant = std::numeric_limits<double>::infinity();
    for (const Face& triangle : g_tris) {
        if (triangle.size() != 3) {
            continue;
        }
        const Vector3d& a = (*metric_planes)[triangle[0]].n;
        const Vector3d& b = (*metric_planes)[triangle[1]].n;
        const Vector3d& c = (*metric_planes)[triangle[2]].n;
        min_plane_determinant = std::min(
            min_plane_determinant,
            std::abs(a.dot(b.cross(c))));
    }
    if (!std::isfinite(min_plane_determinant)) {
        return GlobalMetrics{};
    }
    double min_turn_sine = std::numeric_limits<double>::infinity();
    for (const Face& face : g_polys) {
        for (size_t i = 0; i < face.size(); ++i) {
            const Vector3d a = scratch.verts[face[(i + face.size() - 1) % face.size()]] -
                scratch.verts[face[i]];
            const Vector3d b = scratch.verts[face[(i + 1) % face.size()]] -
                scratch.verts[face[i]];
            const double denominator = a.norm() * b.norm();
            if (denominator > 1e-15) {
                min_turn_sine = std::min(
                    min_turn_sine,
                    a.cross(b).norm() / denominator);
            }
        }
    }
    if (!std::isfinite(min_turn_sine)) {
        return GlobalMetrics{};
    }
    const double determinant_barrier =
        std::log1p(0.02 / std::max(1e-10, min_plane_determinant));
    const double edge_barrier =
        std::log1p(0.002 / std::max(1e-10, relative_min_edge));
    const double turn_barrier =
        std::log1p(0.002 / std::max(1e-10, min_turn_sine));
    const double extent_barrier = 0.10 * std::log1p(max_vertex_norm / 100.0);
    const double worst_degeneracy_barrier = std::max({
        determinant_barrier,
        edge_barrier,
        turn_barrier,
        extent_barrier});
    const double secondary_degeneracy_barriers =
        determinant_barrier + edge_barrier + turn_barrier + extent_barrier -
        worst_degeneracy_barrier;
    metrics.min_plane_determinant = min_plane_determinant;
    metrics.relative_min_edge = relative_min_edge;
    metrics.min_turn_sine = min_turn_sine;
    metrics.max_vertex_norm = max_vertex_norm;
    metrics.degeneracy_penalty = g_global_degeneracy_weight *
        (worst_degeneracy_barrier + 0.05 * secondary_degeneracy_barriers);
    metrics.geometry_penalty =
        0.0010 * std::min(20.0, std::log1p(condition)) +
        0.0002 * std::min(20.0, std::log1p(max_vertex_norm)) +
        metrics.degeneracy_penalty;
    metrics.energy =
        0.05 * static_cast<double>(global_defects(metrics)) +
        0.01 * static_cast<double>(metrics.crossings) +
        0.02 * std::min(10.0, metrics.crossing_loss) +
        metrics.geometry_penalty;
    metrics.canonical = canonical;
    return metrics;
}

bool round_trip_global_state_to_fp32(
    VectorXd& x,
    GlobalMetrics& metrics
) {
    if (x.size() != GLOBAL_PLANE_VALUE_COUNT) {
        return false;
    }
    VectorXd exact(GLOBAL_PLANE_VALUE_COUNT);
    for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
        const float value = static_cast<float>(x[component]);
        if (!std::isfinite(value)) {
            return false;
        }
        exact[component] = static_cast<double>(value);
    }
    PlaneEvaluationScratch scratch;
    GlobalMetrics exact_metrics = evaluate_global_state(exact, true, scratch);
    if (!std::isfinite(exact_metrics.energy)) {
        return false;
    }
    x = std::move(exact);
    metrics = exact_metrics;
    return true;
}

void update_diagonal_cem(
    GlobalTopologyState& state,
    const std::vector<GlobalVerifiedCandidate>& candidates
) {
    if (candidates.empty()) {
        return;
    }
    std::vector<const GlobalVerifiedCandidate*> ranked;
    ranked.reserve(candidates.size());
    const bool has_injected_candidates = std::any_of(
        candidates.begin(),
        candidates.end(),
        [](const GlobalVerifiedCandidate& candidate) {
            return candidate.strategy ==
                static_cast<std::uint32_t>(cuda_search::StrategyKind::Injected);
        });
    if (state.cem_seen_state_hashes.size() > 16384) {
        state.cem_seen_state_hashes.clear();
    }
    for (const GlobalVerifiedCandidate& candidate : candidates) {
        if (has_injected_candidates &&
            candidate.strategy !=
                static_cast<std::uint32_t>(cuda_search::StrategyKind::Injected)) {
            continue;
        }
        if (candidate.x.size() == GLOBAL_PLANE_VALUE_COUNT &&
            std::isfinite(candidate.metrics.energy)) {
            std::uint64_t hash = 1469598103934665603ULL;
            for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
                const float value = static_cast<float>(candidate.x[component]);
                std::uint32_t bits = 0;
                std::memcpy(&bits, &value, sizeof(bits));
                hash ^= static_cast<std::uint64_t>(bits);
                hash *= 1099511628211ULL;
            }
            if (!state.cem_seen_state_hashes.insert(hash).second) {
                continue;
            }
            ranked.push_back(&candidate);
        }
    }
    std::stable_sort(ranked.begin(), ranked.end(), [](const auto* left, const auto* right) {
        return better_global_metrics(left->metrics, right->metrics);
    });
    if (ranked.empty()) {
        return;
    }
    const VectorXd& mode_anchor = ranked.front()->x;
    const int anchor_defects = global_defects(ranked.front()->metrics);
    ranked.erase(
        std::remove_if(
            std::next(ranked.begin()),
            ranked.end(),
            [&](const GlobalVerifiedCandidate* candidate) {
                return global_defects(candidate->metrics) > anchor_defects + 2 ||
                    (candidate->x - mode_anchor).norm() > 1.5;
            }),
        ranked.end());

    // CEM supplies the elite mean; the diagonal covariance update also keeps
    // the displacement of the mean (the useful part of diagonal CMA) so the
    // distribution does not collapse after one unusually tight shortlist.
    const std::size_t elite_count = std::min<std::size_t>(
        ranked.size(),
        std::max<std::size_t>(4, ranked.size() / 6));
    std::array<double, GLOBAL_PLANE_VALUE_COUNT> elite_mean{};
    std::array<double, GLOBAL_PLANE_VALUE_COUNT> elite_variance{};
    double weight_sum = 0.0;
    for (std::size_t rank = 0; rank < elite_count; ++rank) {
        const double weight = std::log(static_cast<double>(elite_count) + 1.5) -
            std::log(static_cast<double>(rank) + 1.0);
        weight_sum += weight;
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            elite_mean[component] += weight * ranked[rank]->x[component];
        }
    }
    for (double& value : elite_mean) {
        value /= weight_sum;
    }
    for (std::size_t rank = 0; rank < elite_count; ++rank) {
        const double weight = std::log(static_cast<double>(elite_count) + 1.5) -
            std::log(static_cast<double>(rank) + 1.0);
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            const double delta = ranked[rank]->x[component] - elite_mean[component];
            elite_variance[component] += weight * delta * delta;
        }
    }
    constexpr double kMinimumVariance = 0.0025 * 0.0025;
    constexpr double kMaximumVariance = 0.35 * 0.35;
    for (double& value : elite_variance) {
        value = std::clamp(value / weight_sum, kMinimumVariance, kMaximumVariance);
    }

    if (!state.cem.initialized) {
        state.cem.mean = elite_mean;
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            state.cem.variance[component] = std::max(
                elite_variance[component],
                0.035 * 0.035);
        }
        state.cem.initialized = true;
    } else {
        constexpr double kLearningRate = 0.18;
        constexpr double kEvolutionWeight = 0.10;
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            const double displacement = elite_mean[component] - state.cem.mean[component];
            state.cem.mean[component] += kLearningRate * displacement;
            const double target_variance = elite_variance[component] +
                kEvolutionWeight * displacement * displacement;
            state.cem.variance[component] = std::clamp(
                (1.0 - kLearningRate) * state.cem.variance[component] +
                    kLearningRate * target_variance,
                kMinimumVariance,
                kMaximumVariance);
        }
    }
    state.cem.updates += 1;
}

std::vector<cuda_search::PlaneState> make_hybrid_seed_pool(
    GlobalTopologyState& state,
    std::uint64_t seed,
    int maximum_states
) {
    std::vector<cuda_search::PlaneState> result;
    if (maximum_states <= 0) {
        return result;
    }
    result.reserve(static_cast<std::size_t>(maximum_states));
    RNG rng(static_cast<RNG::result_type>(seed));
    std::normal_distribution<double> normal(0.0, 1.0);

    const int archive_target = std::min<int>(
        maximum_states / 2,
        static_cast<int>(state.archive.size()));
    std::vector<GlobalTopologyState::DiverseElite*> archive_candidates;
    archive_candidates.reserve(state.archive.size());
    for (auto& item : state.archive) {
        archive_candidates.push_back(&item.second);
    }
    // Sampling the least-used cells first is the MAP-Elites coverage pressure:
    // high-quality cells do not monopolize all descendants.
    std::shuffle(archive_candidates.begin(), archive_candidates.end(), rng);
    std::stable_sort(
        archive_candidates.begin(),
        archive_candidates.end(),
        [](const auto* left, const auto* right) {
            return left->selections < right->selections;
        });
    for (int index = 0; index < archive_target; ++index) {
        GlobalTopologyState::DiverseElite& elite = *archive_candidates[index];
        cuda_search::PlaneState state_value;
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            state_value.values[static_cast<std::size_t>(component)] =
                static_cast<float>(elite.x[component]);
        }
        elite.selections += 1;
        result.push_back(state_value);
    }

    while (state.cem.initialized &&
           static_cast<int>(result.size()) < maximum_states) {
        VectorXd sample(GLOBAL_PLANE_VALUE_COUNT);
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            sample[component] = state.cem.mean[component] +
                std::sqrt(state.cem.variance[component]) * normal(rng);
        }
        if (!normalize_global_plane_state(sample)) {
            continue;
        }
        cuda_search::PlaneState state_value;
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            state_value.values[static_cast<std::size_t>(component)] =
                static_cast<float>(sample[component]);
        }
        result.push_back(state_value);
    }
    return result;
}

double smooth_spsa_objective(const GlobalMetrics& metrics) {
    if (!std::isfinite(metrics.energy)) {
        return 1.0e30;
    }
    // Counts remain a weak guide when a perturbation crosses a discrete
    // boundary; clearance and conditioning provide the differentiable signal.
    return 0.004 * static_cast<double>(global_defects(metrics)) +
        0.035 * std::min(20.0, metrics.crossing_loss) +
        metrics.geometry_penalty;
}

bool refine_with_spsa_adam(
    const VectorXd& start,
    std::uint64_t seed,
    int iterations,
    VectorXd& refined,
    GlobalMetrics& refined_metrics
) {
    if (start.size() != GLOBAL_PLANE_VALUE_COUNT || iterations <= 0) {
        return false;
    }
    RNG rng(static_cast<RNG::result_type>(seed));
    std::uniform_int_distribution<int> sign(0, 1);
    PlaneEvaluationScratch scratch;
    VectorXd current = start;
    std::array<double, GLOBAL_PLANE_VALUE_COUNT> first_moment{};
    std::array<double, GLOBAL_PLANE_VALUE_COUNT> second_moment{};
    constexpr double kBeta1 = 0.82;
    constexpr double kBeta2 = 0.97;

    for (int iteration = 0; iteration < iterations; ++iteration) {
        const double perturbation = 0.022 /
            std::pow(static_cast<double>(iteration + 1), 0.101);
        const double learning_rate = 0.012 /
            std::pow(static_cast<double>(iteration + 5), 0.602);
        VectorXd delta(GLOBAL_PLANE_VALUE_COUNT);
        VectorXd plus = current;
        VectorXd minus = current;
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            delta[component] = sign(rng) == 0 ? -1.0 : 1.0;
            plus[component] += perturbation * delta[component];
            minus[component] -= perturbation * delta[component];
        }
        if (!normalize_global_plane_state(plus) ||
            !normalize_global_plane_state(minus)) {
            continue;
        }
        const double plus_value = smooth_spsa_objective(
            evaluate_global_state(plus, false, scratch));
        const double minus_value = smooth_spsa_objective(
            evaluate_global_state(minus, false, scratch));
        if (!std::isfinite(plus_value) || !std::isfinite(minus_value)) {
            continue;
        }
        const double directional = std::clamp(
            (plus_value - minus_value) / (2.0 * perturbation),
            -100.0,
            100.0);
        const double beta1_power = std::pow(kBeta1, iteration + 1);
        const double beta2_power = std::pow(kBeta2, iteration + 1);
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            const double gradient = directional * delta[component];
            first_moment[component] =
                kBeta1 * first_moment[component] + (1.0 - kBeta1) * gradient;
            second_moment[component] =
                kBeta2 * second_moment[component] + (1.0 - kBeta2) * gradient * gradient;
            const double corrected_first = first_moment[component] / (1.0 - beta1_power);
            const double corrected_second = second_moment[component] / (1.0 - beta2_power);
            current[component] -= learning_rate * corrected_first /
                (std::sqrt(corrected_second) + 1.0e-8);
        }
        normalize_global_plane_state(current);
    }

    refined_metrics = evaluate_global_state(current, true, scratch);
    if (!std::isfinite(refined_metrics.energy)) {
        return false;
    }
    refined = std::move(current);
    return true;
}

GlobalTrialResult run_global_trial(
    const LocalRepairOptions& options,
    const VectorXd& start_x,
    bool start_fresh,
    int seed,
    int iterations,
    const std::atomic<bool>& stop_requested
) {
    GlobalTrialResult result;
    RNG rng(static_cast<RNG::result_type>(seed));
    std::normal_distribution<double> normal(0.0, 1.0);
    std::uniform_real_distribution<double> uniform(0.0, 1.0);
    PlaneEvaluationScratch scratch;

    VectorXd current = start_fresh || start_x.size() == 0
        ? make_random_global_state(rng)
        : start_x;
    if (!start_fresh && start_x.size() != 0) {
        for (int i = 0; i < 3; ++i) {
            apply_global_plane_move(current, options.step * 0.35, false, rng, normal);
        }
    }

    GlobalMetrics current_metrics = evaluate_global_state(current, false, scratch);
    for (int retry = 0; !std::isfinite(current_metrics.energy) && retry < 8; ++retry) {
        current = make_random_global_state(rng);
        current_metrics = evaluate_global_state(current, false, scratch);
    }
    GlobalMetrics best_search = current_metrics;
    VectorXd best_search_x = current;
    result.best = evaluate_global_state(current, true, scratch);
    result.best_x = current;

    RNG parameter_rng(static_cast<RNG::result_type>(seed ^ 0x5bd1e995));
    const double base_step = std::max(
        0.002,
        options.step * random_log_scale(parameter_rng, -1.0, 0.65));
    const double base_temperature = std::max(
        1e-5,
        options.temperature * random_log_scale(parameter_rng, -0.35, 0.65));
    const double beta = std::clamp(
        options.beta + (uniform(parameter_rng) - 0.5) * 0.0012,
        0.995,
        0.99999);
    const double min_step = std::max(1e-9, base_step * options.min_step_ratio);
    double step = base_step;
    double temperature = base_temperature;
    int stagnant = 0;
    int restarts = 0;

    for (int iter = 1; iter <= iterations; ++iter) {
        if ((iter & 255) == 0 && stop_requested.load(std::memory_order_relaxed)) {
            break;
        }
        VectorXd candidate = current;
        const bool large_jump = uniform(rng) < options.jump_chance;
        apply_global_plane_move(candidate, step, large_jump, rng, normal);
        GlobalMetrics candidate_metrics = evaluate_global_state(candidate, false, scratch);

        const bool improves_current = candidate_metrics.energy < current_metrics.energy;
        const double exponent = std::clamp(
            (current_metrics.energy - candidate_metrics.energy) /
                std::max(1e-6, temperature),
            -80.0,
            0.0);
        if (improves_current || uniform(rng) < std::exp(exponent)) {
            current = candidate;
            current_metrics = candidate_metrics;
        }

        bool improved = false;
        if (better_global_metrics(candidate_metrics, best_search)) {
            best_search = candidate_metrics;
            best_search_x = candidate;
            GlobalMetrics canonical = evaluate_global_state(candidate, true, scratch);
            if (better_global_metrics(canonical, result.best)) {
                result.best = canonical;
                result.best_x = candidate;
            }
            improved = true;
        }

        result.iterations = iter;
        stagnant = improved ? 0 : stagnant + 1;
        step = std::max(min_step, step * beta);
        temperature = std::max(base_temperature * 0.04, temperature * beta);
        if (result.best.canonical && result.best.crossings == 0 && result.best.intersections == 0) {
            break;
        }

        if (stagnant >= std::max(256, options.stagnation) &&
            restarts < std::max(4, options.restarts / 8)) {
            restarts += 1;
            if (uniform(rng) < 0.72 && best_search_x.size() != 0) {
                current = best_search_x;
                for (int i = 0; i < 4; ++i) {
                    apply_global_plane_move(current, base_step, true, rng, normal);
                }
            } else {
                current = make_random_global_state(rng);
            }
            current_metrics = evaluate_global_state(current, false, scratch);
            step = base_step * (0.5 + uniform(rng) * 1.5);
            temperature = base_temperature * (0.75 + uniform(rng));
            stagnant = 0;
        }
    }

    GlobalMetrics final_canonical = evaluate_global_state(best_search_x, true, scratch);
    if (better_global_metrics(final_canonical, result.best)) {
        result.best = final_canonical;
        result.best_x = best_search_x;
    }
    return result;
}

std::filesystem::path global_topology_dir(
    const LocalRepairOptions& options,
    int topology
) {
    return std::filesystem::path(options.global_dir) /
        ("topology_" + std::to_string(topology));
}

bool load_global_topology_state(
    const LocalRepairOptions& options,
    GlobalTopologyState& state
) {
    const std::filesystem::path dir = global_topology_dir(options, state.topology);
    const std::filesystem::path state_path = dir / "resume.planes";
    if (!std::filesystem::exists(state_path)) {
        return false;
    }
    VectorXd x;
    if (!load_plane_state(state_path, x) || !valid_plane_state(x)) {
        return false;
    }
    PlaneEvaluationScratch scratch;
    GlobalMetrics metrics = evaluate_global_state(x, true, scratch);
    if (!std::isfinite(metrics.energy) ||
        !round_trip_global_state_to_fp32(x, metrics)) {
        return false;
    }
    state.best_x = x;
    state.best = metrics;
    state.has_state = true;
    const std::filesystem::path obj_path = dir / "resume.obj";
    if (!std::filesystem::exists(obj_path)) {
        export_plane_candidate(obj_path.string().c_str(), state.best_x);
    }

    std::ifstream meta(dir / "resume.meta");
    int version = 0;
    if (meta >> version >> state.visits >> state.trials && version != 1) {
        state.visits = 0;
        state.trials = 0;
    }
    return true;
}

void load_mergeable_checkpoints(
    const LocalRepairOptions& options,
    std::vector<GlobalTopologyState>& states
) {
    using namespace szilassi::checkpoint;
    const CheckpointScan scan = scan_checkpoints(options.global_dir);

    // Preserve legacy previews as archive seeds before a newer mergeable best
    // replaces state.best below.  They are inputs only; they never contribute
    // to mergeable run counters.
    for (GlobalTopologyState& state : states) {
        if (state.has_state && load_global_topology_context(state.topology)) {
            add_global_archive_elite(
                state,
                state.best,
                state.best_x,
                "legacy-preview",
                0,
                0,
                false);
        }
    }

    const LatestCheckpointMap latest = select_latest_per_run_topology(scan);
    for (const auto& entry : latest) {
        const GlobalCheckpoint& checkpoint = entry.second.checkpoint;
        if (checkpoint.topology < 0 || checkpoint.topology >= NUM_TOPOLOGIES ||
            checkpoint.plane_coefficients.size() != static_cast<size_t>(GLOBAL_PLANE_VALUE_COUNT)) {
            continue;
        }
        if (!load_global_topology_context(checkpoint.topology)) {
            continue;
        }
        if (checkpoint.topology_fingerprint != 0 &&
            checkpoint.topology_fingerprint != topology_fingerprint()) {
            std::cerr << "Skipping checkpoint with a different topology fingerprint: "
                      << entry.second.path.string() << std::endl;
            continue;
        }
        VectorXd x(GLOBAL_PLANE_VALUE_COUNT);
        for (int i = 0; i < GLOBAL_PLANE_VALUE_COUNT; ++i) {
            x[i] = static_cast<double>(checkpoint.plane_coefficients[static_cast<size_t>(i)]);
        }
        if (!valid_plane_state(x)) {
            continue;
        }
        PlaneEvaluationScratch scratch;
        GlobalMetrics metrics = evaluate_global_state(x, true, scratch);
        if (!std::isfinite(metrics.energy)) {
            continue;
        }
        GlobalTopologyState& state = states[checkpoint.topology];
        state.visits += checkpoint.visits;
        state.trials += checkpoint.completed_trials;
        state.iterations += checkpoint.completed_iterations;
        if (!state.has_state || better_global_metrics(metrics, state.best)) {
            state.best = metrics;
            state.best_x = x;
            state.has_state = true;
        }
    }

    std::unordered_set<std::string> conflicted_paths;
    for (const CheckpointConflict& conflict : scan.conflicts) {
        for (const std::filesystem::path& path : conflict.paths) {
            conflicted_paths.insert(path.lexically_normal().generic_string());
        }
    }
    std::array<std::vector<const CheckpointRecord*>, NUM_TOPOLOGIES> by_topology;
    for (const CheckpointRecord& record : scan.valid) {
        if (record.checkpoint.topology >= 0 &&
            record.checkpoint.topology < NUM_TOPOLOGIES &&
            conflicted_paths.find(record.path.lexically_normal().generic_string()) ==
                conflicted_paths.end()) {
            by_topology[static_cast<std::size_t>(record.checkpoint.topology)].push_back(
                &record);
        }
    }
    std::size_t imported_unique_states = 0;
    for (int topology = 0; topology < NUM_TOPOLOGIES; ++topology) {
        if (by_topology[static_cast<std::size_t>(topology)].empty() ||
            !load_global_topology_context(topology)) {
            continue;
        }
        const std::uint64_t expected_fingerprint = topology_fingerprint();
        std::unordered_set<std::string> seen_states;
        for (const CheckpointRecord* record : by_topology[static_cast<std::size_t>(topology)]) {
            const GlobalCheckpoint& checkpoint = record->checkpoint;
            if (checkpoint.plane_coefficients.size() !=
                    static_cast<std::size_t>(GLOBAL_PLANE_VALUE_COUNT) ||
                (checkpoint.topology_fingerprint != 0 &&
                 checkpoint.topology_fingerprint != expected_fingerprint)) {
                continue;
            }
            const char* raw = reinterpret_cast<const char*>(
                checkpoint.plane_coefficients.data());
            std::string exact_state(
                raw,
                raw + sizeof(float) * checkpoint.plane_coefficients.size());
            if (!seen_states.insert(std::move(exact_state)).second) {
                continue;
            }
            VectorXd x(GLOBAL_PLANE_VALUE_COUNT);
            for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
                x[component] = static_cast<double>(
                    checkpoint.plane_coefficients[static_cast<std::size_t>(component)]);
            }
            if (!valid_plane_state(x)) {
                continue;
            }
            PlaneEvaluationScratch scratch;
            GlobalMetrics metrics = evaluate_global_state(x, true, scratch);
            if (!std::isfinite(metrics.energy)) {
                continue;
            }
            imported_unique_states += 1;
            add_global_archive_elite(
                states[topology],
                metrics,
                x,
                checkpoint.run_id,
                checkpoint.base_seed,
                checkpoint.sequence,
                false);
        }
    }

    const szilassi::archive::ArchiveScan archive_scan =
        szilassi::archive::scan_archives(options.global_dir);
    const szilassi::archive::ArchiveGenerationMap archive_generations =
        szilassi::archive::select_nonconflicting_deltas(archive_scan);
    std::size_t imported_archive_cells = 0;
    std::size_t imported_archive_deltas = 0;
    std::size_t incompatible_archive_deltas = 0;
    for (int topology = 0; topology < NUM_TOPOLOGIES; ++topology) {
        if (!load_global_topology_context(topology)) {
            continue;
        }
        const std::uint64_t expected_fingerprint = topology_fingerprint();
        std::unordered_set<std::string> seen_archive_states;
        for (const auto& generation : archive_generations) {
            const szilassi::archive::ArchiveRecord& record = generation.second;
            if (record.delta.topology != topology) {
                continue;
            }
            if (record.delta.objective_version <
                    MIN_COMPATIBLE_ARCHIVE_OBJECTIVE_VERSION ||
                record.delta.objective_version > GLOBAL_OBJECTIVE_VERSION ||
                record.delta.topology_fingerprint != expected_fingerprint) {
                incompatible_archive_deltas += 1;
                continue;
            }
            imported_archive_deltas += 1;
            // Re-evaluate every distinct state before per-cell selection.
            // Energy contains configurable weights, so selecting by serialized
            // energy first could discard the current run's true winner.
            for (const szilassi::archive::ArchiveEntry& entry : record.delta.entries) {
                const char* raw = reinterpret_cast<const char*>(
                    entry.plane_coefficients.data());
                std::string exact_state(
                    raw,
                    raw + sizeof(float) * entry.plane_coefficients.size());
                if (!seen_archive_states.insert(std::move(exact_state)).second) {
                    continue;
                }
                VectorXd x(GLOBAL_PLANE_VALUE_COUNT);
                for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
                    x[component] = static_cast<double>(
                        entry.plane_coefficients[static_cast<std::size_t>(component)]);
                }
                if (!valid_plane_state(x)) {
                    continue;
                }
                PlaneEvaluationScratch scratch;
                const GlobalMetrics metrics = evaluate_global_state(x, true, scratch);
                if (!std::isfinite(metrics.energy)) {
                    continue;
                }
                if (add_global_archive_elite(
                        states[topology],
                        metrics,
                        x,
                        entry.origin_run_id,
                        entry.origin_seed,
                        entry.origin_sequence,
                        false)) {
                    imported_archive_cells += 1;
                }
                if (!states[topology].has_state ||
                    better_global_metrics(metrics, states[topology].best)) {
                    states[topology].best = metrics;
                    states[topology].best_x = x;
                    states[topology].has_state = true;
                }
            }
        }
    }
    for (GlobalTopologyState& state : states) {
        std::vector<GlobalVerifiedCandidate> bootstrap_candidates;
        bootstrap_candidates.reserve(state.archive.size());
        for (const auto& item : state.archive) {
            bootstrap_candidates.push_back(GlobalVerifiedCandidate{
                item.second.metrics,
                item.second.x,
                0,
                0});
        }
        state.cem = {};
        update_diagonal_cem(state, bootstrap_candidates);
    }
    std::cout << "Loaded " << imported_unique_states
              << " unique saved FP32 states into the diversity archive."
              << std::endl;
    if (imported_archive_deltas != 0 || imported_archive_cells != 0) {
        std::cout << "Merged " << imported_archive_deltas
                  << " archive delta(s), improving " << imported_archive_cells
                  << " local MAP-Elites cell(s)." << std::endl;
    }
    if (incompatible_archive_deltas != 0) {
        std::cerr << "Ignored " << incompatible_archive_deltas
                  << " archive delta(s) with a different objective version or topology."
                  << std::endl;
    }
    if (!scan.rejected.empty()) {
        std::cerr << "Ignored " << scan.rejected.size()
                  << " incomplete or corrupt checkpoint(s); older generations remain usable."
                  << std::endl;
    }
    if (!scan.conflicts.empty()) {
        std::cerr << "Ignored " << scan.conflicts.size()
                  << " conflicting checkpoint sequence(s)." << std::endl;
    }
    if (!archive_scan.rejected.empty()) {
        std::cerr << "Ignored " << archive_scan.rejected.size()
                  << " incomplete or corrupt archive delta(s)." << std::endl;
    }
    if (!archive_scan.conflicts.empty()) {
        std::cerr << "Ignored " << archive_scan.conflicts.size()
                  << " conflicting archive delta sequence(s)." << std::endl;
    }
}

bool commit_mergeable_checkpoint(
    const LocalRepairOptions& options,
    const szilassi::checkpoint::RunIdentity& identity,
    GlobalTopologyState& state,
    szilassi::checkpoint::CheckpointReason reason
) {
    using namespace szilassi::checkpoint;
    if (!state.has_state || state.best_x.size() != GLOBAL_PLANE_VALUE_COUNT) {
        return false;
    }
    GlobalCheckpoint checkpoint;
    checkpoint.run_id = identity.run_id;
    checkpoint.node_id = identity.node_id;
    checkpoint.producer_id = "Szilassi CUDA global search";
    checkpoint.device_id = g_search_device_id;
    checkpoint.sequence = ++state.checkpoint_sequence;
    checkpoint.topology = state.topology;
    checkpoint.topology_first = options.topology_from;
    checkpoint.topology_last = options.topology_to;
    checkpoint.objective_version = GLOBAL_OBJECTIVE_VERSION;
    checkpoint.topology_fingerprint = topology_fingerprint();
    checkpoint.base_seed = static_cast<std::uint64_t>(
        static_cast<std::uint32_t>(options.seed));
    checkpoint.next_work_unit = state.run_visits + 1;
    checkpoint.completed_work_units = state.run_visits;
    checkpoint.completed_trials = state.run_trials;
    checkpoint.completed_iterations = state.run_iterations;
    checkpoint.visits = state.run_visits;
    checkpoint.plane_coefficients.resize(GLOBAL_PLANE_VALUE_COUNT);
    VectorXd stored_x(GLOBAL_PLANE_VALUE_COUNT);
    for (int i = 0; i < GLOBAL_PLANE_VALUE_COUNT; ++i) {
        const float value = static_cast<float>(state.best_x[i]);
        checkpoint.plane_coefficients[static_cast<size_t>(i)] = value;
        stored_x[i] = static_cast<double>(value);
    }
    PlaneEvaluationScratch scratch;
    GlobalMetrics stored_metrics = evaluate_global_state(stored_x, true, scratch);
    if (!std::isfinite(stored_metrics.energy)) {
        std::cerr << "FP32 checkpoint round-trip is invalid for topology "
                  << state.topology << std::endl;
        return false;
    }
    checkpoint.crossings = stored_metrics.crossings;
    checkpoint.intersections = stored_metrics.intersections;
    checkpoint.crossing_loss = stored_metrics.crossing_loss;
    checkpoint.degeneracy_penalty = stored_metrics.degeneracy_penalty;
    checkpoint.energy = stored_metrics.energy;
    checkpoint.verification = stored_metrics.precise
        ? VerificationPrecision::DoubleDouble
        : VerificationPrecision::Double;
    checkpoint.reason = reason;
    checkpoint.flags = CheckpointFlagCanonical |
        (stored_metrics.precise ? CheckpointFlagDdVerified : CheckpointFlagNone) |
        ((stored_metrics.crossings == 0 && stored_metrics.intersections == 0)
            ? CheckpointFlagFound
            : CheckpointFlagNone);
    const CommitResult result = commit_checkpoint(options.global_dir, std::move(checkpoint));
    if (!result) {
        std::cerr << "Durable checkpoint failed for topology " << state.topology
                  << ": " << result.error << std::endl;
        return false;
    }
    return true;
}

bool commit_mergeable_archive_delta(
    const LocalRepairOptions& options,
    const szilassi::checkpoint::RunIdentity& identity,
    GlobalTopologyState& state
) {
    if (state.archive_dirty.empty()) {
        return true;
    }
    using namespace szilassi::archive;
    std::map<std::uint64_t, ArchiveEntry> exact_entries;
    std::vector<std::uint64_t> persisted_descriptors;
    std::vector<std::uint64_t> stale_descriptors;
    const std::int64_t created_at = szilassi::checkpoint::unix_time_ns_now();
    for (std::uint64_t dirty_descriptor : state.archive_dirty) {
        const auto found = state.archive.find(dirty_descriptor);
        if (found == state.archive.end()) {
            stale_descriptors.push_back(dirty_descriptor);
            continue;
        }
        const GlobalTopologyState::DiverseElite& elite = found->second;
        ArchiveEntry entry;
        VectorXd exact_x(GLOBAL_PLANE_VALUE_COUNT);
        for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) {
            const float value = static_cast<float>(elite.x[component]);
            entry.plane_coefficients[static_cast<std::size_t>(component)] = value;
            exact_x[component] = static_cast<double>(value);
        }
        PlaneEvaluationScratch scratch;
        const GlobalMetrics metrics = evaluate_global_state(exact_x, true, scratch);
        if (!std::isfinite(metrics.energy)) {
            std::cerr << "Archive FP32 round-trip is invalid for topology "
                      << state.topology << std::endl;
            continue;
        }
        entry.descriptor_key = global_archive_descriptor(metrics);
        entry.quality.crossings = metrics.crossings;
        entry.quality.intersections = metrics.intersections;
        entry.quality.crossing_loss = metrics.crossing_loss;
        entry.quality.geometry_penalty = metrics.geometry_penalty;
        entry.quality.degeneracy_penalty = metrics.degeneracy_penalty;
        entry.quality.energy = metrics.energy;
        entry.quality.min_abs_determinant = metrics.min_plane_determinant;
        entry.quality.min_edge_ratio = metrics.relative_min_edge;
        entry.quality.min_turn_sine = metrics.min_turn_sine;
        entry.quality.extent = metrics.max_vertex_norm;
        entry.verification = metrics.precise
            ? VerificationPrecision::DoubleDouble
            : VerificationPrecision::Double;
        entry.flags = ArchiveEntryFlagCanonical |
            (metrics.precise ? ArchiveEntryFlagDdVerified : ArchiveEntryFlagNone);
        entry.origin_run_id = elite.source_run_id.empty()
            ? identity.run_id
            : elite.source_run_id;
        entry.origin_seed = elite.source_seed != 0
            ? elite.source_seed
            : static_cast<std::uint64_t>(static_cast<std::uint32_t>(options.seed));
        entry.origin_sequence = elite.source_sequence;
        entry.discovered_unix_ns = created_at;
        const auto existing = exact_entries.find(entry.descriptor_key);
        if (existing == exact_entries.end()) {
            exact_entries.emplace(entry.descriptor_key, std::move(entry));
        } else if (quality_is_better(entry, existing->second)) {
            existing->second = std::move(entry);
        }
        persisted_descriptors.push_back(dirty_descriptor);
    }
    for (std::uint64_t descriptor : stale_descriptors) {
        state.archive_dirty.erase(descriptor);
    }
    if (exact_entries.empty()) {
        return state.archive_dirty.empty();
    }

    ArchiveDelta delta;
    delta.run_id = identity.run_id;
    delta.node_id = identity.node_id;
    // Reserve the sequence before publication.  A POSIX directory fsync can
    // fail after rename has already made the file visible; reusing that
    // sequence with different bytes would manufacture a merge conflict.
    delta.sequence = ++state.archive_sequence;
    delta.created_unix_ns = created_at;
    delta.topology = state.topology;
    delta.objective_version = GLOBAL_OBJECTIVE_VERSION;
    delta.topology_fingerprint = topology_fingerprint();
    delta.base_seed = static_cast<std::uint64_t>(
        static_cast<std::uint32_t>(options.seed));
    delta.entries.reserve(exact_entries.size());
    for (auto& item : exact_entries) {
        delta.entries.push_back(std::move(item.second));
    }
    const CommitResult result = commit_delta(options.global_dir, std::move(delta));
    if (!result) {
        std::cerr << "Durable archive delta failed for topology " << state.topology
                  << ": " << result.error << std::endl;
        return false;
    }
    for (std::uint64_t descriptor : persisted_descriptors) {
        state.archive_dirty.erase(descriptor);
    }
    if (!state.archive_dirty.empty()) {
        std::cerr << "Some archive cells could not be represented durably for topology "
                  << state.topology << std::endl;
        return false;
    }
    return true;
}

void save_global_topology_state(
    const LocalRepairOptions& options,
    const GlobalTopologyState& state,
    bool,
    int
) {
    if (!state.has_state || state.best_x.size() == 0) {
        return;
    }
    const std::filesystem::path dir = global_topology_dir(options, state.topology);
    std::filesystem::create_directories(dir);
    save_plane_state(dir / "resume.planes", state.best_x);
    export_plane_candidate((dir / "resume.obj").string().c_str(), state.best_x);
}

void write_global_leaderboard(
    const LocalRepairOptions& options,
    const std::vector<GlobalTopologyState>& states
) {
    std::vector<int> order(states.size());
    std::iota(order.begin(), order.end(), 0);
    std::sort(order.begin(), order.end(), [&](int a, int b) {
        if (states[a].has_state != states[b].has_state) {
            return states[a].has_state;
        }
        if (!states[a].has_state) {
            return states[a].topology < states[b].topology;
        }
        return better_global_metrics(states[a].best, states[b].best);
    });

    const std::filesystem::path root(options.global_dir);
    std::filesystem::create_directories(root);
    std::ofstream out(root / "leaderboard.tsv");
    out << "rank\ttopology\tC\tI\tdefects\tprecision\tvisits\ttrials"
        << "\tarchive_cells\tenergy\tdegeneracy\n";
    int rank = 1;
    for (int ix : order) {
        const GlobalTopologyState& state = states[ix];
        out << rank++ << "\t" << state.topology << "\t";
        if (state.has_state) {
            out << state.best.crossings << "\t" << state.best.intersections << "\t"
                << global_defects(state.best) << "\t"
                << (state.best.precise ? "dd31" : "double") << "\t"
                << state.visits << "\t"
                << state.trials << "\t" << state.archive.size() << "\t"
                << std::setprecision(17) << state.best.energy << "\t"
                << state.best.degeneracy_penalty;
        } else {
            out << "-\t-\t-\t-\t0\t0\t0\t-\t-";
        }
        out << "\n";
    }
}

bool write_run_manifest(
    const std::filesystem::path& run_directory,
    const LocalRepairOptions& options,
    const szilassi::checkpoint::RunIdentity& identity,
    int effective_cuda_chains,
    int effective_cuda_session_cache
) {
    const std::filesystem::path manifest_path = run_directory / "run.tsv";
    std::ofstream out(manifest_path, std::ios::binary | std::ios::trunc);
    if (!out) {
        std::cerr << "Cannot write run manifest: " << manifest_path << std::endl;
        return false;
    }
    out << "format\tszilassi-global-search-v4\n"
        << "run_id\t" << identity.run_id << "\n"
        << "node_id\t" << identity.node_id << "\n"
        << "seed\t" << options.seed << "\n"
        << "topology_from\t" << options.topology_from << "\n"
        << "topology_to\t" << options.topology_to << "\n"
        << "backend\t" << (options.use_cuda ? "cuda-fp32" : "cpu-double") << "\n"
        << "objective_version\t" << GLOBAL_OBJECTIVE_VERSION << "\n"
        << "cuda_math\tstandard-fp32\n"
        << "cuda_chains_requested\t" << options.cuda_chains << "\n"
        << "cuda_chains_effective\t" << effective_cuda_chains << "\n"
        << "cuda_iterations_per_batch\t" << options.cuda_iterations << "\n"
        << "cuda_depth_batches\t6\n"
        << "cuda_session_cache\t" << effective_cuda_session_cache << "\n"
        << "algorithm\thybrid-quality-diversity-v1\n"
        << "control_baseline_fraction\t0.25\n"
        << "strategy_weights\tbaseline:4,replica:3,adaptive:3,pbt:3,injected:3\n"
        << "fresh_fractions\tdepth:1/4,breadth:7/8,injected-protected\n"
        << "replica_exchange\tgroup:8,temperature-ratio:16\n"
        << "pbt\trotating-pairs,depth-chance:0.08,breadth-chance:0.04\n"
        << "verification_quotas\toverall:128,per-strategy:max(8,overall/5),per-injected-seed:1\n"
        << "accounted_fp32_steps\tproposal-iterations-plus-injection-pbt;setup-retries-excluded\n"
        << "map_elites_bins\tdeterminant:8,edge:8,turn:8,extent:8\n"
        << "map_elites_max_cells_per_topology\t4096\n"
        << "cem\tdiagonal-weighted,covariance-adaptation:0.18,new-injected-only\n"
        << "spsa\tadam:6,cpu-double,fp32-roundtrip,final-canonical-dd-gate\n"
        << "topology_scheduler\tucb-plus-reward-plus-staleness,full-refresh-every-5\n"
        << "topology_scheduler_mode\t"
        << (options.prioritize_worst ? "worst-first" : "quality-first") << "\n"
        << "worst_priority_coefficient\t0.65\n"
        << "cpu_iterations_per_trial\t" << options.iterations << "\n"
        << "degeneracy_formula\tworst-plus-0.05-rest\n"
        << "degeneracy_weight\t" << std::setprecision(17)
        << options.degeneracy_weight << "\n";
    out.flush();
    if (!out) {
        std::cerr << "Cannot flush run manifest: " << manifest_path << std::endl;
        return false;
    }
    out.close();
    if (!out) {
        std::cerr << "Cannot close run manifest: " << manifest_path << std::endl;
        return false;
    }

#ifdef _WIN32
    const HANDLE handle = CreateFileW(
        manifest_path.c_str(),
        GENERIC_WRITE,
        FILE_SHARE_READ | FILE_SHARE_WRITE | FILE_SHARE_DELETE,
        nullptr,
        OPEN_EXISTING,
        FILE_ATTRIBUTE_NORMAL,
        nullptr);
    if (handle == INVALID_HANDLE_VALUE || !FlushFileBuffers(handle)) {
        std::cerr << "Cannot durably flush run manifest: " << manifest_path << std::endl;
        if (handle != INVALID_HANDLE_VALUE) {
            CloseHandle(handle);
        }
        return false;
    }
    CloseHandle(handle);
#endif
    return true;
}

cuda_search::Topology make_cuda_topology() {
    cuda_search::Topology topology;
    for (int vertex = 0; vertex < cuda_search::kVertexCount; ++vertex) {
        for (int component = 0; component < 3; ++component) {
            topology.vertex_planes[vertex][component] = static_cast<std::uint8_t>(
                g_tris[static_cast<size_t>(vertex)][static_cast<size_t>(component)]);
        }
    }
    for (int face = 0; face < cuda_search::kPolygonCount; ++face) {
        for (int vertex = 0; vertex < cuda_search::kPolygonVertexCount; ++vertex) {
            topology.polygons[face][vertex] = static_cast<std::uint8_t>(
                g_polys[static_cast<size_t>(face)][static_cast<size_t>(vertex)]);
        }
    }
    for (int edge = 0; edge < cuda_search::kEdgeCount; ++edge) {
        topology.edges[edge][0] = static_cast<std::uint8_t>(g_edges[edge].first);
        topology.edges[edge][1] = static_cast<std::uint8_t>(g_edges[edge].second);
    }
    return topology;
}

struct CudaSessionSlot {
    int topology = -1;
    std::uint64_t last_used = 0;
    cuda_search::BatchSession session;
};

bool run_global_topology_round_cuda(
    const LocalRepairOptions& options,
    GlobalTopologyState& state,
    bool depth,
    const std::atomic<bool>& stop_requested,
    int round_seed,
    int effective_chain_count,
    cuda_search::BatchSession& session
) {
    const auto round_started_at = std::chrono::steady_clock::now();
    state.last_round = {};
    const GlobalMetrics previous_best = state.best;
    const bool previous_has_state = state.has_state;
    const bool session_was_initialized = session.initialized();
    if (!session_was_initialized) {
        cuda_search::SearchConfig config;
        config.chain_count = effective_chain_count;
        config.iterations_per_kernel = std::clamp(options.cuda_iterations, 1, 16);
        config.iterations_per_batch = options.cuda_iterations;
        config.shortlist_size = std::min(128, config.chain_count);
        config.seed = static_cast<std::uint64_t>(static_cast<std::uint32_t>(round_seed));
        config.initial_step = static_cast<float>(options.step);
        config.minimum_step = static_cast<float>(
            std::max(1e-7, options.step * options.min_step_ratio));
        config.cooling = static_cast<float>(options.beta);
        config.initial_temperature = static_cast<float>(options.temperature);
        config.minimum_temperature = static_cast<float>(options.temperature * 0.04);
        config.jump_chance = static_cast<float>(options.jump_chance);
        config.initial_state_jitter = static_cast<float>(options.step * 0.10);
        config.stagnation_iterations = std::max(256, options.stagnation);
        config.degeneracy_weight = static_cast<float>(options.degeneracy_weight);

        std::vector<cuda_search::PlaneState> initial_states = make_hybrid_seed_pool(
            state,
            static_cast<std::uint64_t>(static_cast<std::uint32_t>(round_seed)) ^
                0x6a09e667f3bcc909ULL,
            std::min(128, effective_chain_count));
        if (initial_states.empty() && state.has_state &&
            state.best_x.size() == GLOBAL_PLANE_VALUE_COUNT) {
            cuda_search::PlaneState initial;
            for (int i = 0; i < GLOBAL_PLANE_VALUE_COUNT; ++i) {
                initial.values[static_cast<size_t>(i)] =
                    static_cast<float>(state.best_x[i]);
            }
            initial_states.push_back(initial);
        }
        std::string initialize_error;
        if (!session.initialize(
                make_cuda_topology(), config, initial_states, initialize_error)) {
            std::cerr << "CUDA session initialization failed for topology "
                      << state.topology << ": " << initialize_error << std::endl;
            state.last_round.backend_error = true;
            return false;
        }
    }

    const int requested_batches = depth ? 6 : 1;
    int completed_batches = 0;
    double kernel_milliseconds = 0.0;
    double transfer_milliseconds = 0.0;
    std::uint64_t evaluated_states = 0;
    std::uint64_t new_trials = session_was_initialized
        ? 0
        : static_cast<std::uint64_t>(effective_chain_count);
    std::string device_name;
    PlaneEvaluationScratch scratch;
    for (int batch = 0; batch < requested_batches; ++batch) {
        if (stop_requested.load(std::memory_order_relaxed)) {
            break;
        }
        cuda_search::BatchRunConfig run_config;
        run_config.step_scale = depth ? 1.0f : 1.5f;
        // Exactly one quarter remains the unchanged simulated-annealing
        // control. The other disjoint cohorts test complementary mechanisms.
        run_config.baseline_weight = 4;
        run_config.replica_exchange_weight = 3;
        run_config.adaptive_move_weight = 3;
        run_config.pbt_weight = 3;
        run_config.injected_weight = 3;
        run_config.replica_group_size = 8;
        run_config.replica_temperature_ratio = 16.0f;
        run_config.pbt_exploit_chance = depth ? 0.08f : 0.04f;
        run_config.pbt_state_jitter = static_cast<float>(
            std::max(0.002, options.step * (depth ? 0.025 : 0.05)));
        run_config.injected_state_jitter = static_cast<float>(
            std::max(0.002, options.step * (depth ? 0.018 : 0.04)));
        run_config.injected_states = make_hybrid_seed_pool(
            state,
            static_cast<std::uint64_t>(static_cast<std::uint32_t>(round_seed)) +
                static_cast<std::uint64_t>(batch + 1) * 0x9e3779b97f4a7c15ULL,
            64);
        if (batch == 0 && (previous_has_state || session_was_initialized)) {
            run_config.fresh_numerator = depth ? 1 : 7;
            run_config.fresh_denominator = depth ? 4 : 8;
        }
        const cuda_search::BatchResult gpu = session.run(run_config);
        if (!gpu.success) {
            std::cerr << "CUDA batch failed for topology " << state.topology
                      << ": " << gpu.error << std::endl;
            state.last_round.backend_error = true;
            session.reset();
            break;
        }
        device_name = gpu.device_name;
        g_search_device_id = gpu.device_name;
        kernel_milliseconds += gpu.kernel_milliseconds;
        transfer_milliseconds += gpu.transfer_milliseconds;
        evaluated_states += gpu.evaluated_states;
        state.last_round.replica_swaps_attempted +=
            gpu.replica_exchange_attempts;
        state.last_round.replica_swaps_accepted +=
            gpu.replica_exchange_accepts;
        constexpr std::array<int, 5> kStrategyWeights = {4, 3, 3, 3, 3};
        constexpr int kTotalStrategyWeight = 16;
        for (int chain = 0; chain < effective_chain_count; ++chain) {
            const int residue = chain % kTotalStrategyWeight;
            int boundary = 0;
            for (int strategy = 0; strategy < 5; ++strategy) {
                boundary += kStrategyWeights[static_cast<std::size_t>(strategy)];
                if (residue < boundary) {
                    state.last_round.strategy_evaluated[static_cast<std::size_t>(strategy)] +=
                        static_cast<std::uint64_t>(options.cuda_iterations);
                    break;
                }
            }
        }
        state.last_round.strategy_evaluated[3] += gpu.pbt_exploits;
        state.last_round.strategy_evaluated[4] += gpu.injected_chains;
        new_trials += gpu.fresh_chains + gpu.injected_chains;
        completed_batches += 1;

        std::vector<GlobalVerifiedCandidate> batch_verified;
        for (const cuda_search::Candidate& candidate : gpu.shortlist) {
            VectorXd x(GLOBAL_PLANE_VALUE_COUNT);
            for (int i = 0; i < GLOBAL_PLANE_VALUE_COUNT; ++i) {
                x[i] = static_cast<double>(
                    candidate.state.values[static_cast<size_t>(i)]);
            }
            if (!valid_plane_state(x)) {
                continue;
            }
            GlobalMetrics verified = evaluate_global_state(x, true, scratch);
            if (!std::isfinite(verified.energy)) {
                continue;
            }
            const std::size_t strategy = std::min<std::size_t>(
                static_cast<std::size_t>(candidate.strategy),
                state.last_round.strategy_verified.size() - 1);
            state.last_round.verified_candidates += 1;
            state.last_round.strategy_verified[strategy] += 1;
            const bool archive_improved = add_global_archive_elite(
                state,
                verified,
                x,
                {},
                static_cast<std::uint64_t>(static_cast<std::uint32_t>(round_seed)),
                candidate.iterations,
                true);
            if (archive_improved) {
                state.last_round.archive_improvements += 1;
                state.last_round.strategy_archive_improvements[strategy] += 1;
            }
            batch_verified.push_back(GlobalVerifiedCandidate{
                verified,
                x,
                static_cast<std::uint32_t>(strategy),
                candidate.chain_id});
            if (!state.has_state || better_global_metrics(verified, state.best)) {
                state.best = verified;
                state.best_x = x;
                state.has_state = true;
                state.last_round.strategy_global_improvements[strategy] += 1;
            }
        }
        update_diagonal_cem(state, batch_verified);
    }
    if (completed_batches == 0) {
        return false;
    }

    if (depth && state.has_state && !state.last_round.backend_error &&
        !stop_requested.load(std::memory_order_relaxed)) {
        state.last_round.spsa_attempted = true;
        VectorXd refined;
        GlobalMetrics refined_metrics;
        if (refine_with_spsa_adam(
                state.best_x,
                static_cast<std::uint64_t>(static_cast<std::uint32_t>(round_seed)) ^
                    0xbb67ae8584caa73bULL,
                6,
                refined,
                refined_metrics) &&
            round_trip_global_state_to_fp32(refined, refined_metrics)) {
            const bool archive_improved = add_global_archive_elite(
                state,
                refined_metrics,
                refined,
                {},
                static_cast<std::uint64_t>(static_cast<std::uint32_t>(round_seed)),
                0,
                true);
            if (archive_improved) {
                state.last_round.archive_improvements += 1;
            }
            if (better_global_metrics(refined_metrics, state.best)) {
                state.best = refined_metrics;
                state.best_x = refined;
                state.last_round.spsa_accepted = true;
            } else if (archive_improved) {
                state.last_round.spsa_accepted = true;
            }
        }
    }

    state.visits += 1;
    state.run_visits += 1;
    state.trials += new_trials;
    state.run_trials += new_trials;
    state.iterations += evaluated_states;
    state.run_iterations += evaluated_states;
    state.last_round.completed = true;
    state.last_round.new_trials = new_trials;
    state.last_round.evaluated_states = evaluated_states;
    state.last_round.kernel_milliseconds = kernel_milliseconds;
    state.last_round.transfer_milliseconds = transfer_milliseconds;
    state.last_round.wall_milliseconds = std::chrono::duration<double, std::milli>(
        std::chrono::steady_clock::now() - round_started_at).count();
    std::cout << "  CUDA " << device_name
              << ": " << std::fixed << std::setprecision(1)
              << kernel_milliseconds << " ms kernel, "
              << transfer_milliseconds << " ms transfer, "
              << evaluated_states << " accounted FP32 search steps in "
              << completed_batches << " batch(es)"
              << std::defaultfloat << std::setprecision(6) << std::endl;
    state.last_round.improved = state.has_state &&
        (!previous_has_state || better_global_metrics(state.best, previous_best));
    return state.last_round.improved;
}

bool run_global_topology_round(
    const LocalRepairOptions& options,
    GlobalTopologyState& state,
    bool depth,
    int worker_count,
    const std::atomic<bool>& stop_requested,
    int round_seed,
    int effective_cuda_chains,
    cuda_search::BatchSession* cuda_session
) {
    if (options.use_cuda) {
        if (cuda_session == nullptr) {
            return false;
        }
        return run_global_topology_round_cuda(
            options,
            state,
            depth,
            stop_requested,
            round_seed,
            effective_cuda_chains,
            *cuda_session);
    }
    const auto round_started_at = std::chrono::steady_clock::now();
    state.last_round = {};
    const GlobalMetrics previous_best = state.best;
    const bool previous_has_state = state.has_state;
    GlobalMetrics shared_best = state.best;
    VectorXd shared_best_x = state.best_x;
    bool shared_has_state = state.has_state;

    const int trial_count = depth
        ? std::max(8, worker_count * 2)
        : std::max(8, worker_count);
    const int iterations = depth
        ? std::max(1000, options.iterations)
        : std::max(1000, options.iterations / 3);
    std::atomic<int> next_trial{0};
    std::atomic<int> completed{0};
    std::atomic<std::uint64_t> completed_iterations{0};
    std::mutex best_mutex;
    std::vector<std::thread> workers;
    workers.reserve(worker_count);

    auto worker = [&]() {
        while (!stop_requested.load(std::memory_order_relaxed)) {
            const int trial = next_trial.fetch_add(1, std::memory_order_relaxed);
            if (trial >= trial_count) {
                break;
            }
            VectorXd start_x;
            bool has_start = false;
            {
                std::lock_guard<std::mutex> lock(best_mutex);
                has_start = shared_has_state;
                if (has_start) {
                    start_x = shared_best_x;
                }
            }
            const bool start_fresh = !has_start || (depth ? trial % 4 == 0 : trial % 8 != 0);
            const std::uint64_t seed_value =
                static_cast<std::uint64_t>(static_cast<std::uint32_t>(round_seed)) +
                static_cast<std::uint64_t>(trial + 1) * 1000003ULL +
                static_cast<std::uint64_t>(state.topology + 1) * 104729ULL;
            const int seed = static_cast<int>(1 + seed_value % 2147483646ULL);
            GlobalTrialResult result = run_global_trial(
                options, start_x, start_fresh, seed, iterations, stop_requested);
            if (std::isfinite(result.best.energy) &&
                !round_trip_global_state_to_fp32(result.best_x, result.best)) {
                result.best = GlobalMetrics{};
            }
            {
                std::lock_guard<std::mutex> lock(best_mutex);
                if (better_global_metrics(result.best, shared_best)) {
                    shared_best = result.best;
                    shared_best_x = result.best_x;
                    shared_has_state = true;
                }
            }
            completed.fetch_add(1, std::memory_order_relaxed);
            completed_iterations.fetch_add(
                static_cast<std::uint64_t>(result.iterations),
                std::memory_order_relaxed);
        }
    };

    const int actual_workers = std::min(worker_count, trial_count);
    for (int i = 0; i < actual_workers; ++i) {
        workers.emplace_back(worker);
    }
    for (std::thread& thread : workers) {
        thread.join();
    }

    state.visits += 1;
    state.run_visits += 1;
    const std::uint64_t round_trials =
        static_cast<std::uint64_t>(completed.load(std::memory_order_relaxed));
    const std::uint64_t round_iterations =
        completed_iterations.load(std::memory_order_relaxed);
    state.trials += round_trials;
    state.iterations += round_iterations;
    state.run_trials += round_trials;
    state.run_iterations += round_iterations;
    state.has_state = shared_has_state;
    state.best = shared_best;
    state.best_x = shared_best_x;
    if (state.has_state && add_global_archive_elite(
            state,
            state.best,
            state.best_x,
            {},
            static_cast<std::uint64_t>(static_cast<std::uint32_t>(round_seed)),
            state.iterations,
            true)) {
        state.last_round.archive_improvements = 1;
    }
    state.last_round.completed = true;
    state.last_round.new_trials = round_trials;
    state.last_round.evaluated_states = round_iterations;
    state.last_round.verified_candidates = round_trials;
    state.last_round.strategy_evaluated[0] = round_iterations;
    state.last_round.strategy_verified[0] = round_trials;
    state.last_round.wall_milliseconds = std::chrono::duration<double, std::milli>(
        std::chrono::steady_clock::now() - round_started_at).count();
    state.last_round.improved = state.has_state &&
        (!previous_has_state || better_global_metrics(state.best, previous_best));
    return state.last_round.improved;
}

int global_search_all(const LocalRepairOptions& options) {
    if (!wide_real_self_test()) {
        std::cerr << "WideReal self-test failed; refusing high-precision search." << std::endl;
        return 2;
    }
    const std::filesystem::path root(options.global_dir);
    std::filesystem::create_directories(root);
    const szilassi::checkpoint::RunIdentity run_identity =
        szilassi::checkpoint::make_run_identity();
    g_global_degeneracy_weight = options.degeneracy_weight;
#ifdef _WIN32
    SetPriorityClass(GetCurrentProcess(), BELOW_NORMAL_PRIORITY_CLASS);
    SetThreadPriority(GetCurrentThread(), THREAD_PRIORITY_BELOW_NORMAL);
#endif

    cuda_search::BackendInfo cuda_backend;
    int effective_cuda_chains = 0;
    if (options.use_cuda) {
        cuda_backend = cuda_search::query_backend(0);
        if (!cuda_backend.available) {
            std::cerr << "CUDA backend is unavailable: " << cuda_backend.error << std::endl;
            std::cerr << "Install CUDA Toolkit 13.3 with Visual Studio integration, then rebuild x64."
                      << std::endl;
            return 3;
        }
        effective_cuda_chains = options.cuda_chains > 0
            ? options.cuda_chains
            : cuda_backend.recommended_chain_count;
        if (effective_cuda_chains <= 0) {
            std::cerr << "CUDA backend did not provide a valid chain count." << std::endl;
            return 3;
        }
        g_search_device_id = cuda_backend.device_name;
        std::cout << "CUDA device: " << cuda_backend.device_name << std::endl;
        std::cout << "CUDA SMs: " << cuda_backend.multiprocessor_count
                  << ", resident blocks/SM: "
                  << cuda_backend.active_blocks_per_multiprocessor << std::endl;
        std::cout << "GPU scheduling: low-priority stream, no artificial throttle"
                  << std::endl;
    }

    const unsigned int hardware_threads = std::thread::hardware_concurrency();
    const int worker_count = std::max(
        1,
        options.threads > 0
            ? options.threads
            : static_cast<int>(hardware_threads == 0 ? 1 : hardware_threads));
    const auto started_at = std::chrono::steady_clock::now();
    const auto deadline = options.time_limit_seconds > 0
        ? started_at + std::chrono::seconds(options.time_limit_seconds)
        : std::chrono::steady_clock::time_point::max();
    std::atomic<bool> stop_requested{false};
    std::atomic<bool> watcher_done{false};
    std::thread watcher([&]() {
        while (!watcher_done.load(std::memory_order_relaxed)) {
            if (std::chrono::steady_clock::now() >= deadline) {
                stop_requested.store(true, std::memory_order_relaxed);
            }
            if (!options.stop_file_path.empty() &&
                std::filesystem::exists(options.stop_file_path)) {
                stop_requested.store(true, std::memory_order_relaxed);
            }
            std::this_thread::sleep_for(std::chrono::milliseconds(250));
        }
    });

    std::vector<GlobalTopologyState> states(NUM_TOPOLOGIES);
    for (int topology = 0; topology < NUM_TOPOLOGIES; ++topology) {
        states[topology].topology = topology;
        if (load_global_topology_context(topology)) {
            load_global_topology_state(options, states[topology]);
        }
    }
    load_mergeable_checkpoints(options, states);
    write_global_leaderboard(options, states);

    std::vector<int> active_topologies;
    active_topologies.reserve(options.topology_to - options.topology_from + 1);
    for (int topology = options.topology_from; topology <= options.topology_to; ++topology) {
        active_topologies.push_back(topology);
        states[topology].scheduler_pulls = std::min<std::uint64_t>(
            8,
            static_cast<std::uint64_t>(std::floor(
                std::log2(static_cast<double>(states[topology].visits) + 1.0))));
    }

    const std::filesystem::path run_directory =
        root / "runs" / run_identity.run_id;
    std::filesystem::create_directories(run_directory);
    const int effective_cuda_session_cache = options.use_cuda
        ? std::min(
            CUDA_SESSION_CACHE_LIMIT,
            static_cast<int>(active_topologies.size()))
        : 0;
    if (!write_run_manifest(
            run_directory,
            options,
            run_identity,
            effective_cuda_chains,
            effective_cuda_session_cache)) {
        stop_requested.store(true, std::memory_order_relaxed);
        watcher_done.store(true, std::memory_order_relaxed);
        watcher.join();
        return 2;
    }
    std::ofstream run_log(run_directory / "run.log", std::ios::app);
    run_log << "\n=== global-search seed " << options.seed
            << ", run " << run_identity.run_id
            << ", node " << run_identity.node_id
            << ", threads " << worker_count
            << ", scheduler "
            << (options.prioritize_worst ? "worst-first" : "quality-first")
            << ", minutes " << (options.time_limit_seconds / 60.0)
            << " ===\n";
    std::ofstream metrics_log(
        run_directory / "metrics.tsv",
        std::ios::binary | std::ios::trunc);
    if (!metrics_log) {
        std::cerr << "Cannot create strategy telemetry: "
                  << (run_directory / "metrics.tsv") << std::endl;
        stop_requested.store(true, std::memory_order_relaxed);
        watcher_done.store(true, std::memory_order_relaxed);
        watcher.join();
        return 2;
    }
    metrics_log
        << "round\tunix_ns\tphase\ttopology\tdepth\tbandit_score\tworstness"
        << "\tbefore_C\tbefore_I\tafter_C\tafter_I\timproved\tbackend_error"
        << "\tarchive_cells\tarchive_improvements\tverified\taccounted_fp32_steps"
        << "\tnew_trials\tkernel_ms\ttransfer_ms\twall_ms"
        << "\trex_attempts\trex_accepts\tspsa_attempted\tspsa_accepted"
        << "\tsteps_baseline\tsteps_replica\tsteps_adaptive\tsteps_pbt\tsteps_injected"
        << "\tverified_baseline\tverified_replica\tverified_adaptive"
        << "\tverified_pbt\tverified_injected"
        << "\tarchive_baseline\tarchive_replica\tarchive_adaptive"
        << "\tarchive_pbt\tarchive_injected"
        << "\tglobal_baseline\tglobal_replica\tglobal_adaptive"
        << "\tglobal_pbt\tglobal_injected\n";
    metrics_log.flush();

    std::cout << "===================" << std::endl;
    std::cout << "Mode      : global-search" << std::endl;
    std::cout << "Seed      : " << options.seed
              << " (saved in " << (run_directory / "run.tsv").string() << ")"
              << std::endl;
    std::cout << "Topologies: " << options.topology_from << ".."
              << options.topology_to << " (" << active_topologies.size() << ")" << std::endl;
    std::cout << "Scheduler : "
              << (options.prioritize_worst
                    ? "priority to the worst current topologies"
                    : "priority to the most promising topologies")
              << std::endl;
    std::cout << "Degenerate: worst barrier + 5% of the rest, weight "
              << options.degeneracy_weight << std::endl;
    std::cout << "Start     : saved MAP-Elites/CEM seeds plus independent random starts"
              << std::endl;
    std::cout << "Coordinates: "
              << (options.use_cuda ? "CUDA FP32 search" : "CPU double search")
              << std::endl;
    std::cout << "Near goal : WideReal double-double (~"
              << WideReal::decimal_digits << " digits)" << std::endl;
    std::cout << "Resume dir: " << options.global_dir << std::endl;
    if (options.use_cuda) {
        std::cout << "CUDA chains: " << effective_cuda_chains
                  << (options.cuda_chains == 0 ? " (automatic)" : " (manual)")
                  << std::endl;
        std::cout << "CUDA iters/batch: " << options.cuda_iterations << std::endl;
    } else {
        std::cout << "Threads   : " << worker_count << std::endl;
        std::cout << "Depth iters/trial: " << options.iterations << std::endl;
    }
    if (options.time_limit_seconds > 0) {
        std::cout << "Time limit: " << options.time_limit_seconds << " sec" << std::endl;
    }
    std::cout << "===================" << std::endl;

    int completed_rounds = 0;
    bool found = false;
    bool fatal_search_error = false;
    std::vector<CudaSessionSlot> cuda_sessions(effective_cuda_session_cache);
    std::uint64_t cuda_session_use_counter = 0;
    auto acquire_cuda_session = [&](int topology) -> cuda_search::BatchSession* {
        if (!options.use_cuda || cuda_sessions.empty()) {
            return nullptr;
        }
        ++cuda_session_use_counter;
        for (CudaSessionSlot& slot : cuda_sessions) {
            if (slot.topology == topology) {
                slot.last_used = cuda_session_use_counter;
                return &slot.session;
            }
        }
        CudaSessionSlot* selected = nullptr;
        for (CudaSessionSlot& slot : cuda_sessions) {
            if (slot.topology < 0) {
                selected = &slot;
                break;
            }
            if (selected == nullptr || slot.last_used < selected->last_used) {
                selected = &slot;
            }
        }
        selected->session.reset();
        selected->topology = topology;
        selected->last_used = cuda_session_use_counter;
        return &selected->session;
    };
    auto persist_topology_if_due = [&](GlobalTopologyState& state,
                                       std::chrono::steady_clock::time_point now,
                                       bool force,
                                       szilassi::checkpoint::CheckpointReason forced_reason) {
        if (state.run_visits == 0) {
            return true;
        }
        const bool checkpoint_changed =
            state.run_visits > state.last_checkpointed_run_visits;
        const bool archive_changed = !state.archive_dirty.empty();
        if (!checkpoint_changed && !archive_changed) {
            return true;
        }
        const bool due = force ||
            state.last_checkpoint_at.time_since_epoch().count() == 0 ||
            now - state.last_checkpoint_at >=
                std::chrono::seconds(options.checkpoint_seconds);
        if (!due) {
            return true;
        }
        if (!load_global_topology_context(state.topology)) {
            std::cerr << "Cannot restore topology context for durable save: "
                      << state.topology << std::endl;
            return false;
        }

        bool checkpoint_saved = true;
        if (checkpoint_changed) {
            const auto reason = force
                ? forced_reason
                : (state.pending_checkpoint_improvement
                    ? szilassi::checkpoint::CheckpointReason::Improvement
                    : szilassi::checkpoint::CheckpointReason::Periodic);
            checkpoint_saved = commit_mergeable_checkpoint(
                options,
                run_identity,
                state,
                reason);
            if (checkpoint_saved) {
                state.last_checkpointed_run_visits = state.run_visits;
                state.pending_checkpoint_improvement = false;
                save_global_topology_state(options, state, true, options.seed);
            }
        }
        const bool archive_saved = commit_mergeable_archive_delta(
            options,
            run_identity,
            state);
        if (checkpoint_saved && archive_saved && state.archive_dirty.empty()) {
            state.last_checkpoint_at = now;
            return true;
        }
        return false;
    };
    auto sweep_due_topologies = [&](int restore_topology,
                                    bool force,
                                    szilassi::checkpoint::CheckpointReason forced_reason) {
        const auto now = std::chrono::steady_clock::now();
        bool success = true;
        for (int item : active_topologies) {
            if (!persist_topology_if_due(states[item], now, force, forced_reason)) {
                success = false;
            }
        }
        if (restore_topology >= 0) {
            if (!load_global_topology_context(restore_topology)) {
                std::cerr << "Cannot restore active topology context after durable sweep: "
                          << restore_topology << std::endl;
                success = false;
            }
        }
        return success;
    };
    auto run_round = [&](int topology, bool depth, const char* phase, int position, int total) {
        if (stop_requested.load(std::memory_order_relaxed)) {
            return;
        }
        if (!load_global_topology_context(topology)) {
            std::cout << "Topology " << topology << " failed combinatorial setup." << std::endl;
            return;
        }
        GlobalTopologyState& state = states[topology];
        const bool had_before = state.has_state;
        const GlobalMetrics before = state.best;
        const std::uint64_t total_scheduler_pulls = std::accumulate(
            active_topologies.begin(),
            active_topologies.end(),
            std::uint64_t{0},
            [&](std::uint64_t sum, int item) {
                return sum + states[item].scheduler_pulls;
            });
        double current_best_severity = std::numeric_limits<double>::infinity();
        double current_worst_severity = -std::numeric_limits<double>::infinity();
        int current_best_defects = std::numeric_limits<int>::max() / 4;
        for (int item : active_topologies) {
            if (states[item].has_state) {
                const double severity = topology_scheduler_severity(states[item].best);
                current_best_severity = std::min(current_best_severity, severity);
                current_worst_severity = std::max(current_worst_severity, severity);
                current_best_defects = std::min(
                    current_best_defects,
                    global_defects(states[item].best));
            }
        }
        if (!std::isfinite(current_best_severity)) {
            current_best_severity = 0.0;
            current_worst_severity = 0.0;
            current_best_defects = 0;
        }
        const double selection_score = topology_bandit_score(
            state,
            total_scheduler_pulls,
            static_cast<std::uint64_t>(completed_rounds + 1),
            current_best_severity,
            current_worst_severity,
            current_best_defects,
            options.prioritize_worst);
        const double selection_worstness = topology_scheduler_worstness(
            state,
            current_best_severity,
            current_worst_severity);
        const int round_seed = static_cast<int>(
            1 + (static_cast<std::uint64_t>(static_cast<std::uint32_t>(options.seed)) +
                 static_cast<std::uint64_t>(state.visits + 1) * 15485863ULL +
                 static_cast<std::uint64_t>(topology + 1) * 32452843ULL) %
                2147483646ULL);
        std::cout << phase << " " << position << "/" << total
                  << " | topology " << topology
                  << " | " << (depth ? "depth" : "breadth") << std::endl;
        cuda_search::BatchSession* cuda_session = acquire_cuda_session(topology);
        const bool improved = run_global_topology_round(
            options,
            state,
            depth,
            worker_count,
            stop_requested,
            round_seed,
            effective_cuda_chains,
            cuda_session);
        if (!state.last_round.completed) {
            if (state.last_round.backend_error) {
                fatal_search_error = true;
                stop_requested.store(true, std::memory_order_relaxed);
                std::cerr << "CUDA search failed; stopping instead of retrying an invalid round."
                          << std::endl;
            }
            return;
        }
        completed_rounds += 1;
        update_topology_bandit_reward(
            state,
            had_before,
            before,
            state.last_round.archive_improvements,
            static_cast<std::uint64_t>(completed_rounds));
        state.pending_checkpoint_improvement =
            state.pending_checkpoint_improvement || improved;
        if (state.last_round.backend_error) {
            fatal_search_error = true;
            stop_requested.store(true, std::memory_order_relaxed);
        }
        if (!sweep_due_topologies(
                topology,
                false,
                szilassi::checkpoint::CheckpointReason::Periodic)) {
            std::cerr << "Periodic durable save failed; stopping search safely."
                      << std::endl;
            stop_requested.store(true, std::memory_order_relaxed);
        }
        write_global_leaderboard(options, states);

        const double elapsed_minutes = std::chrono::duration<double>(
            std::chrono::steady_clock::now() - started_at).count() / 60.0;
        std::cout << "  result topology " << topology << ": C/I "
                  << state.best.crossings << "/" << state.best.intersections
                  << ", defects " << global_defects(state.best)
                  << ", visits " << state.visits
                  << (state.best.precise ? ", DD31" : ", double")
                  << (improved ? " NEW BEST" : "")
                  << ", archive " << state.archive.size()
                  << ", bandit " << std::fixed << std::setprecision(3)
                  << selection_score
                  << ", elapsed " << std::fixed << std::setprecision(1)
                  << elapsed_minutes << " min" << std::defaultfloat << std::setprecision(6)
                  << std::endl;
        run_log << phase << "\t" << topology << "\t"
                << state.best.crossings << "\t" << state.best.intersections << "\t"
                << state.visits << "\t" << state.trials << "\t"
                << std::setprecision(17) << state.best.energy << "\n";
        run_log.flush();
        metrics_log << completed_rounds << "\t"
                    << szilassi::checkpoint::unix_time_ns_now() << "\t"
                    << phase << "\t" << topology << "\t" << (depth ? 1 : 0) << "\t"
                    << std::setprecision(17) << selection_score << "\t"
                    << selection_worstness << "\t"
                    << (had_before ? before.crossings : -1) << "\t"
                    << (had_before ? before.intersections : -1) << "\t"
                    << (state.has_state ? state.best.crossings : -1) << "\t"
                    << (state.has_state ? state.best.intersections : -1) << "\t"
                    << (improved ? 1 : 0) << "\t"
                    << (state.last_round.backend_error ? 1 : 0) << "\t"
                    << state.archive.size() << "\t"
                    << state.last_round.archive_improvements << "\t"
                    << state.last_round.verified_candidates << "\t"
                    << state.last_round.evaluated_states << "\t"
                    << state.last_round.new_trials << "\t"
                    << state.last_round.kernel_milliseconds << "\t"
                    << state.last_round.transfer_milliseconds << "\t"
                    << state.last_round.wall_milliseconds << "\t"
                    << state.last_round.replica_swaps_attempted << "\t"
                    << state.last_round.replica_swaps_accepted << "\t"
                    << (state.last_round.spsa_attempted ? 1 : 0) << "\t"
                    << (state.last_round.spsa_accepted ? 1 : 0);
        for (int strategy = 0; strategy < 5; ++strategy) {
            metrics_log << "\t" << state.last_round.strategy_evaluated[strategy];
        }
        for (int strategy = 0; strategy < 5; ++strategy) {
            metrics_log << "\t" << state.last_round.strategy_verified[strategy];
        }
        for (int strategy = 0; strategy < 5; ++strategy) {
            metrics_log << "\t"
                        << state.last_round.strategy_archive_improvements[strategy];
        }
        for (int strategy = 0; strategy < 5; ++strategy) {
            metrics_log << "\t"
                        << state.last_round.strategy_global_improvements[strategy];
        }
        metrics_log << "\n";
        metrics_log.flush();

        if (state.best.crossings == 0 && state.best.intersections == 0) {
            int saved_crossings = 0;
            int saved_intersections = 0;
            const std::filesystem::path found_path =
                root / ("FOUND_topology_" + std::to_string(topology) + ".obj");
            if (export_and_validate_found_candidate(
                    found_path, state.best_x, saved_crossings, saved_intersections)) {
                std::cout << "FOUND robust 0/0: " << found_path.string() << std::endl;
                run_log << "FOUND\t" << topology << "\t" << found_path.string() << "\n";
                run_log.flush();
                found = true;
                stop_requested.store(true, std::memory_order_relaxed);
            }
        }
    };

    std::vector<int> breadth_order = active_topologies;
    const int active_count = static_cast<int>(breadth_order.size());
    const int rotation = ((options.seed % active_count) + active_count) % active_count;
    std::rotate(breadth_order.begin(), breadth_order.begin() + rotation, breadth_order.end());
    std::vector<int> missing;
    for (int topology : breadth_order) {
        if (!states[topology].has_state) {
            missing.push_back(topology);
        }
    }
    for (size_t i = 0; i < missing.size() && !stop_requested.load(std::memory_order_relaxed); ++i) {
        run_round(missing[i], false, "BREADTH", static_cast<int>(i + 1), static_cast<int>(missing.size()));
    }

    int cycle = 0;
    while (!stop_requested.load(std::memory_order_relaxed)) {
        std::vector<int> ranked = active_topologies;
        const std::uint64_t total_scheduler_pulls = std::accumulate(
            active_topologies.begin(),
            active_topologies.end(),
            std::uint64_t{0},
            [&](std::uint64_t sum, int topology) {
                return sum + states[topology].scheduler_pulls;
            });
        double best_severity = std::numeric_limits<double>::infinity();
        double worst_severity = -std::numeric_limits<double>::infinity();
        int best_defects = std::numeric_limits<int>::max() / 4;
        for (int topology : active_topologies) {
            if (states[topology].has_state) {
                const double severity = topology_scheduler_severity(states[topology].best);
                best_severity = std::min(best_severity, severity);
                worst_severity = std::max(worst_severity, severity);
                best_defects = std::min(
                    best_defects,
                    global_defects(states[topology].best));
            }
        }
        if (!std::isfinite(best_severity)) {
            best_severity = 0.0;
            worst_severity = 0.0;
            best_defects = 0;
        }
        std::unordered_map<int, double> bandit_scores;
        for (int topology : active_topologies) {
            bandit_scores[topology] = topology_bandit_score(
                states[topology],
                total_scheduler_pulls,
                static_cast<std::uint64_t>(completed_rounds + 1),
                best_severity,
                worst_severity,
                best_defects,
                options.prioritize_worst);
        }
        std::sort(ranked.begin(), ranked.end(), [&](int a, int b) {
            if (bandit_scores[a] != bandit_scores[b]) {
                return bandit_scores[a] > bandit_scores[b];
            }
            if (states[a].has_state != states[b].has_state) {
                return states[a].has_state;
            }
            if (!states[a].has_state) {
                return a < b;
            }
            if (better_global_metrics(states[a].best, states[b].best)) {
                return true;
            }
            if (better_global_metrics(states[b].best, states[a].best)) {
                return false;
            }
            return a < b;
        });

        const bool exploration_cycle = cycle % 5 == 4;
        int keep = active_count;
        if (!exploration_cycle) {
            keep = cycle == 0 ? 30 : (cycle == 1 ? 16 : 8);
        } else {
            RNG shuffle_rng(static_cast<RNG::result_type>(options.seed + cycle * 7919));
            std::shuffle(ranked.begin(), ranked.end(), shuffle_rng);
        }
        ranked.resize(std::min<int>(keep, static_cast<int>(ranked.size())));
        const char* phase = exploration_cycle ? "REFRESH" : "DEPTH";
        for (size_t i = 0; i < ranked.size() && !stop_requested.load(std::memory_order_relaxed); ++i) {
            run_round(
                ranked[i],
                !exploration_cycle,
                phase,
                static_cast<int>(i + 1),
                static_cast<int>(ranked.size()));
        }
        cycle += 1;
    }

    watcher_done.store(true, std::memory_order_relaxed);
    watcher.join();
    bool persistence_succeeded = sweep_due_topologies(
        -1,
        true,
        szilassi::checkpoint::CheckpointReason::Stop);
    if (!persistence_succeeded) {
        // A sequence is never reused, so one immediate retry is safe even if
        // an OS reported an error after publishing an immutable generation.
        persistence_succeeded = sweep_due_topologies(
            -1,
            true,
            szilassi::checkpoint::CheckpointReason::Stop);
    }
    if (persistence_succeeded && !options.stop_file_path.empty()) {
        std::error_code remove_error;
        std::filesystem::remove(options.stop_file_path, remove_error);
    } else if (!persistence_succeeded) {
        std::cerr << "Final durable save failed. The stop marker was retained; "
                  << "the process returns an error instead of claiming a clean stop."
                  << std::endl;
    }
    write_global_leaderboard(options, states);

    std::vector<int> ranked = active_topologies;
    std::sort(ranked.begin(), ranked.end(), [&](int a, int b) {
        if (states[a].has_state != states[b].has_state) {
            return states[a].has_state;
        }
        return states[a].has_state
            ? better_global_metrics(states[a].best, states[b].best)
            : a < b;
    });
    std::cout << "===================" << std::endl;
    std::cout << "Global search finished after " << completed_rounds << " topology rounds." << std::endl;
    if (!ranked.empty() && states[ranked.front()].has_state) {
        const GlobalTopologyState& leader = states[ranked.front()];
        std::cout << "Leader: topology " << leader.topology << ", C/I "
                  << leader.best.crossings << "/" << leader.best.intersections << std::endl;
    }
    std::cout << "Leaderboard: " << (root / "leaderboard.tsv").string() << std::endl;
    std::cout << "===================" << std::endl;
    if (!persistence_succeeded) {
        return 4;
    }
    return fatal_search_error ? 3 : 0;
}

int repair_local_shape(const LocalRepairOptions& options) {
    g_topology = options.topology;
    set_rand_seed(options.seed);

    const std::filesystem::path out_prefix_path(options.out_prefix);
    if (out_prefix_path.has_parent_path()) {
        std::filesystem::create_directories(out_prefix_path.parent_path());
    }
    const std::filesystem::path report_path(options.report_path);
    if (report_path.has_parent_path()) {
        std::filesystem::create_directories(report_path.parent_path());
    }

    Verts3D original;
    import_obj(options.obj_path.c_str(), original, g_polys);
    if (original.empty() || g_polys.empty()) {
        std::cerr << "Failed to load OBJ: " << options.obj_path << std::endl;
        return 2;
    }
    make_edges(g_polys, g_edges);
    Edges dual_edges;
    dual_graph(g_polys, g_tris, dual_edges);

    Verts3D current = original;
    Verts3D best = original;
    int current_c = 0;
    int current_i = 0;
    int best_c = 0;
    int best_i = 0;
    double current_loss = local_repair_objective(current, original, options.movable_vertices, current_c, current_i);
    double best_loss = local_repair_objective(best, original, options.movable_vertices, best_c, best_i);
    double step = options.step;
    double temperature = options.temperature;

    std::normal_distribution<double> normal(0.0, 1.0);
    std::uniform_real_distribution<double> uniform(0.0, 1.0);
    std::uniform_int_distribution<int> vertex_pick(0, (int)options.movable_vertices.size() - 1);

    std::ofstream report(options.report_path, std::ios::app);
    report << "# C++ local repair run\n\n";
    report << "- OBJ: `" << options.obj_path << "`\n";
    report << "- Seed: `" << options.seed << "`\n";
    report << "- Topology: `" << options.topology << "`\n";
    report << "- Iterations: `" << options.iterations << "`\n";
    report << "- Initial step: `" << options.step << "`\n";
    report << "- Movable OBJ vertices: `[33, 34, 28, 5, 13, 24]`\n";
    report << "- Initial C/I: `" << best_c << "/" << best_i << "`\n\n";

    std::cout << "===================" << std::endl;
    std::cout << "Mode      : repair-local" << std::endl;
    std::cout << "Loaded    : " << options.obj_path << std::endl;
    std::cout << "Seed      : " << options.seed << std::endl;
    std::cout << "Topology  : " << g_topology << std::endl;
    std::cout << "Iters     : " << options.iterations << std::endl;
    std::cout << "Step      : " << options.step << std::endl;
    std::cout << "Initial C/I/loss: " << best_c << "/" << best_i << "/" << best_loss << std::endl;
    std::cout << "===================" << std::endl;

    for (int iter = 1; iter <= options.iterations; ++iter) {
        Verts3D candidate = current;
        const int vertex_ix = options.movable_vertices[vertex_pick(eng)];
        candidate[vertex_ix].x() += normal(eng) * step;
        candidate[vertex_ix].y() += normal(eng) * step;
        candidate[vertex_ix].z() += normal(eng) * step;

        int candidate_c = 0;
        int candidate_i = 0;
        const double candidate_loss = local_repair_objective(
            candidate,
            original,
            options.movable_vertices,
            candidate_c,
            candidate_i
        );

        const bool improves_current = candidate_loss < current_loss;
        const double accept_probability = std::exp((current_loss - candidate_loss) / std::max(1e-6, temperature));
        if (improves_current || uniform(eng) < accept_probability) {
            current = candidate;
            current_loss = candidate_loss;
            current_c = candidate_c;
            current_i = candidate_i;
        }

        const bool better_counts = candidate_c < best_c || (candidate_c == best_c && candidate_i < best_i);
        if (better_counts || (candidate_c == best_c && candidate_i == best_i && candidate_loss < best_loss)) {
            best = candidate;
            best_loss = candidate_loss;
            best_c = candidate_c;
            best_i = candidate_i;
            std::cout << "Best iter " << iter << ": C/I/loss "
                      << best_c << "/" << best_i << "/" << best_loss
                      << " step " << step << std::endl;
            report << "- Best iter `" << iter << "`: C/I/loss `"
                   << best_c << "/" << best_i << "/" << best_loss << "`, step `" << step << "`\n";
        }

        step *= options.beta;
        temperature *= options.beta;

        if (iter % options.report_every == 0) {
            std::cout << "Iter " << iter << ": current "
                      << current_c << "/" << current_i << "/" << current_loss
                      << ", best " << best_c << "/" << best_i << "/" << best_loss
                      << ", step " << step << std::endl;
        }
        if (best_c == 0 && best_i == 0) {
            break;
        }
    }

    Planes best_planes;
    v3ds_to_planes(best, g_polys, best_planes);
    save_sample(options.out_prefix.c_str(), best_planes, best, options.seed, true);
    if (best_c == 0 && best_i == 0) {
        export_obj("runtime/candidates/FOUND_top4_candidate.obj", best, g_polys);
    }
    report << "\n- Final best C/I: `" << best_c << "/" << best_i << "`\n";
    report << "- Final best loss: `" << best_loss << "`\n";
    report << "- Output prefix: `" << options.out_prefix << "`\n\n";

    std::cout << "===================" << std::endl;
    std::cout << "Local repair finished." << std::endl;
    std::cout << "Best C/I: " << best_c << "/" << best_i << std::endl;
    if (best_c == 0 && best_i == 0) {
        std::cout << "FOUND candidate saved to runtime/candidates/FOUND_top4_candidate.obj" << std::endl;
    }
    std::cout << "===================" << std::endl;
    return 0;
}

int main(int argc, char* argv[]) {
    if (argc > 1) {
        if (std::string(argv[1]) == "--help" || std::string(argv[1]) == "-h") {
            print_usage(argv[0]);
            return 0;
        }

        StudyOptions options;
        LocalRepairOptions repair_options;
        bool run_study = false;
        bool run_repair_local = false;
        bool run_hunt_local = false;
        bool run_batch_hunt = false;
        bool run_global_search = false;
        for (int i = 1; i < argc; ++i) {
            const std::string arg = argv[i];
            if (arg == "--study") {
                run_study = true;
                if (i + 1 < argc && std::string(argv[i + 1]).rfind("--", 0) != 0) {
                    options.obj_path = argv[++i];
                }
            } else if (arg == "--repair-local") {
                run_repair_local = true;
                if (i + 1 < argc && std::string(argv[i + 1]).rfind("--", 0) != 0) {
                    repair_options.obj_path = argv[++i];
                }
            } else if (arg == "--hunt-local") {
                run_hunt_local = true;
                if (i + 1 < argc && std::string(argv[i + 1]).rfind("--", 0) != 0) {
                    repair_options.obj_path = argv[++i];
                }
            } else if (arg == "--batch-hunt") {
                run_batch_hunt = true;
                if (i + 1 < argc && std::string(argv[i + 1]).rfind("--", 0) != 0) {
                    repair_options.obj_path = argv[++i];
                }
            } else if (arg == "--global-search") {
                run_global_search = true;
            } else if (arg == "--obj" && i + 1 < argc) {
                options.obj_path = argv[++i];
                repair_options.obj_path = options.obj_path;
            } else if (arg == "--topology" && i + 1 < argc) {
                options.topology = std::atoi(argv[++i]);
                repair_options.topology = options.topology;
            } else if (arg == "--seed" && i + 1 < argc) {
                options.seed = std::atoi(argv[++i]);
                repair_options.seed = options.seed;
            } else if (arg == "--iters" && i + 1 < argc) {
                options.max_iters = std::atoi(argv[++i]);
                repair_options.iterations = options.max_iters;
            } else if (arg == "--clusters" && i + 1 < argc) {
                options.clusters = std::atoi(argv[++i]);
                repair_options.trials = options.clusters;
            } else if (arg == "--sigma" && i + 1 < argc) {
                options.sigma = (double)std::atof(argv[++i]);
                repair_options.step = options.sigma;
            } else if (arg == "--beta" && i + 1 < argc) {
                options.beta = (double)std::atof(argv[++i]);
                repair_options.beta = options.beta;
            } else if (arg == "--objective" && i + 1 < argc) {
                options.objective = argv[++i];
            } else if (arg == "--symmetry") {
                options.use_symmetry = true;
            } else if (arg == "--out" && i + 1 < argc) {
                options.out_prefix = argv[++i];
                repair_options.out_prefix = options.out_prefix;
            } else if (arg == "--report" && i + 1 < argc) {
                options.report_path = argv[++i];
                repair_options.report_path = options.report_path;
            } else if (arg == "--start-planes" && i + 1 < argc) {
                repair_options.start_planes_path = argv[++i];
            } else if (arg == "--temperature" && i + 1 < argc) {
                repair_options.temperature = (double)std::atof(argv[++i]);
            } else if (arg == "--report-every" && i + 1 < argc) {
                repair_options.report_every = std::atoi(argv[++i]);
            } else if (arg == "--trials" && i + 1 < argc) {
                repair_options.trials = std::atoi(argv[++i]);
            } else if (arg == "--threads" && i + 1 < argc) {
                repair_options.threads = std::atoi(argv[++i]);
            } else if (arg == "--minutes" && i + 1 < argc) {
                const double minutes = std::max(0.0, std::atof(argv[++i]));
                repair_options.time_limit_seconds = minutes > 0.0
                    ? static_cast<int>(std::ceil(std::min(
                        minutes * 60.0,
                        static_cast<double>(std::numeric_limits<int>::max()))))
                    : 0;
            } else if (arg == "--stop-file" && i + 1 < argc) {
                repair_options.stop_file_path = argv[++i];
            } else if (arg == "--global-dir" && i + 1 < argc) {
                repair_options.global_dir = argv[++i];
            } else if (arg == "--topology-from" && i + 1 < argc) {
                repair_options.topology_from = std::atoi(argv[++i]);
            } else if (arg == "--topology-to" && i + 1 < argc) {
                repair_options.topology_to = std::atoi(argv[++i]);
            } else if (arg == "--cuda") {
                repair_options.use_cuda = true;
            } else if (arg == "--cuda-chains" && i + 1 < argc) {
                repair_options.cuda_chains = std::max(0, std::atoi(argv[++i]));
            } else if (arg == "--cuda-iters" && i + 1 < argc) {
                repair_options.cuda_iterations = std::clamp(std::atoi(argv[++i]), 1, 64);
            } else if (arg == "--checkpoint-seconds" && i + 1 < argc) {
                repair_options.checkpoint_seconds = std::clamp(std::atoi(argv[++i]), 5, 3600);
            } else if (arg == "--degeneracy-weight" && i + 1 < argc) {
                repair_options.degeneracy_weight = std::clamp(
                    std::atof(argv[++i]), 0.0, 1.0);
            } else if (arg == "--prioritize-worst") {
                repair_options.prioritize_worst = true;
            } else if (arg == "--restarts" && i + 1 < argc) {
                repair_options.restarts = std::atoi(argv[++i]);
            } else if (arg == "--stagnation" && i + 1 < argc) {
                repair_options.stagnation = std::atoi(argv[++i]);
            } else if (arg == "--jump-chance" && i + 1 < argc) {
                repair_options.jump_chance = (double)std::atof(argv[++i]);
            } else if (arg == "--min-step-ratio" && i + 1 < argc) {
                repair_options.min_step_ratio = std::clamp(
                    (double)std::atof(argv[++i]), 1e-12, 1.0);
            } else {
                std::cerr << "Unknown or incomplete argument: " << arg << std::endl;
                print_usage(argv[0]);
                return 2;
            }
        }

        if (run_study) {
            return study_shape(options);
        }
        if (run_global_search) {
            if (repair_options.topology_from < 0 ||
                repair_options.topology_to >= NUM_TOPOLOGIES ||
                repair_options.topology_from > repair_options.topology_to) {
                std::cerr << "Invalid topology range: "
                          << repair_options.topology_from << ".."
                          << repair_options.topology_to << std::endl;
                return 2;
            }
            return global_search_all(repair_options);
        }
        if (run_batch_hunt) {
            return batch_hunt_shape(repair_options);
        }
        if (run_hunt_local) {
            return hunt_local_shape(repair_options);
        }
        if (run_repair_local) {
            return repair_local_shape(repair_options);
        }

        print_usage(argv[0]);
        return 2;
    }

    //validate_files(true);
    //return 0;

    int seed = 123;
    std::cout << "Seed: ";
    std::cin >> seed;
    std::cout << std::endl;
    
    std::cout << "Topology: ";
    std::cin >> g_topology;
    std::cout << std::endl;
    
    set_rand_seed(seed);

    //Run the optimizer (choose one)
    main_solver();
    //quality_solver();
    //explore_shape("results/topologies/topology_42/shape_c0_i4_optimalsymmetric.obj");

    return 0;
}