//#define USE_CAIRO #include "util.h" #include "solver.h" #include "wide_real.h" #include "GlobalCheckpoint.h" #include "SearchArchive.h" #include "OnlineSurrogate.h" #include "TrainingArchive.h" #include "cuda_search.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifdef _WIN32 #ifndef NOMINMAX #define NOMINMAX #endif #include #else #include #include #include #endif #ifdef USE_CAIRO #include #endif #define NUM_TOPOLOGIES 59 #define DUAL_PROBLEM 0 constexpr int GLOBAL_PLANE_VALUE_COUNT = 36; constexpr std::uint32_t GLOBAL_OBJECTIVE_VERSION = 5; constexpr std::uint32_t GLOBAL_NEURAL_SCHEMA_VERSION = 2; constexpr std::uint32_t MIN_COMPATIBLE_ARCHIVE_OBJECTIVE_VERSION = 3; constexpr int CUDA_SESSION_CACHE_LIMIT = NUM_TOPOLOGIES; constexpr std::uint64_t TRAINING_CACHE_LIMIT_BYTES = 200000000000ULL; constexpr std::uint64_t NEURAL_MODEL_BUDGET_PER_TOPOLOGY = 4ULL * 1024ULL * 1024ULL; constexpr std::uint32_t TRAINING_ARCHIVE_SCHEMA_VERSION = 3; 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 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); bool load_global_topology_context(int topology); 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 [options]\n" << " Continue optimization from an existing OBJ candidate.\n\n" << " " << exe_name << " --repair-local [options]\n" << " Move local vertices [33,34,28,5,13,24] with a smooth crossing surrogate.\n\n" << " " << exe_name << " --hunt-local [options]\n" << " Fast focused search around the two known crossings on edge 33-34.\n\n" << " " << exe_name << " --batch-hunt [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 Default: 4\n" << " --seed Default: 30000157\n" << " --iters Default: 360\n" << " --clusters Default: 10000\n" << " --sigma Default: 0.01\n" << " --beta Default: 0.9\n" << " --objective 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 Default: runtime/candidates/top4_study\n" << " --report Default: runtime/reports/02_study_cpp.md\n\n" << "Local repair options:\n" << " --iters Default: 50000\n" << " --sigma Initial coordinate step. Default: 0.5\n" << " --beta Step cooling. Default: 0.9995\n" << " --temperature Annealing temperature. Default: 0.02\n" << " --report-every Default: 1000\n" << " --trials batch-hunt attempts. Default: 2000\n" << " --threads batch-hunt worker threads. 0 = all cores.\n" << " --minutes Clean time limit for batch-hunt. 0 = unlimited.\n" << " --stop-file Stop cleanly when this file appears.\n" << " --global-dir Checkpoints for --global-search. Default: results/search\n" << " --topology-from First topology to search, inclusive. Default: 0\n" << " --topology-to Last topology to search, inclusive. Default: 58\n" << " --cuda Require the FP32 CUDA search backend.\n" << " --cuda-chains Parallel GPU chains. 0 = automatic (default).\n" << " --cuda-iters Iterations per short GPU batch. Default: 64\n" << " --checkpoint-seconds Durable checkpoint period. Default: 30\n" << " --degeneracy-weight Worst-barrier degeneracy penalty. Default: 0.01\n" << " --prioritize-worst Give depth priority to the worst current topologies.\n" << " --start-planes

Continue batch-hunt from a saved .planes sidecar.\n" << " --restarts hunt-local restarts. Default: 256\n" << " --stagnation Iterations before hunt-local reheat. Default: 2000\n" << " --jump-chance Large local jump probability. Default: 0.08\n" << " --min-step-ratio 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 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 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 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 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(0, (int)paths.size() - 1)(eng)]; Verts3D obj_verts; Edges dual_edges; Planes obj_planes; VectorXd obj_x; std::vector 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 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 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 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(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(1, g_edges.size()); } double anchor_loss(const Verts3D& verts, const Verts3D& original, const std::vector& 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(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& movable, int& crossings, int& intersections ) { if (!is_finite(verts)) { return std::numeric_limits::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& known_targets() { static const std::array 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(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 collect_affected_faces(const std::vector& movable) { std::unordered_set movable_set(movable.begin(), movable.end()); std::vector 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(face_ix)); break; } } } return faces; } Edges collect_affected_edges(const std::vector& movable) { std::unordered_set 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& 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(1, affected_edges.size()); } HuntMetrics evaluate_hunt_candidate( const Verts3D& verts, const Verts3D& original, const std::vector& movable, const std::vector& affected_faces, const Edges& affected_edges, double scale, bool verify_full ) { HuntMetrics metrics; if (!is_finite(verts)) { metrics.target_crossings = std::numeric_limits::max() / 4; metrics.crossings = std::numeric_limits::max() / 4; metrics.intersections = std::numeric_limits::max() / 4; metrics.local_loss = std::numeric_limits::infinity(); metrics.total_loss = std::numeric_limits::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& normal) { return Vector3d(normal(eng), normal(eng), normal(eng)); } Vector3d random_vec3(std::normal_distribution& normal, RNG& rng) { return Vector3d(normal(rng), normal(rng), normal(rng)); } void apply_hunt_move( Verts3D& candidate, const std::vector& movable, double step, bool large_jump, std::normal_distribution& normal, std::uniform_real_distribution& uniform ) { const double s = step * (large_jump ? 4.0 : 1.0); const int move_type = std::uniform_int_distribution(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(0, static_cast(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& movable, double step, std::normal_distribution& 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(x.size()); } struct PlaneEvaluationScratch { Planes planes; Planes canonical_planes; Verts3D verts; }; HuntMetrics invalid_hunt_metrics() { HuntMetrics metrics; metrics.target_crossings = std::numeric_limits::max() / 4; metrics.crossings = std::numeric_limits::max() / 4; metrics.intersections = std::numeric_limits::max() / 4; metrics.target_loss = std::numeric_limits::infinity(); metrics.crossing_loss = std::numeric_limits::infinity(); metrics.regularization_loss = std::numeric_limits::infinity(); metrics.local_loss = std::numeric_limits::infinity(); metrics.total_loss = std::numeric_limits::infinity(); return metrics; } bool valid_plane_state(const VectorXd& x) { if (x.size() != static_cast(g_polys.size() * 3) || !x.allFinite()) { return false; } for (int i = 0; i < x.size(); i += 3) { if (Eigen::Map(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(defects) + 0.01 * static_cast(metrics.crossings) + 0.0125 * static_cast(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& movable_faces, double step, bool large_jump, std::normal_distribution& normal ) { const double s = step * (large_jump ? 4.0 : 1.0); const int move_type = std::uniform_int_distribution(0, 5)(eng); auto add_plane = [&](int face_ix, const Vector3d& delta) { Eigen::Map(x.data() + face_ix * 3) += delta; }; if (move_type == 0) { const int face_ix = movable_faces[std::uniform_int_distribution(0, static_cast(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& movable_faces, double step, bool large_jump, std::normal_distribution& normal, RNG& rng ) { const double s = step * (large_jump ? 4.0 : 1.0); const int move_type = std::uniform_int_distribution(0, 5)(rng); auto add_plane = [&](int face_ix, const Vector3d& delta) { Eigen::Map(x.data() + face_ix * 3) += delta; }; if (move_type == 0) { const int face_ix = movable_faces[std::uniform_int_distribution(0, static_cast(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& movable_faces, double step, std::normal_distribution& normal ) { for (int face_ix : movable_faces) { Eigen::Map(x.data() + face_ix * 3) += random_vec3(normal) * (step * 0.35); } } void jitter_seed_planes( VectorXd& x, const std::vector& movable_faces, double step, std::normal_distribution& normal, RNG& rng ) { for (int face_ix : movable_faces) { Eigen::Map(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 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 normal(0.0, 1.0); std::uniform_real_distribution 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 unique_faces(std::vector faces) { std::sort(faces.begin(), faces.end()); faces.erase(std::unique(faces.begin(), faces.end()), faces.end()); return faces; } std::vector all_face_indices() { std::vector faces; faces.reserve(g_polys.size()); for (int i = 0; i < static_cast(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 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(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(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(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(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::max() / 4; intersections = std::numeric_limits::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& movable_faces, int trial_ix, const std::atomic& 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(movable_faces.size()); RNG trial_rng(trial_options.seed); std::normal_distribution normal(0.0, 1.0); std::uniform_real_distribution 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 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 target_faces = unique_faces({4, 5}); std::vector 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(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& 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( static_cast(options.seed) + static_cast(trial * 104729 + strategy * 7919) )); std::uniform_real_distribution uniform(0.0, 1.0); std::uniform_int_distribution seed_noise(0, 999999); LocalRepairOptions trial_options = options; const long long seed_value = static_cast(options.seed) + 1000003LL * trial + 7919LL * strategy + seed_noise(param_rng); trial_options.seed = static_cast( 1 + static_cast(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(uniform(param_rng) * options.stagnation)); trial_options.restarts = std::max(1, options.restarts / 4); return trial_options; }; std::atomic next_trial{1}; std::atomic completed_trials{0}; std::atomic active_workers{worker_count}; std::atomic found{false}; std::atomic stop_requested{false}; std::mutex result_mutex; std::vector 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 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 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 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(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 lock(result_mutex); std::cout << "Progress : " << completed << "/" << options.trials << " complete, " << issued - completed << " active, " << std::fixed << std::setprecision(2) << (static_cast(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::max() / 4; int intersections = std::numeric_limits::max() / 4; double crossing_loss = std::numeric_limits::infinity(); double intersection_loss = std::numeric_limits::infinity(); double geometry_penalty = std::numeric_limits::infinity(); double degeneracy_penalty = std::numeric_limits::infinity(); double worst_degeneracy = std::numeric_limits::infinity(); double condition_number = std::numeric_limits::infinity(); double condition_penalty = std::numeric_limits::infinity(); double energy = std::numeric_limits::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::infinity(); std::uint16_t crossing_face_mask = 0; std::uint16_t intersection_face_mask = 0; 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(static_cast(value)); hash *= 1099511628211ULL; }; for (const Face& triangle : g_tris) { add(static_cast(triangle.size())); for (int value : triangle) { add(value); } } for (const Face& polygon : g_polys) { add(static_cast(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 archive; std::unordered_set archive_dirty; struct DiagonalCemState { bool initialized = false; std::array mean{}; std::array variance{}; std::uint64_t updates = 0; } cem; std::unordered_set cem_seen_state_hashes; std::unique_ptr neural; bool neural_dirty = false; std::uint64_t neural_sequence = 0; // 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; // Run-local strategy portfolio. Durable geometry remains independent of // these values, so two machines can merge ordinary checkpoints without // contending for a learned scheduler state. std::array strategy_reward_ema{{0.15, 0.02, 0.05, 0.03, 0.25}}; std::array strategy_observations{}; 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 strategy_evaluated{}; std::array strategy_verified{}; std::array strategy_archive_improvements{}; std::array strategy_global_improvements{}; std::array strategy_weights{{4, 3, 3, 3, 3}}; std::uint64_t neural_samples_added = 0; std::uint64_t neural_training_steps = 0; std::uint64_t neural_seed_count = 0; std::uint64_t training_verified_records = 0; std::uint64_t training_fp32_records = 0; std::uint64_t training_seed_records = 0; std::uint64_t training_rollout_records = 0; std::uint64_t training_refinement_records = 0; } last_round; }; struct GlobalTrialResult { GlobalMetrics best; VectorXd best_x; int iterations = 0; std::uint64_t evaluated_states = 0; }; struct GlobalVerifiedCandidate { GlobalMetrics metrics; VectorXd x; std::uint32_t strategy = 0; std::uint64_t chain_id = 0; }; struct TrainingArchiveSession { std::unique_ptr writer; szilassi::training::RunId run_id{}; std::filesystem::path run_directory; std::uint64_t next_record_sequence = 0; std::uint64_t legacy_records = 0; bool fatal_error = false; bool limit_warning_emitted = false; }; std::uint64_t training_hash_bytes(const void* data, std::size_t size) { const auto* bytes = static_cast(data); std::uint64_t hash = 1469598103934665603ULL; for (std::size_t index = 0; index < size; ++index) { hash ^= static_cast(bytes[index]); hash *= 1099511628211ULL; } return hash; } std::uint64_t training_hash_text(const std::string& text) { return training_hash_bytes(text.data(), text.size()); } std::uint64_t training_splitmix64(std::uint64_t value) { value += 0x9e3779b97f4a7c15ULL; value = (value ^ (value >> 30U)) * 0xbf58476d1ce4e5b9ULL; value = (value ^ (value >> 27U)) * 0x94d049bb133111ebULL; return value ^ (value >> 31U); } std::uint64_t cuda_chain_initial_rng_state( std::uint64_t session_seed, std::uint64_t chain_id ) { std::uint64_t state = training_splitmix64( session_seed ^ (chain_id + 1ULL) * 0x9e3779b97f4a7c15ULL); if (state == 0) state = 0x2545f4914f6cdd1dULL; return state; } szilassi::training::RunId training_run_id_from_string(const std::string& text) { szilassi::training::RunId result; std::string hexadecimal; hexadecimal.reserve(32); for (char value : text) { if (value != '-') { hexadecimal.push_back(value); } } auto digit = [](char value) -> int { if (value >= '0' && value <= '9') return value - '0'; if (value >= 'a' && value <= 'f') return 10 + value - 'a'; if (value >= 'A' && value <= 'F') return 10 + value - 'A'; return -1; }; bool valid = hexadecimal.size() == result.bytes.size() * 2; for (std::size_t index = 0; valid && index < result.bytes.size(); ++index) { const int high = digit(hexadecimal[index * 2]); const int low = digit(hexadecimal[index * 2 + 1]); valid = high >= 0 && low >= 0; if (valid) { result.bytes[index] = static_cast((high << 4) | low); } } if (!valid) { const std::uint64_t first = training_hash_text(text); const std::uint64_t second = training_hash_text("szilassi-run:" + text); std::memcpy(result.bytes.data(), &first, sizeof(first)); std::memcpy(result.bytes.data() + sizeof(first), &second, sizeof(second)); } return result; } szilassi::training::PlaneState training_plane_state( const cuda_search::PlaneState& source ) { szilassi::training::PlaneState result; result.values = source.values; return result; } szilassi::training::PlaneState training_plane_state(const VectorXd& source) { szilassi::training::PlaneState result; for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) { result.values[static_cast(component)] = component < source.size() && std::isfinite(source[component]) ? static_cast(source[component]) : 0.0f; } return result; } szilassi::training::PrecisePlaneState training_precise_plane_state( const VectorXd& source ) { szilassi::training::PrecisePlaneState result; for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) { result.values[static_cast(component)] = component < source.size() && std::isfinite(source[component]) ? source[component] : 0.0; } return result; } std::uint64_t training_plane_hash(const szilassi::training::PlaneState& state) { std::array bytes{}; for (std::size_t component = 0; component < state.values.size(); ++component) { std::uint32_t bits = 0; std::memcpy(&bits, &state.values[component], sizeof(bits)); for (std::size_t byte = 0; byte < 4; ++byte) { bytes[component * 4 + byte] = static_cast( bits >> (byte * 8)); } } return training_hash_bytes(bytes.data(), bytes.size()); } szilassi::training::Metrics training_metrics(const GlobalMetrics& source) { szilassi::training::Metrics result; result.crossings = source.crossings; result.intersections = source.intersections; result.crossing_loss = source.crossing_loss; result.intersection_loss = source.intersection_loss; result.geometry_penalty = source.geometry_penalty; result.degeneracy_penalty = source.degeneracy_penalty; result.worst_degeneracy = source.worst_degeneracy; result.energy = source.energy; result.min_plane_determinant = source.min_plane_determinant; result.relative_min_edge = source.relative_min_edge; result.min_turn_sine = source.min_turn_sine; result.max_vertex_norm = source.max_vertex_norm; result.condition_number = source.condition_number; result.condition_penalty = source.condition_penalty; result.crossing_face_mask = source.crossing_face_mask; result.intersection_face_mask = source.intersection_face_mask; result.flags = szilassi::training::MetricFinite | szilassi::training::MetricHasSmoothLosses | szilassi::training::MetricHasGeometryPenalty | szilassi::training::MetricHasWorstDegeneracy | szilassi::training::MetricHasShapeDescriptors | szilassi::training::MetricHasCondition | szilassi::training::MetricHasFaceMasks | szilassi::training::MetricHasEnergy; if (source.canonical) { result.flags |= szilassi::training::MetricCanonical; } if (source.precise) { result.flags |= szilassi::training::MetricPrecise | szilassi::training::MetricDdVerified; } return result; } szilassi::training::ApproximateMetrics training_cuda_metrics( const cuda_search::Candidate& source ) { szilassi::training::ApproximateMetrics result; result.crossings = std::max(0, source.crossings); result.intersections = std::max(0, source.intersections); result.flags = szilassi::training::ApproximateHasCounts; if (std::isfinite(source.crossing_loss) && source.crossing_loss >= 0.0f && std::isfinite(source.intersection_loss) && source.intersection_loss >= 0.0f) { result.crossing_loss = source.crossing_loss; result.intersection_loss = source.intersection_loss; result.flags |= szilassi::training::ApproximateHasCudaLosses; } if (std::isfinite(source.geometry_penalty) && source.geometry_penalty >= 0.0f && std::isfinite(source.degeneracy_penalty) && source.degeneracy_penalty >= 0.0f) { result.geometry_penalty = source.geometry_penalty; result.degeneracy_penalty = source.degeneracy_penalty; result.flags |= szilassi::training::ApproximateHasGeometryPenalty; } if (std::isfinite(source.energy)) { result.energy = source.energy; result.flags |= szilassi::training::ApproximateHasEnergy; } if (std::isfinite(source.min_plane_determinant) && source.min_plane_determinant >= 0.0f && std::isfinite(source.relative_min_edge) && source.relative_min_edge >= 0.0f && std::isfinite(source.min_turn_sine) && source.min_turn_sine >= 0.0f && std::isfinite(source.max_vertex_norm) && source.max_vertex_norm >= 0.0f) { result.min_plane_determinant = source.min_plane_determinant; result.relative_min_edge = source.relative_min_edge; result.min_turn_sine = source.min_turn_sine; result.max_vertex_norm = source.max_vertex_norm; result.flags |= szilassi::training::ApproximateHasShapeDescriptors; } result.ambiguity_flags = source.ambiguity_flags; result.flags |= szilassi::training::ApproximateHasAmbiguity; return result; } szilassi::training::BehaviorPrediction training_behavior_prediction( const szilassi::surrogate::Prediction& source ) { szilassi::training::BehaviorPrediction result; result.improvement_logit = source.improvement_logit; result.improvement_probability = source.improvement_probability; result.expected_defect_gain = source.expected_defect_gain; result.uncertainty = source.uncertainty; result.plane_probabilities = source.plane_probabilities; result.move_probabilities = source.move_probabilities; result.scale_probabilities = source.scale_probabilities; result.flags = source.finite ? 1U : 0U; return result; } std::string training_status_system_error(const char* operation, int code) { #ifdef _WIN32 return std::string(operation) + ": " + std::system_category().message(code); #else return std::string(operation) + ": " + std::generic_category().message(code); #endif } void set_training_status_error(std::string* error, std::string message) { if (error != nullptr) { *error = std::move(message); } } bool durably_flush_file( const std::filesystem::path& path, std::string& error ) { #ifdef _WIN32 const HANDLE handle = CreateFileW( 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) { error = training_status_system_error( "Cannot open training status file for durable flush", static_cast(GetLastError())); return false; } const bool flushed = FlushFileBuffers(handle) != 0; const int flush_error = flushed ? 0 : static_cast(GetLastError()); const bool closed = CloseHandle(handle) != 0; const int close_error = closed ? 0 : static_cast(GetLastError()); if (!flushed) { error = training_status_system_error( "Cannot flush training status file", flush_error); return false; } if (!closed) { error = training_status_system_error( "Cannot close training status file", close_error); return false; } return true; #else int flags = O_RDONLY; #ifdef O_CLOEXEC flags |= O_CLOEXEC; #endif const int descriptor = ::open(path.c_str(), flags); if (descriptor < 0) { error = training_status_system_error( "Cannot open training status file for fsync", errno); return false; } const bool flushed = ::fsync(descriptor) == 0; const int flush_error = flushed ? 0 : errno; const bool closed = ::close(descriptor) == 0; const int close_error = closed ? 0 : errno; if (!flushed) { error = training_status_system_error( "Cannot fsync training status file", flush_error); return false; } if (!closed) { error = training_status_system_error( "Cannot close training status file", close_error); return false; } return true; #endif } bool durably_flush_parent_directory( const std::filesystem::path& path, std::string& error ) { #ifdef _WIN32 (void)path; (void)error; // MOVEFILE_WRITE_THROUGH below is the Windows publication boundary. return true; #else const std::filesystem::path directory = path.parent_path().empty() ? std::filesystem::path(".") : path.parent_path(); int flags = O_RDONLY; #ifdef O_CLOEXEC flags |= O_CLOEXEC; #endif #ifdef O_DIRECTORY flags |= O_DIRECTORY; #endif const int descriptor = ::open(directory.c_str(), flags); if (descriptor < 0) { error = training_status_system_error( "Cannot open training status directory for fsync", errno); return false; } const bool flushed = ::fsync(descriptor) == 0; const int flush_error = flushed ? 0 : errno; const bool closed = ::close(descriptor) == 0; const int close_error = closed ? 0 : errno; if (!flushed) { error = training_status_system_error( "Cannot fsync training status directory", flush_error); return false; } if (!closed) { error = training_status_system_error( "Cannot close training status directory", close_error); return false; } return true; #endif } bool write_training_status_file( const std::filesystem::path& path, const std::string& contents, std::string& error ) { std::ofstream out(path, std::ios::binary | std::ios::trunc); if (!out) { error = "Cannot create training status temporary file: " + path.string(); return false; } out << contents; out.flush(); if (!out) { error = "Cannot write training status temporary file: " + path.string(); return false; } out.close(); if (!out) { error = "Cannot close training status temporary file: " + path.string(); return false; } return durably_flush_file(path, error); } bool write_training_status_file_atomic( const std::filesystem::path& path, const std::string& contents, std::string* error = nullptr ) { set_training_status_error(error, {}); std::error_code directory_error; if (path.has_parent_path()) { std::filesystem::create_directories(path.parent_path(), directory_error); if (directory_error) { set_training_status_error( error, "Cannot create training status directory: " + directory_error.message()); return false; } } const std::uint64_t nonce = static_cast( std::chrono::high_resolution_clock::now().time_since_epoch().count()); std::ostringstream suffix; #ifdef _WIN32 suffix << ".tmp." << GetCurrentProcessId() << "." << std::hex << nonce; #else suffix << ".tmp." << ::getpid() << "." << std::hex << nonce; #endif const std::filesystem::path temporary = path.string() + suffix.str(); std::string operation_error; if (!write_training_status_file(temporary, contents, operation_error)) { std::error_code remove_error; std::filesystem::remove(temporary, remove_error); set_training_status_error(error, std::move(operation_error)); return false; } #ifdef _WIN32 const bool published = MoveFileExW( temporary.c_str(), path.c_str(), MOVEFILE_REPLACE_EXISTING | MOVEFILE_WRITE_THROUGH) != 0; if (!published) { const int publish_error = static_cast(GetLastError()); std::error_code remove_error; std::filesystem::remove(temporary, remove_error); set_training_status_error( error, training_status_system_error( "Cannot atomically publish training status file", publish_error)); return false; } #else if (::rename(temporary.c_str(), path.c_str()) != 0) { const int publish_error = errno; std::error_code remove_error; std::filesystem::remove(temporary, remove_error); set_training_status_error( error, training_status_system_error( "Cannot atomically publish training status file", publish_error)); return false; } if (!durably_flush_parent_directory(path, operation_error)) { // The rename may already be visible. Do not delete the only complete // marker; report that its directory entry could not be guaranteed // durable so the caller can retry or stop explicitly. set_training_status_error(error, std::move(operation_error)); return false; } #endif return true; } void poll_training_archive_notices(TrainingArchiveSession& archive) { if (archive.writer == nullptr) return; std::string detail; const bool transitioned = archive.writer->take_limit_notice(detail); if ((transitioned || archive.writer->limit_reached()) && !archive.limit_warning_emitted) { archive.limit_warning_emitted = true; const std::string warning = "ВНИМАНИЕ: достигнут лимит обучающего кэша 200 000 000 000 байт. " "Сбор новых данных и онлайн-обучение остановлены; геометрический поиск " "продолжается с уже обученной моделью. Перед дальнейшим накоплением " "нужно пересмотреть программу и формат архива."; std::cerr << warning << std::endl; std::ostringstream marker; marker << "format\tszilassi-training-limit-v1\n" << "limit_bytes\t" << TRAINING_CACHE_LIMIT_BYTES << "\n" << "accounted_bytes\t" << archive.writer->committed_bytes() << "\n" << "collection_enabled\t0\n" << "required_action\trewrite-or-revise-training-archive\n" << "detail\t" << detail << "\n"; const std::filesystem::path marker_path = archive.run_directory / "TRAINING_DATA_LIMIT_REACHED_REWRITE_REQUIRED.tsv"; std::string marker_error; if (!write_training_status_file_atomic( marker_path, marker.str(), &marker_error)) { std::cerr << "Cannot durably write training-limit warning marker: " << marker_path << ": " << marker_error << std::endl; } } detail.clear(); if (archive.writer->take_write_error_notice(detail)) { archive.fatal_error = true; std::cerr << "Training archive write failed; stopping search safely: " << detail << std::endl; } } template bool append_training_record(TrainingArchiveSession* archive, const Record& record) { if (archive == nullptr || archive->writer == nullptr || archive->fatal_error || !archive->writer->collection_enabled()) { return false; } const szilassi::training::AppendResult result = archive->writer->append(record); poll_training_archive_notices(*archive); if (result.status == szilassi::training::AppendStatus::InvalidRecord) { archive->fatal_error = true; std::cerr << "Internal error: rejected invalid training archive record; " << "stopping instead of losing training data." << std::endl; } else if (result.status == szilassi::training::AppendStatus::WriteError) { archive->fatal_error = true; } return result.accepted; } bool checkpoint_training_archive(TrainingArchiveSession* archive) { if (archive == nullptr || archive->writer == nullptr) return true; std::string error; if (!archive->writer->flush(&error)) { archive->fatal_error = true; std::cerr << "Cannot durably checkpoint training archive: " << error << std::endl; poll_training_archive_notices(*archive); return false; } poll_training_archive_notices(*archive); return true; } bool seal_training_archive(TrainingArchiveSession* archive) { if (archive == nullptr || archive->writer == nullptr) return true; std::string error; if (!archive->writer->seal(&error)) { archive->fatal_error = true; std::cerr << "Cannot seal immutable training shard: " << error << std::endl; poll_training_archive_notices(*archive); return false; } poll_training_archive_notices(*archive); return true; } szilassi::training::CommonContext make_training_context( TrainingArchiveSession& archive, const LocalRepairOptions& options, const GlobalTopologyState& state, std::uint64_t round_seed, const szilassi::training::BatchConfig& batch, std::uint64_t session_seed = 0, std::uint64_t session_batch_ordinal = 0, std::uint64_t completed_trials_offset = 0, std::uint64_t completed_iterations_offset = 0 ) { szilassi::training::CommonContext context; context.run_id = archive.run_id; context.record_sequence = ++archive.next_record_sequence; context.unix_time_ns = szilassi::checkpoint::unix_time_ns_now(); context.topology_fingerprint = topology_fingerprint(); context.algorithm_fingerprint = training_hash_text( options.use_cuda ? "hybrid-quality-diversity-neural-v3-cuda-fp32" : "cpu-parallel-simulated-annealing-double"); static const std::uint64_t build_fingerprint = training_hash_text( std::string(__DATE__) + " " + __TIME__); context.build_fingerprint = build_fingerprint; context.device_fingerprint = training_hash_text(g_search_device_id); context.run_seed = static_cast( static_cast(options.seed)); context.round_seed = round_seed; context.session_seed = session_seed; context.session_batch_ordinal = session_batch_ordinal; context.completed_trials = state.trials + completed_trials_offset; context.completed_iterations = state.iterations + completed_iterations_offset; context.topology = static_cast(state.topology); context.objective_version = GLOBAL_OBJECTIVE_VERSION; context.search_mode = options.use_cuda ? 1U : 2U; context.flags = TRAINING_ARCHIVE_SCHEMA_VERSION; context.batch = batch; return context; } szilassi::training::BatchConfig make_cuda_training_batch_config( const LocalRepairOptions& options, int effective_chain_count, const std::array& strategy_weights, const cuda_search::BatchRunConfig& run_config, const cuda_search::BatchResult& result ) { szilassi::training::BatchConfig batch; batch.chain_count = static_cast(std::max(0, effective_chain_count)); batch.iterations_per_chain = static_cast( std::max(0, options.cuda_iterations)); batch.iterations_per_kernel = static_cast( std::clamp(options.cuda_iterations, 1, 16)); batch.shortlist_size = static_cast( std::min(128, std::max(0, effective_chain_count))); const int total_weight = std::accumulate( strategy_weights.begin(), strategy_weights.end(), 0); std::array counts{}; if (total_weight > 0) { for (int chain = 0; chain < effective_chain_count; ++chain) { const int residue = chain % total_weight; int boundary = 0; for (std::size_t strategy = 0; strategy < counts.size(); ++strategy) { boundary += strategy_weights[strategy]; if (residue < boundary) { counts[strategy] += 1; break; } } } } batch.baseline_chains = counts[0]; batch.replica_chains = counts[1]; batch.adaptive_chains = counts[2]; batch.pbt_chains = counts[3]; // This field describes the disjoint strategy cohort, not the number of // host seeds actually copied into that cohort during this batch. batch.injected_chains = counts[4]; batch.fresh_numerator = static_cast( std::max(0, run_config.fresh_numerator)); batch.fresh_denominator = static_cast( std::max(1, run_config.fresh_denominator)); batch.replica_group_size = static_cast( std::max(0, run_config.replica_group_size)); batch.injected_pool_size = static_cast( std::min( run_config.injected_states.size(), std::numeric_limits::max())); batch.evaluated_states = result.evaluated_states; batch.initialization_evaluated_states = result.initialization_evaluated_states; batch.additional_evaluated_states = result.strategy_evaluated_states; batch.stagnation_iterations = static_cast( std::max(256, options.stagnation)); batch.initial_temperature = static_cast(options.temperature); batch.final_temperature = static_cast(options.temperature * 0.04); batch.proposal_scale = static_cast(options.step) * run_config.step_scale; batch.minimum_step = static_cast( std::max(1e-7, options.step * options.min_step_ratio)); batch.cooling = static_cast(options.beta); batch.minimum_temperature = static_cast(options.temperature * 0.04); batch.jump_chance = static_cast(options.jump_chance); batch.initial_state_jitter = static_cast(options.step * 0.10); batch.replica_temperature_ratio = run_config.replica_temperature_ratio; batch.pbt_exploit_chance = run_config.pbt_exploit_chance; batch.pbt_state_jitter = run_config.pbt_state_jitter; batch.injected_state_jitter = run_config.injected_state_jitter; batch.degeneracy_weight = static_cast(options.degeneracy_weight); for (std::size_t strategy = 0; strategy < strategy_weights.size(); ++strategy) { // Preserve the exact integer cohort pattern. Counts above provide the // realized fractions on a particular GPU, so normalization later is // lossless while the chain-id mapping remains reproducible. batch.strategy_weights[strategy] = static_cast( strategy_weights[strategy]); } return batch; } szilassi::training::BatchConfig make_cpu_training_batch_config( const LocalRepairOptions& options, int trial_count, int iterations ) { szilassi::training::BatchConfig batch; batch.chain_count = static_cast(std::max(0, trial_count)); batch.iterations_per_chain = static_cast(std::max(0, iterations)); batch.iterations_per_kernel = 1; batch.shortlist_size = static_cast(std::max(0, trial_count)); batch.baseline_chains = batch.chain_count; batch.evaluated_states = static_cast(batch.chain_count) * static_cast(batch.iterations_per_chain); batch.fresh_denominator = 1; batch.stagnation_iterations = static_cast( std::max(0, options.stagnation)); batch.initial_temperature = static_cast(options.temperature); batch.final_temperature = static_cast( options.temperature * std::pow( options.beta, static_cast(std::max(0, iterations)))); batch.proposal_scale = static_cast(options.step); batch.minimum_step = static_cast( std::max(1e-12, options.step * options.min_step_ratio)); batch.cooling = static_cast(options.beta); batch.minimum_temperature = 0.0f; batch.jump_chance = static_cast(options.jump_chance); batch.degeneracy_weight = static_cast(options.degeneracy_weight); batch.strategy_weights[0] = 1.0f; return batch; } std::filesystem::path neural_model_relative_path(int topology) { return std::filesystem::path( "topology_" + std::to_string(topology)) / "neural" / ("online_surrogate_v" + std::to_string(GLOBAL_OBJECTIVE_VERSION) + "_s" + std::to_string(GLOBAL_NEURAL_SCHEMA_VERSION) + ".szonn"); } bool scan_neural_cache_bytes( const std::filesystem::path& root, std::uint64_t& bytes, std::array& final_snapshot_bytes, std::string& error ) { bytes = 0; final_snapshot_bytes.fill(0); std::error_code exists_error; const bool exists = std::filesystem::exists(root, exists_error); if (exists_error) { error = "Cannot inspect cache root " + root.string() + ": " + exists_error.message(); return false; } if (!exists) return true; std::error_code iterator_error; std::filesystem::recursive_directory_iterator iterator( root, std::filesystem::directory_options::none, iterator_error); const std::filesystem::recursive_directory_iterator end; while (!iterator_error && iterator != end) { std::error_code status_error; const bool regular = iterator->is_regular_file(status_error); if (status_error) { error = "Cannot inspect cache entry " + iterator->path().string() + ": " + status_error.message(); return false; } bool inside_neural_directory = false; const std::filesystem::path relative = iterator->path().lexically_relative(root); for (const auto& component : relative) { if (component == "neural") { inside_neural_directory = true; break; } } if (regular && inside_neural_directory) { const std::uint64_t size = iterator->file_size(status_error); if (status_error || size > std::numeric_limits::max() - bytes) { error = "Cannot account cache file " + iterator->path().string(); return false; } bytes += size; for (int topology = 0; topology < NUM_TOPOLOGIES; ++topology) { if (relative == neural_model_relative_path(topology)) { final_snapshot_bytes[static_cast(topology)] = size; break; } } } iterator.increment(iterator_error); } if (iterator_error) { error = "Cannot scan cache root " + root.string() + ": " + iterator_error.message(); return false; } return true; } bool scan_training_archive_tree_bytes( const std::filesystem::path& runs_root, std::uint64_t& bytes, std::string& error ) { bytes = 0; std::error_code exists_error; if (!std::filesystem::exists(runs_root, exists_error)) { if (exists_error) { error = "Cannot inspect run archive root: " + exists_error.message(); return false; } return true; } std::error_code iterator_error; std::filesystem::directory_iterator iterator(runs_root, iterator_error); const std::filesystem::directory_iterator end; while (!iterator_error && iterator != end) { std::error_code type_error; if (iterator->is_directory(type_error) && !type_error) { const std::filesystem::path training_directory = iterator->path() / "training"; std::uint64_t run_bytes = 0; std::string scan_error; if (!szilassi::training::scan_archive_bytes( training_directory, run_bytes, &scan_error)) { error = std::move(scan_error); return false; } if (run_bytes > std::numeric_limits::max() - bytes) { error = "Training archive byte count overflow"; return false; } bytes += run_bytes; } else if (type_error) { error = "Cannot inspect run archive entry: " + type_error.message(); return false; } iterator.increment(iterator_error); } if (iterator_error) { error = "Cannot enumerate run archive root: " + iterator_error.message(); return false; } return true; } bool recover_training_archive_tree( const std::filesystem::path& runs_root, std::uint64_t archive_cap_bytes, szilassi::training::RecoveryReport& aggregate, std::string& error ) { aggregate = {}; std::uint64_t total_bytes = 0; if (!scan_training_archive_tree_bytes(runs_root, total_bytes, error)) { return false; } std::error_code exists_error; if (!std::filesystem::exists(runs_root, exists_error)) { if (exists_error) { error = "Cannot inspect run archive root: " + exists_error.message(); return false; } return true; } std::vector training_directories; std::error_code iterator_error; std::filesystem::directory_iterator iterator(runs_root, iterator_error); const std::filesystem::directory_iterator end; while (!iterator_error && iterator != end) { std::error_code type_error; if (iterator->is_directory(type_error) && !type_error) { training_directories.push_back(iterator->path() / "training"); } else if (type_error) { error = "Cannot inspect run archive entry: " + type_error.message(); return false; } iterator.increment(iterator_error); } if (iterator_error) { error = "Cannot enumerate run archive root: " + iterator_error.message(); return false; } std::sort(training_directories.begin(), training_directories.end()); for (const std::filesystem::path& directory : training_directories) { const std::filesystem::path wal = directory / "active.sztd.wal"; const std::filesystem::path commit = directory / "active.sztd.wal.commit"; const std::filesystem::path commit_temporary = directory / "active.sztd.wal.commit.tmp"; std::error_code wal_error; const bool has_wal = std::filesystem::exists(wal, wal_error); if (wal_error) { error = "Cannot inspect training WAL " + wal.string() + ": " + wal_error.message(); return false; } std::error_code commit_error; std::error_code commit_temporary_error; const bool has_commit = std::filesystem::exists(commit, commit_error); const bool has_commit_temporary = std::filesystem::exists(commit_temporary, commit_temporary_error); if (commit_error || commit_temporary_error) { error = "Cannot inspect training commit marker in " + directory.string() + ": " + (commit_error ? commit_error.message() : commit_temporary_error.message()); return false; } if (!has_wal && !has_commit && !has_commit_temporary) continue; if (total_bytes > archive_cap_bytes && has_wal) { // A Git merge can combine two independently capped clones. Keep // the local WAL byte-for-byte for a future archive revision, but // do not make an already-over-limit cache stop geometric search. ++aggregate.wals_found; ++aggregate.wals_deferred_at_limit; continue; } std::uint64_t local_before = 0; std::string local_error; if (!szilassi::training::scan_archive_bytes( directory, local_before, &local_error)) { error = std::move(local_error); return false; } const std::uint64_t remaining = total_bytes <= archive_cap_bytes ? archive_cap_bytes - total_bytes : 0; if (local_before > std::numeric_limits::max() - remaining) { error = "Training recovery byte-count overflow"; return false; } szilassi::training::RecoveryReport report; if (!szilassi::training::recover_archive_wals( directory, local_before + remaining, report, &local_error)) { error = std::move(local_error); return false; } std::uint64_t local_after = 0; if (!szilassi::training::scan_archive_bytes( directory, local_after, &local_error)) { error = std::move(local_error); return false; } if (local_after >= local_before) { const std::uint64_t added = local_after - local_before; if (added > archive_cap_bytes - total_bytes) { error = "Recovered training WAL exceeded the exact cache allocation"; return false; } total_bytes += added; } else { const std::uint64_t removed = local_before - local_after; if (removed > total_bytes) { error = "Recovered training WAL produced an invalid byte count"; return false; } total_bytes -= removed; } aggregate.wals_found += report.wals_found; aggregate.shards_sealed += report.shards_sealed; aggregate.wals_deferred_at_limit += report.wals_deferred_at_limit; aggregate.empty_wals_removed += report.empty_wals_removed; aggregate.torn_bytes_discarded += report.torn_bytes_discarded; aggregate.records_recovered += report.records_recovered; } return true; } std::filesystem::path legacy_replay_export_marker( const LocalRepairOptions& options, int topology ) { return std::filesystem::path(options.global_dir) / ("topology_" + std::to_string(topology)) / "neural" / "long_term_archive_v1_exported.tsv"; } struct LegacyReplayExportState { bool exists = false; bool complete = false; std::uint64_t expected_records = 0; std::vector run_ids; }; std::string serialize_legacy_replay_export_state( int topology, const LegacyReplayExportState& state ) { std::ostringstream marker; marker << "format\tszilassi-legacy-replay-export-v2\n" << "status\t" << (state.complete ? "complete" : "pending") << "\n" << "topology\t" << topology << "\n" << "objective_version\t" << GLOBAL_OBJECTIVE_VERSION << "\n" << "neural_schema\t" << GLOBAL_NEURAL_SCHEMA_VERSION << "\n" << "expected_records\t" << state.expected_records << "\n"; for (const std::string& run_id : state.run_ids) { marker << "run_id\t" << run_id << "\n"; } return marker.str(); } bool read_legacy_replay_export_state( const std::filesystem::path& path, int expected_topology, LegacyReplayExportState& state, std::string& error ) { state = {}; std::error_code exists_error; if (!std::filesystem::exists(path, exists_error)) { if (exists_error) { error = "Cannot inspect legacy replay marker: " + exists_error.message(); return false; } return true; } std::ifstream input(path, std::ios::binary); if (!input) { error = "Cannot open legacy replay marker: " + path.string(); return false; } state.exists = true; int topology = -1; int objective_version = -1; int neural_schema = -1; bool has_expected = false; bool old_complete_marker = false; std::string line; while (std::getline(input, line)) { if (!line.empty() && line.back() == '\r') line.pop_back(); const std::size_t separator = line.find('\t'); if (separator == std::string::npos) continue; const std::string key = line.substr(0, separator); const std::string value = line.substr(separator + 1); try { if (key == "status") { state.complete = value == "complete"; } else if (key == "topology") { topology = std::stoi(value); } else if (key == "objective_version") { objective_version = std::stoi(value); } else if (key == "neural_schema") { neural_schema = std::stoi(value); } else if (key == "expected_records") { state.expected_records = std::stoull(value); has_expected = true; } else if (key == "records") { state.expected_records = std::stoull(value); has_expected = true; old_complete_marker = true; } else if (key == "run_id") { if (value.empty() || value.size() > 64 || value.find_first_of("/\\") != std::string::npos) { error = "Legacy replay marker contains an invalid run UUID"; return false; } if (std::find(state.run_ids.begin(), state.run_ids.end(), value) == state.run_ids.end()) { state.run_ids.push_back(value); } } } catch (const std::exception&) { error = "Legacy replay marker contains an invalid numeric field"; return false; } } if (!input.eof() || topology != expected_topology || objective_version != static_cast(GLOBAL_OBJECTIVE_VERSION) || neural_schema != static_cast(GLOBAL_NEURAL_SCHEMA_VERSION) || !has_expected) { error = "Legacy replay marker is incomplete or belongs to another schema"; return false; } if (old_complete_marker) state.complete = true; return true; } bool count_published_legacy_records( const LocalRepairOptions& options, int topology, const LegacyReplayExportState& marker, std::uint64_t& count, std::string& error ) { count = 0; for (const std::string& run_id : marker.run_ids) { const std::filesystem::path directory = std::filesystem::path(options.global_dir) / "runs" / run_id / "training"; std::error_code exists_error; if (!std::filesystem::exists(directory, exists_error)) { if (exists_error) { error = "Cannot inspect legacy replay shard directory: " + exists_error.message(); return false; } continue; } std::vector shards; std::error_code iterator_error; std::filesystem::directory_iterator iterator(directory, iterator_error); const std::filesystem::directory_iterator end; while (!iterator_error && iterator != end) { std::error_code type_error; if (iterator->is_regular_file(type_error) && !type_error && iterator->path().extension() == ".sztd") { shards.push_back(iterator->path()); } else if (type_error) { error = "Cannot inspect legacy replay shard: " + type_error.message(); return false; } iterator.increment(iterator_error); } if (iterator_error) { error = "Cannot enumerate legacy replay shards: " + iterator_error.message(); return false; } std::sort(shards.begin(), shards.end()); for (const std::filesystem::path& shard : shards) { szilassi::training::StreamCallbacks callbacks; callbacks.legacy_replay = [&](const auto& record) { if (record.context.topology == static_cast(topology)) { ++count; } return true; }; std::string read_error; if (!szilassi::training::stream_read_shard( shard, callbacks, nullptr, &read_error)) { error = "Cannot validate legacy replay shard " + shard.string() + ": " + read_error; return false; } } } return true; } bool export_legacy_neural_replay( const LocalRepairOptions& options, const std::vector& active_topologies, std::vector& states, const std::string& run_id_text, TrainingArchiveSession& archive, const std::atomic& stop_requested ) { for (int topology : active_topologies) { if (stop_requested.load(std::memory_order_relaxed)) break; if (archive.writer == nullptr || !archive.writer->collection_enabled()) { poll_training_archive_notices(archive); break; } const std::filesystem::path marker_path = legacy_replay_export_marker(options, topology); LegacyReplayExportState marker; std::string marker_error; if (!read_legacy_replay_export_state( marker_path, topology, marker, marker_error)) { std::cerr << marker_error << std::endl; return false; } if (marker.complete) continue; GlobalTopologyState& state = states[topology]; if (state.neural == nullptr) continue; if (!load_global_topology_context(topology)) { std::cerr << "Cannot load topology while exporting legacy neural replay: " << topology << std::endl; return false; } const std::vector samples = state.neural->replay_samples(); if (!marker.exists) { marker.exists = true; marker.expected_records = static_cast(samples.size()); } else if (marker.expected_records != samples.size()) { std::cerr << "Pending legacy replay export no longer matches the durable " << "model snapshot for topology " << topology << std::endl; return false; } std::uint64_t already_exported = 0; if (!count_published_legacy_records( options, topology, marker, already_exported, marker_error)) { std::cerr << marker_error << std::endl; return false; } if (already_exported > marker.expected_records) { std::cerr << "Legacy replay export contains more records than its durable " << "receipt for topology " << topology << std::endl; return false; } if (already_exported < marker.expected_records && std::find(marker.run_ids.begin(), marker.run_ids.end(), run_id_text) == marker.run_ids.end()) { marker.run_ids.push_back(run_id_text); } if (!write_training_status_file_atomic( marker_path, serialize_legacy_replay_export_state(topology, marker), &marker_error)) { std::cerr << "Cannot durably prepare legacy replay export receipt: " << marker_path << ": " << marker_error << std::endl; return false; } std::uint64_t exported = 0; bool complete = true; for (std::size_t index = static_cast(already_exported); index < samples.size(); ++index) { if (stop_requested.load(std::memory_order_relaxed)) { complete = false; break; } const auto& sample = samples[index]; if (!archive.writer->collection_enabled()) { complete = false; break; } szilassi::training::LegacyReplaySample record; record.context = make_training_context( archive, options, state, 0, szilassi::training::BatchConfig{}, 0, 0); record.input = sample.input; record.improved = sample.improved; record.defect_gain = sample.defect_gain; record.plane_target = sample.plane_target; record.move_target = sample.move_target; record.scale_target = sample.scale_target; record.value_weight = sample.value_weight; record.plane_weight = sample.plane_weight; record.move_weight = sample.move_weight; record.scale_weight = sample.scale_weight; record.sequence = sample.sequence; if (!append_training_record(&archive, record)) { complete = false; if (archive.fatal_error) return false; break; } ++exported; } if (!checkpoint_training_archive(&archive)) return false; if (!complete) { poll_training_archive_notices(archive); break; } // Publishing before the marker keeps the marker's meaning simple: // every marked topology has a complete immutable legacy export. if (!seal_training_archive(&archive)) return false; marker.complete = true; if (!write_training_status_file_atomic( marker_path, serialize_legacy_replay_export_state(topology, marker), &marker_error)) { std::cerr << "Cannot durably write legacy replay export marker: " << marker_path << ": " << marker_error << std::endl; return false; } archive.legacy_records += exported; } return !archive.fatal_error; } 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(std::floor(scaled * static_cast(bin_count))), 0, bin_count - 1); } std::uint64_t global_archive_descriptor(const GlobalMetrics& metrics) { if (!std::isfinite(metrics.energy)) { return std::numeric_limits::max(); } const std::uint64_t determinant = static_cast( logarithmic_archive_bin(metrics.min_plane_determinant, -10.0, 0.0, 8)); const std::uint64_t edge = static_cast( logarithmic_archive_bin(metrics.relative_min_edge, -8.0, 0.0, 8)); const std::uint64_t turn = static_cast( logarithmic_archive_bin(metrics.min_turn_sine, -8.0, 0.0, 8)); const std::uint64_t extent = static_cast( logarithmic_archive_bin(1.0 + metrics.max_vertex_norm, 0.0, 6.0, 8)); // Geometry conditioning describes the basin; face masks describe which // local constraints still fail. Counts and energy remain quality, so a // cell still improves monotonically after ordinary Git merging. // High-defect random states have almost arbitrary masks and would fill the // archive with noise. Failure localization becomes a descriptor only in // the useful near-goal regime. const bool localize_failures = global_defects(metrics) <= 16; const std::uint64_t crossing_signature = localize_failures ? static_cast(metrics.crossing_face_mask & 0x0fffU) : 0U; const std::uint64_t intersection_signature = localize_failures ? static_cast(metrics.intersection_face_mask & 0x0fffU) : 0U; return determinant | (edge << 3U) | (turn << 6U) | (extent << 9U) | (crossing_signature << 12U) | (intersection_signature << 24U); } 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::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(global_defects(metrics)) + 0.10 * static_cast(std::max(metrics.crossings, metrics.intersections)) + 0.01 * static_cast(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(total_pulls) + 2.0) / (static_cast(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(current_round + 1) : static_cast(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(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(std::max(0, balance_gain)); if (balance_gain == 0) { reward += 0.25 * static_cast( 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(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; } } std::array choose_strategy_weights(const GlobalTopologyState& state) { // Baseline retains exactly the original 25% safety floor. Every other // strategy remains reachable, and injected MAP/CEM/neural repair keeps a // larger floor because it has produced most verified defect reductions. std::array result{{4, 1, 1, 1, 3}}; constexpr int kTotalWeight = 16; const std::uint64_t total_observations = std::accumulate( state.strategy_observations.begin(), state.strategy_observations.end(), std::uint64_t{0}); while (std::accumulate(result.begin(), result.end(), 0) < kTotalWeight) { int selected = 0; double selected_score = -std::numeric_limits::infinity(); for (int strategy = 0; strategy < 5; ++strategy) { if (result[strategy] >= 6) { continue; } const double exploration = 0.08 * std::sqrt( std::log(static_cast(total_observations) + 2.0) / (static_cast(state.strategy_observations[strategy]) + 1.0)); const double dilution = std::sqrt(static_cast(result[strategy])); const double score = (state.strategy_reward_ema[strategy] + exploration) / dilution; if (score > selected_score) { selected = strategy; selected_score = score; } } result[selected] += 1; } return result; } void update_strategy_portfolio( GlobalTopologyState& state, const std::array& rewards, const std::array& observations ) { for (int strategy = 0; strategy < 5; ++strategy) { if (observations[strategy] == 0) { continue; } // Reward is the best verified result from this cohort in the batch: // exact defect drops dominate; balance/energy-only progress is tiny. const double observed = rewards[strategy]; state.strategy_reward_ema[strategy] = 0.94 * state.strategy_reward_ema[strategy] + 0.06 * observed; state.strategy_observations[strategy] += 1; } } 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(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 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(x.data() + face_ix * 3) = direction * std::exp(log_distance); } VectorXd make_random_global_state(RNG& rng) { VectorXd x(static_cast(g_polys.size() * 3)); for (int face_ix = 0; face_ix < static_cast(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& normal ) { const int face_count = static_cast(g_polys.size()); const double s = step * (large_jump ? 3.5 : 1.0); const int move_type = std::uniform_int_distribution(0, 6)(rng); const int a = std::uniform_int_distribution(0, face_count - 1)(rng); const int b = std::uniform_int_distribution(0, face_count - 1)(rng); auto block = [&](int face_ix) { return Eigen::Map(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::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.crossing_face_mask); metrics.intersections = count_edge_face_intersections_strict( scratch.verts, *metric_planes, 1e-8, &metrics.intersection_loss, &metrics.intersection_face_mask); metrics.precise = false; if (canonical && metrics.crossings + metrics.intersections <= PRECISE_DEFECT_THRESHOLD) { metrics.crossings = count_self_crossings_precise( scratch.verts, *metric_planes, 1e-12, &metrics.crossing_face_mask); metrics.intersections = count_edge_face_intersections_precise( scratch.verts, *metric_planes, 1e-12, &metrics.intersection_loss, &metrics.intersection_face_mask); metrics.precise = true; } 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::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::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.worst_degeneracy = worst_degeneracy_barrier; metrics.condition_number = condition; metrics.condition_penalty = 0.0010 * std::min(20.0, std::log1p(condition)); metrics.degeneracy_penalty = g_global_degeneracy_weight * (worst_degeneracy_barrier + 0.05 * secondary_degeneracy_barriers); metrics.geometry_penalty = metrics.condition_penalty + 0.0002 * std::min(20.0, std::log1p(max_vertex_norm)) + metrics.degeneracy_penalty; metrics.energy = 0.05 * static_cast(global_defects(metrics)) + 0.01 * static_cast(metrics.crossings) + 0.02 * std::min(10.0, metrics.crossing_loss) + // The smooth I guide is bounded below the exact 0.05 cost of one // defect. It only directs motion inside a fixed exact-I stratum. 0.01 * std::min(4.0, metrics.intersection_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(x[component]); if (!std::isfinite(value)) { return false; } exact[component] = static_cast(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& candidates ) { if (candidates.empty()) { return; } std::vector 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(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(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(candidate.x[component]); std::uint32_t bits = 0; std::memcpy(&bits, &value, sizeof(bits)); hash ^= static_cast(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( ranked.size(), std::max(4, ranked.size() / 6)); std::array elite_mean{}; std::array elite_variance{}; double weight_sum = 0.0; for (std::size_t rank = 0; rank < elite_count; ++rank) { const double weight = std::log(static_cast(elite_count) + 1.5) - std::log(static_cast(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(elite_count) + 1.5) - std::log(static_cast(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 make_base_hybrid_seed_pool( GlobalTopologyState& state, std::uint64_t seed, int maximum_states ) { std::vector result; if (maximum_states <= 0) { return result; } result.reserve(static_cast(maximum_states)); RNG rng(static_cast(seed)); std::normal_distribution normal(0.0, 1.0); const int archive_target = std::min( maximum_states / 2, static_cast(state.archive.size())); std::vector 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(component)] = static_cast(elite.x[component]); } elite.selections += 1; result.push_back(state_value); } while (state.cem.initialized && static_cast(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(component)] = static_cast(sample[component]); } result.push_back(state_value); } return result; } using NeuralInput = szilassi::surrogate::Input; using NeuralPrediction = szilassi::surrogate::Prediction; struct HybridSeedPool { std::vector states; std::vector training_states; std::vector inputs; std::vector metrics; std::vector injected_metrics; std::vector behavior_predictions; std::vector neural_guided; std::vector plane_actions; std::vector move_actions; std::vector scale_actions; std::vector plane_propensities; std::vector move_propensities; std::vector scale_propensities; std::vector> plane_sampling_probabilities; std::vector> move_sampling_probabilities; std::vector> scale_sampling_probabilities; }; NeuralInput make_neural_input( const VectorXd& x, const GlobalMetrics& metrics, bool guided_context ) { NeuralInput input{}; for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) { const double value = component < x.size() ? x[component] : 0.0; input[static_cast(component)] = static_cast( std::clamp(std::isfinite(value) ? value : 0.0, -8.0, 8.0)); } const auto finite_or_zero = [](double value) { return std::isfinite(value) ? value : 0.0; }; input[36] = static_cast(std::clamp(metrics.crossings / 16.0, 0.0, 4.0)); input[37] = static_cast(std::clamp(metrics.intersections / 32.0, 0.0, 4.0)); input[38] = static_cast(std::log1p(std::max(0.0, finite_or_zero(metrics.crossing_loss)))); input[39] = static_cast(std::log1p(std::max(0.0, finite_or_zero(metrics.intersection_loss)))); input[40] = static_cast(std::clamp( (std::log10(std::max(1.0e-12, finite_or_zero(metrics.min_plane_determinant))) + 12.0) / 12.0, 0.0, 1.0)); input[41] = static_cast(std::clamp( finite_or_zero(metrics.relative_min_edge), 0.0, 1.0)); input[42] = static_cast(std::clamp( finite_or_zero(metrics.min_turn_sine), 0.0, 1.0)); input[43] = static_cast(std::clamp( std::log1p(std::max(0.0, finite_or_zero(metrics.max_vertex_norm))) / 12.0, 0.0, 4.0)); input[44] = guided_context ? 1.0f : 0.0f; return input; } template int sample_neural_action( const std::array& probabilities, double exploration, RNG& rng, double* selected_probability = nullptr, const std::array* extra_bias = nullptr, std::array* sampling_probabilities = nullptr ) { std::array weights{}; for (std::size_t index = 0; index < N; ++index) { const double learned = std::isfinite(probabilities[index]) ? std::max(0.0, static_cast(probabilities[index])) : 0.0; weights[index] = exploration / static_cast(N) + (1.0 - exploration) * learned; if (extra_bias != nullptr) { weights[index] += std::max(0.0, (*extra_bias)[index]); } } double total_weight = std::accumulate(weights.begin(), weights.end(), 0.0); if (total_weight <= 0.0) { weights.fill(1.0); total_weight = static_cast(N); } if (sampling_probabilities != nullptr) { for (std::size_t index = 0; index < N; ++index) { (*sampling_probabilities)[index] = static_cast(std::clamp( weights[index] / total_weight, 0.0, 1.0)); } } const int selected = std::discrete_distribution(weights.begin(), weights.end())(rng); if (selected_probability != nullptr) { *selected_probability = std::clamp( weights[static_cast(selected)] / total_weight, 0.0, 1.0); } return selected; } void apply_neural_repair_move( VectorXd& x, int selected_plane, int move_type, int scale_index, RNG& rng, std::normal_distribution& normal ) { constexpr std::array kStepScales{{0.003, 0.012, 0.05}}; const int face_count = static_cast(g_polys.size()); const int a = std::clamp(selected_plane, 0, face_count - 1); int b = std::uniform_int_distribution(0, face_count - 2)(rng); if (b >= a) { ++b; } const double step = kStepScales[static_cast( std::clamp(scale_index, 0, static_cast(kStepScales.size()) - 1))]; auto block = [&](int face) { return Eigen::Map(x.data() + face * 3); }; const auto noise = [&]() { return Vector3d(normal(rng), normal(rng), normal(rng)); }; if (move_type == 0) { block(a) += noise() * step; } else if (move_type == 1) { const Vector3d delta = noise() * step; block(a) += delta; block(b) -= delta; } else if (move_type == 2) { const Vector3d delta = noise() * step; block(a) += delta; block(b) += delta; } else if (move_type == 3) { for (int face = 0; face < face_count; ++face) { block(face) += noise() * (step * 0.18); } } else { block(a) *= std::exp(std::clamp(normal(rng) * step * 2.0, -0.25, 0.25)); } normalize_global_plane_state(x); } HybridSeedPool make_hybrid_seed_pool( GlobalTopologyState& topology_state, std::uint64_t seed, int maximum_states ) { HybridSeedPool pool; if (maximum_states <= 0) { return pool; } pool.states.reserve(static_cast(maximum_states)); pool.training_states.reserve(static_cast(maximum_states)); pool.inputs.reserve(static_cast(maximum_states)); pool.metrics.reserve(static_cast(maximum_states)); pool.injected_metrics.reserve(static_cast(maximum_states)); pool.behavior_predictions.reserve(static_cast(maximum_states)); pool.neural_guided.reserve(static_cast(maximum_states)); pool.plane_actions.reserve(static_cast(maximum_states)); pool.move_actions.reserve(static_cast(maximum_states)); pool.scale_actions.reserve(static_cast(maximum_states)); pool.plane_propensities.reserve(static_cast(maximum_states)); pool.move_propensities.reserve(static_cast(maximum_states)); pool.scale_propensities.reserve(static_cast(maximum_states)); pool.plane_sampling_probabilities.reserve(static_cast(maximum_states)); pool.move_sampling_probabilities.reserve(static_cast(maximum_states)); pool.scale_sampling_probabilities.reserve(static_cast(maximum_states)); auto append = [&](const VectorXd& injected_x, const GlobalMetrics& injected_metrics, const VectorXd& training_x, const NeuralInput& training_input, const GlobalMetrics& training_metrics, const NeuralPrediction& behavior_prediction, bool guided, int plane_action, int move_action, int scale_action, double plane_propensity, double move_propensity, double scale_propensity, const std::array& plane_sampling_probabilities, const std::array& move_sampling_probabilities, const std::array& scale_sampling_probabilities) { if (static_cast(pool.states.size()) >= maximum_states || injected_x.size() != GLOBAL_PLANE_VALUE_COUNT || training_x.size() != GLOBAL_PLANE_VALUE_COUNT || !injected_x.allFinite() || !std::isfinite(injected_metrics.energy) || !training_x.allFinite() || !std::isfinite(training_metrics.energy)) { return; } cuda_search::PlaneState value; cuda_search::PlaneState training_value; for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) { value.values[static_cast(component)] = static_cast(injected_x[component]); training_value.values[static_cast(component)] = static_cast(training_x[component]); } pool.states.push_back(value); pool.training_states.push_back(training_value); pool.inputs.push_back(training_input); pool.metrics.push_back(training_metrics); pool.injected_metrics.push_back(injected_metrics); pool.behavior_predictions.push_back(behavior_prediction); pool.neural_guided.push_back(guided ? 1u : 0u); pool.plane_actions.push_back(plane_action); pool.move_actions.push_back(move_action); pool.scale_actions.push_back(scale_action); pool.plane_propensities.push_back(static_cast(plane_propensity)); pool.move_propensities.push_back(static_cast(move_propensity)); pool.scale_propensities.push_back(static_cast(scale_propensity)); pool.plane_sampling_probabilities.push_back(plane_sampling_probabilities); pool.move_sampling_probabilities.push_back(move_sampling_probabilities); pool.scale_sampling_probabilities.push_back(scale_sampling_probabilities); }; // Base seeds do not sample an action. A zero distribution records that // fact explicitly; a uniform distribution would falsely imply behavior // policy data where no policy decision occurred. const std::array no_plane_sampling_distribution{}; const std::array no_move_sampling_distribution{}; const std::array no_scale_sampling_distribution{}; const int guided_target = topology_state.has_state ? std::min(16, std::max(1, maximum_states / 4)) : 0; const std::vector base = make_base_hybrid_seed_pool( topology_state, seed ^ 0x243f6a8885a308d3ULL, maximum_states - guided_target); PlaneEvaluationScratch scratch; for (const cuda_search::PlaneState& value : base) { VectorXd x(GLOBAL_PLANE_VALUE_COUNT); for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) { x[component] = static_cast(value.values[static_cast(component)]); } const GlobalMetrics metrics = evaluate_global_state(x, true, scratch); append( x, metrics, x, make_neural_input(x, metrics, false), metrics, NeuralPrediction{}, false, -1, -1, -1, 1.0, 1.0, 1.0, no_plane_sampling_distribution, no_move_sampling_distribution, no_scale_sampling_distribution); } if (guided_target == 0) { return pool; } struct Anchor { VectorXd x; GlobalMetrics metrics; NeuralInput input{}; NeuralPrediction prediction{}; double acquisition = 0.0; }; std::vector anchors; anchors.push_back(Anchor{ topology_state.best_x, topology_state.best, make_neural_input(topology_state.best_x, topology_state.best, true), {}, 0.0}); std::vector elites; elites.reserve(topology_state.archive.size()); for (const auto& item : topology_state.archive) { elites.push_back(&item.second); } std::stable_sort(elites.begin(), elites.end(), [](const auto* left, const auto* right) { return better_global_metrics(left->metrics, right->metrics); }); for (std::size_t index = 0; index < std::min(31, elites.size()); ++index) { anchors.push_back(Anchor{ elites[index]->x, elites[index]->metrics, make_neural_input(elites[index]->x, elites[index]->metrics, true), {}, 0.0}); } const bool model_ready = topology_state.neural != nullptr && topology_state.neural->replay_size() >= 128; for (Anchor& anchor : anchors) { if (model_ready) { anchor.prediction = topology_state.neural->predict(anchor.input); } if (!anchor.prediction.finite) { anchor.prediction.plane_probabilities.fill( 1.0f / static_cast(cuda_search::kPlaneCount)); anchor.prediction.move_probabilities.fill( 1.0f / static_cast(szilassi::surrogate::kMoveCount)); anchor.prediction.scale_probabilities = {{0.25f, 0.55f, 0.20f}}; } const double uncertainty = std::sqrt(std::max( 0.0, static_cast(anchor.prediction.uncertainty))); anchor.acquisition = model_ready ? static_cast(anchor.prediction.expected_defect_gain) + 0.75 * uncertainty + 0.10 * static_cast(anchor.prediction.improvement_probability) : 0.0; } std::stable_sort(anchors.begin(), anchors.end(), [](const Anchor& left, const Anchor& right) { if (left.acquisition != right.acquisition) { return left.acquisition > right.acquisition; } return better_global_metrics(left.metrics, right.metrics); }); RNG rng(static_cast(seed)); std::normal_distribution normal(0.0, 1.0); int attempts = 0; while (std::count(pool.neural_guided.begin(), pool.neural_guided.end(), 1u) < guided_target && attempts < guided_target * 12) { ++attempts; std::size_t anchor_index = 0; if (attempts % 4 == 0) { anchor_index = std::uniform_int_distribution( 0, anchors.size() - 1)(rng); } else { anchor_index = static_cast(attempts - 1) % std::min(8, anchors.size()); } const Anchor& anchor = anchors[anchor_index]; std::array plane_bias{}; const std::uint16_t failure_mask = static_cast( anchor.metrics.crossing_face_mask | anchor.metrics.intersection_face_mask); for (std::size_t plane = 0; plane < plane_bias.size(); ++plane) { if ((failure_mask & static_cast(1u << plane)) != 0) { plane_bias[plane] = 0.35; } } double plane_propensity = 1.0; double move_propensity = 1.0; double scale_propensity = 1.0; std::array plane_sampling_probabilities{}; std::array move_sampling_probabilities{}; std::array scale_sampling_probabilities{}; const int plane = sample_neural_action( anchor.prediction.plane_probabilities, model_ready ? 0.25 : 0.65, rng, &plane_propensity, &plane_bias, &plane_sampling_probabilities); const int move = sample_neural_action( anchor.prediction.move_probabilities, model_ready ? 0.25 : 0.65, rng, &move_propensity, static_cast*>(nullptr), &move_sampling_probabilities); const int scale = sample_neural_action( anchor.prediction.scale_probabilities, model_ready ? 0.25 : 0.50, rng, &scale_propensity, static_cast*>(nullptr), &scale_sampling_probabilities); VectorXd candidate = anchor.x; apply_neural_repair_move(candidate, plane, move, scale, rng, normal); GlobalMetrics candidate_metrics; if (!round_trip_global_state_to_fp32(candidate, candidate_metrics)) { continue; } append( candidate, candidate_metrics, anchor.x, anchor.input, anchor.metrics, anchor.prediction, true, plane, move, scale, plane_propensity, move_propensity, scale_propensity, plane_sampling_probabilities, move_sampling_probabilities, scale_sampling_probabilities); } return pool; } bool complete_hybrid_seed(const HybridSeedPool& pool, std::size_t index) { return index < pool.states.size() && index < pool.training_states.size() && index < pool.inputs.size() && index < pool.metrics.size() && index < pool.injected_metrics.size() && index < pool.behavior_predictions.size() && index < pool.neural_guided.size() && index < pool.plane_actions.size() && index < pool.move_actions.size() && index < pool.scale_actions.size() && index < pool.plane_propensities.size() && index < pool.move_propensities.size() && index < pool.scale_propensities.size() && index < pool.plane_sampling_probabilities.size() && index < pool.move_sampling_probabilities.size() && index < pool.scale_sampling_probabilities.size(); } std::uint64_t archive_hybrid_seed_pool( TrainingArchiveSession* archive, const LocalRepairOptions& options, const GlobalTopologyState& state, std::uint64_t round_seed, const szilassi::training::BatchConfig& source_batch, std::uint64_t session_seed, std::uint64_t session_batch_ordinal, std::uint64_t completed_trials_offset, std::uint64_t completed_iterations_offset, const HybridSeedPool& pool, szilassi::training::SeedProposalPurpose purpose, std::uint64_t proposal_seed, std::uint64_t assigned_chain_count, const std::vector* include_mask = nullptr ) { if (archive == nullptr || archive->writer == nullptr || archive->fatal_error || !archive->writer->collection_enabled() || pool.states.empty()) { return 0; } const std::uint32_t pool_size = static_cast( std::min( pool.states.size(), std::numeric_limits::max())); szilassi::training::BatchConfig batch = source_batch; batch.injected_pool_size = pool_size; std::uint64_t archived = 0; for (std::size_t index = 0; index < pool.states.size(); ++index) { if (include_mask != nullptr && (index >= include_mask->size() || (*include_mask)[index] == 0)) { continue; } if (!complete_hybrid_seed(pool, index) || !archive->writer->collection_enabled()) { break; } szilassi::training::SeedProposal record; record.context = make_training_context( *archive, options, state, round_seed, batch, session_seed, session_batch_ordinal, completed_trials_offset, completed_iterations_offset); record.purpose = purpose; record.anchor_state = training_plane_state(pool.training_states[index]); record.anchor_metrics = training_metrics(pool.metrics[index]); record.seed_state = training_plane_state(pool.states[index]); record.seed_metrics = training_metrics(pool.injected_metrics[index]); record.neural_input = pool.inputs[index]; record.behavior = training_behavior_prediction( pool.behavior_predictions[index]); record.plane_sampling_distribution = pool.plane_sampling_probabilities[index]; record.move_sampling_distribution = pool.move_sampling_probabilities[index]; record.scale_sampling_distribution = pool.scale_sampling_probabilities[index]; record.plane_action = pool.plane_actions[index]; record.move_action = pool.move_actions[index]; record.scale_action = pool.scale_actions[index]; record.plane_propensity = pool.plane_propensities[index]; record.move_propensity = pool.move_propensities[index]; record.scale_propensity = pool.scale_propensities[index]; record.proposal_seed = proposal_seed; record.anchor_hash = training_plane_hash(record.anchor_state); record.seed_hash = training_plane_hash(record.seed_state); record.seed_index = static_cast(index); record.pool_size = pool_size; record.assigned_chains = static_cast(std::min( assigned_chain_count / pool.states.size() + (index < assigned_chain_count % pool.states.size() ? 1ULL : 0ULL), std::numeric_limits::max())); record.flags |= szilassi::training::SeedProposalHasNeuralInput; if (pool.behavior_predictions[index].finite) { record.flags |= szilassi::training::SeedProposalPredictionSupplied; } if (pool.neural_guided[index] != 0) { record.flags |= szilassi::training::SeedProposalNeuralGuided; } if (!append_training_record(archive, record)) { if (archive->fatal_error || !archive->writer->collection_enabled()) { break; } continue; } ++archived; } return archived; } szilassi::surrogate::TrainingSample make_neural_training_sample( const NeuralInput& input, const GlobalMetrics& start_metrics, const GlobalMetrics& result_metrics, int plane_action, int move_action, int scale_action, float plane_propensity, float move_propensity, float scale_propensity, std::uint64_t sequence ) { szilassi::surrogate::TrainingSample sample; sample.input = input; const int defect_gain = global_defects(start_metrics) - global_defects(result_metrics); sample.improved = defect_gain > 0 ? 1.0f : 0.0f; sample.defect_gain = static_cast(std::clamp(defect_gain, 0, 4)); sample.value_weight = 1.0f; sample.sequence = sequence; const bool useful_policy = defect_gain > 0 || (defect_gain == 0 && better_global_metrics(result_metrics, start_metrics)); const bool valid_action = plane_action >= 0 && plane_action < cuda_search::kPlaneCount && move_action >= 0 && move_action < szilassi::surrogate::kMoveCount && scale_action >= 0 && scale_action < szilassi::surrogate::kScaleCount && std::isfinite(plane_propensity) && plane_propensity > 0.0f && std::isfinite(move_propensity) && move_propensity > 0.0f && std::isfinite(scale_propensity) && scale_propensity > 0.0f; if (!useful_policy || !valid_action) { sample.plane_weight = 0.0f; sample.move_weight = 0.0f; sample.scale_weight = 0.0f; return sample; } sample.plane_target[static_cast(plane_action)] = 1.0f; sample.move_target[static_cast(move_action)] = 1.0f; sample.scale_target[static_cast(scale_action)] = 1.0f; const auto propensity_correction = [](float probability, std::size_t choices) { return std::clamp( 1.0f / (static_cast(choices) * probability), 0.25f, 4.0f); }; const float plane_correction = propensity_correction( plane_propensity, szilassi::surrogate::kPlaneCount); const float move_correction = propensity_correction( move_propensity, szilassi::surrogate::kMoveCount); const float scale_correction = propensity_correction( scale_propensity, szilassi::surrogate::kScaleCount); // Global noise perturbs every plane and does not use selected_plane, so it // must not teach the plane head an arbitrary self-reinforcing label. sample.plane_weight = move_action == 3 ? 0.0f : (defect_gain > 0 ? 3.0f : 1.0f) * plane_correction; sample.move_weight = (defect_gain > 0 ? 2.0f : 0.5f) * move_correction; sample.scale_weight = (defect_gain > 0 ? 2.0f : 0.5f) * scale_correction; return sample; } 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(global_defects(metrics)) + 0.035 * std::min(20.0, metrics.crossing_loss) + 0.035 * std::min(20.0, metrics.intersection_loss) + metrics.geometry_penalty; } bool refine_with_spsa_adam( const VectorXd& start, std::uint64_t seed, int iterations, VectorXd& refined, std::vector& trace ) { if (start.size() != GLOBAL_PLANE_VALUE_COUNT || iterations <= 0) { return false; } RNG rng(static_cast(seed)); std::uniform_int_distribution sign(0, 1); PlaneEvaluationScratch scratch; VectorXd current = start; std::array first_moment{}; std::array second_moment{}; constexpr double kBeta1 = 0.82; constexpr double kBeta2 = 0.97; trace.clear(); trace.reserve(static_cast(iterations)); for (int iteration = 0; iteration < iterations; ++iteration) { const double perturbation = 0.022 / std::pow(static_cast(iteration + 1), 0.101); const double learning_rate = 0.012 / std::pow(static_cast(iteration + 5), 0.602); VectorXd delta(GLOBAL_PLANE_VALUE_COUNT); VectorXd plus = current; VectorXd minus = current; szilassi::training::RefinementStep trace_step; trace_step.iteration = static_cast(iteration); trace_step.perturbation = perturbation; trace_step.learning_rate = learning_rate; trace_step.center_state = training_precise_plane_state(current); for (int component = 0; component < GLOBAL_PLANE_VALUE_COUNT; ++component) { delta[component] = sign(rng) == 0 ? -1.0 : 1.0; trace_step.direction[static_cast(component)] = delta[component] < 0.0 ? -1 : 1; plus[component] += perturbation * delta[component]; minus[component] -= perturbation * delta[component]; } if (!normalize_global_plane_state(plus) || !normalize_global_plane_state(minus)) { trace_step.plus_state = training_precise_plane_state(plus); trace_step.minus_state = training_precise_plane_state(minus); trace_step.updated_state = training_precise_plane_state(current); trace.push_back(std::move(trace_step)); continue; } const GlobalMetrics plus_metrics = evaluate_global_state(plus, false, scratch); const GlobalMetrics minus_metrics = evaluate_global_state(minus, false, scratch); trace_step.plus_state = training_precise_plane_state(plus); trace_step.minus_state = training_precise_plane_state(minus); trace_step.flags |= szilassi::training::RefinementStepPlusEvaluated | szilassi::training::RefinementStepMinusEvaluated; if (std::isfinite(plus_metrics.energy)) { trace_step.plus_metrics = training_metrics(plus_metrics); } if (std::isfinite(minus_metrics.energy)) { trace_step.minus_metrics = training_metrics(minus_metrics); } const double plus_value = smooth_spsa_objective(plus_metrics); const double minus_value = smooth_spsa_objective(minus_metrics); if (!std::isfinite(plus_value) || !std::isfinite(minus_value)) { trace_step.updated_state = training_precise_plane_state(current); trace.push_back(std::move(trace_step)); 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); trace_step.updated_state = training_precise_plane_state(current); trace_step.flags |= szilassi::training::RefinementStepAccepted; trace.push_back(std::move(trace_step)); } if (!normalize_global_plane_state(current)) { 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& stop_requested ) { GlobalTrialResult result; RNG rng(static_cast(seed)); std::normal_distribution normal(0.0, 1.0); std::uniform_real_distribution uniform(0.0, 1.0); PlaneEvaluationScratch scratch; const auto evaluate = [&](const VectorXd& state, bool canonical) { ++result.evaluated_states; return evaluate_global_state(state, canonical, 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(current, false); for (int retry = 0; !std::isfinite(current_metrics.energy) && retry < 8; ++retry) { current = make_random_global_state(rng); current_metrics = evaluate(current, false); } GlobalMetrics best_search = current_metrics; VectorXd best_search_x = current; result.best = evaluate(current, true); result.best_x = current; RNG parameter_rng(static_cast(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(candidate, false); 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(candidate, true); 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(current, false); step = base_step * (0.5 + uniform(rng) * 1.5); temperature = base_temperature * (0.75 + uniform(rng)); stagnant = 0; } } GlobalMetrics final_canonical = evaluate(best_search_x, true); 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& 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(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(checkpoint.plane_coefficients[static_cast(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 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, 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(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(topology)].empty() || !load_global_topology_context(topology)) { continue; } const std::uint64_t expected_fingerprint = topology_fingerprint(); std::unordered_set seen_states; for (const CheckpointRecord* record : by_topology[static_cast(topology)]) { const GlobalCheckpoint& checkpoint = record->checkpoint; if (checkpoint.plane_coefficients.size() != static_cast(GLOBAL_PLANE_VALUE_COUNT) || (checkpoint.topology_fingerprint != 0 && checkpoint.topology_fingerprint != expected_fingerprint)) { continue; } const char* raw = reinterpret_cast( 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( checkpoint.plane_coefficients[static_cast(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 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( 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( entry.plane_coefficients[static_cast(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 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( static_cast(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(state.best_x[i]); checkpoint.plane_coefficients[static_cast(i)] = value; stored_x[i] = static_cast(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 exact_entries; std::vector persisted_descriptors; std::vector 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(elite.x[component]); entry.plane_coefficients[static_cast(component)] = value; exact_x[component] = static_cast(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(static_cast(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( static_cast(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& states ) { std::vector 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\tcrossing_loss\tintersection_loss\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.crossing_loss << "\t" << state.best.intersection_loss << "\t" << state.best.degeneracy_penalty; } else { out << "-\t-\t-\t-\t0\t0\t0\t-\t-\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, std::uint64_t training_cache_accounted_bytes, std::uint64_t neural_model_reserve_bytes, bool training_collection_enabled, const szilassi::training::RecoveryReport& training_recovery ) { 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-v5\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\t" << (options.use_cuda ? "hybrid-quality-diversity-neural-v3" : "cpu-parallel-simulated-annealing") << "\n" << "control_baseline_min_fraction\t0.25\n" << "strategy_weights\tadaptive-total:16,floors:baseline4/replica1/adaptive1/pbt1/injected3,max6\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" << "fp32_population_sample\tper-producing-strategy:64," "finite-chain-bests-before-host-selection,deterministic-circular-offset," "exact-numerator-denominator-recorded\n" << "accounted_fp32_steps\titerations-plus-initialization-fresh-stagnation-injection-pbt-and-retries\n" << "map_elites_bins\tdeterminant:8,edge:8,turn:8,extent:8,crossing-face-mask:12,intersection-face-mask:12\n" << "map_elites_max_cells_per_topology\t4096\n" << "cem\tdiagonal-weighted,covariance-adaptation:0.18,new-injected-only\n" << "intersection_loss\tsegment-depth-times-boundary-depth,weight:0.01,cap:4\n" << "repair\t" << (options.use_cuda ? "injected-quarter,offending-face-biased,neural-policy-with-exploration" : "disabled-on-cpu-backend") << "\n" << "neural_enabled\t" << (options.use_cuda ? 1 : 0) << "\n" << "neural_schema\t" << GLOBAL_NEURAL_SCHEMA_VERSION << "\n" << "neural_model_format\t" << szilassi::surrogate::kModelFormatVersion << "\n" << "neural\t" << (options.use_cuda ? "per-topology-ensemble:5,residual-mlp:45-128-128-128,value-policy,online-adamw,replay:4096,recent:25%" : "disabled") << "\n" << "neural_safety\tbaseline-floor:25%,policy-exploration:25%,exact-cuda-plus-cpu-dd-authoritative\n" << "training_archive_schema\t" << TRAINING_ARCHIVE_SCHEMA_VERSION << "\n" << "training_archive_format\t" << szilassi::training::kTrainingShardFormatVersion << "\n" << "training_archive_layout\trun-uuid/immutable-64MiB-crc-shards-plus-durable-wal\n" << "training_archive_records\tseed-proposals,stratified-fp32-chain-bests," "cpu-verified,anchor-injected-result,spsa-steps,legacy-replay\n" << "training_cache_limit_bytes\t" << TRAINING_CACHE_LIMIT_BYTES << "\n" << "training_cache_accounted_bytes_at_start\t" << training_cache_accounted_bytes << "\n" << "training_neural_model_reserve_bytes\t" << neural_model_reserve_bytes << "\n" << "training_collection_enabled_at_start\t" << (training_collection_enabled ? 1 : 0) << "\n" << "training_recovered_wals\t" << training_recovery.wals_found << "\n" << "training_wals_deferred_at_limit\t" << training_recovery.wals_deferred_at_limit << "\n" << "training_recovered_records\t" << training_recovery.records_recovered << "\n" << "training_recovered_torn_bytes_discarded\t" << training_recovery.torn_bytes_discarded << "\n" << "training_limit_behavior\tfreeze-collection-and-online-learning;search-continues;rewrite-warning\n" << "spsa\tadam:6,cpu-double,all-plus-minus-and-updated-double-states," "c-plus-i-smooth-loss,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; } std::filesystem::path neural_model_path( const LocalRepairOptions& options, int topology ) { return std::filesystem::path(options.global_dir) / neural_model_relative_path(topology); } szilassi::surrogate::Config expected_neural_config(int topology) { szilassi::surrogate::Config config; config.seed = 0x5a17a551c0ffeeULL ^ (static_cast(topology + 1) * 0x9e3779b97f4a7c15ULL) ^ static_cast(GLOBAL_OBJECTIVE_VERSION); config.replay_capacity = 4096; return config; } bool same_neural_config( const szilassi::surrogate::Config& left, const szilassi::surrogate::Config& right) { return left.seed == right.seed && left.replay_capacity == right.replay_capacity && left.learning_rate == right.learning_rate && left.weight_decay == right.weight_decay && left.beta1 == right.beta1 && left.beta2 == right.beta2 && left.adam_epsilon == right.adam_epsilon && left.gradient_clip == right.gradient_clip && left.input_clip == right.input_clip && left.gain_loss_weight == right.gain_loss_weight; } void initialize_neural_model( const LocalRepairOptions& options, GlobalTopologyState& state ) { const szilassi::surrogate::Config config = expected_neural_config(state.topology); state.neural = std::make_unique(config); const std::filesystem::path path = neural_model_path(options, state.topology); if (!std::filesystem::exists(path)) { return; } std::string error; bool cache_valid = state.neural->load(path, &error); if (cache_valid && !same_neural_config(state.neural->config(), config)) { cache_valid = false; error = "surrogate topology/objective/configuration mismatch"; } if (!cache_valid) { std::cerr << "Ignoring invalid derived neural cache for topology " << state.topology << ": " << error << std::endl; state.neural->reset(config.seed); state.neural_dirty = true; return; } // Search trials dominate the number of replay samples, so this keeps // provenance monotonic across restarts without persisting another mutable // counter in the authoritative geometry checkpoint. state.neural_sequence = std::max( state.trials, state.neural->max_replay_sequence()); } bool persist_neural_model( const LocalRepairOptions& options, GlobalTopologyState& state ) { if (state.neural == nullptr || !state.neural_dirty) { return true; } const std::filesystem::path path = neural_model_path(options, state.topology); std::error_code directory_error; std::filesystem::create_directories(path.parent_path(), directory_error); if (directory_error) { std::cerr << "Cannot create neural checkpoint directory for topology " << state.topology << ": " << directory_error.message() << std::endl; return false; } std::string error; if (!state.neural->save_atomic(path, true, &error)) { std::cerr << "Cannot durably save neural checkpoint for topology " << state.topology << ": " << error << std::endl; return false; } state.neural_dirty = false; 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( g_tris[static_cast(vertex)][static_cast(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( g_polys[static_cast(face)][static_cast(vertex)]); } } for (int edge = 0; edge < cuda_search::kEdgeCount; ++edge) { topology.edges[edge][0] = static_cast(g_edges[edge].first); topology.edges[edge][1] = static_cast(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& stop_requested, int round_seed, int effective_chain_count, cuda_search::BatchSession& session, TrainingArchiveSession* training_archive ) { 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(); HybridSeedPool initial_pool_for_archive; std::uint64_t initial_pool_seed = 0; 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(static_cast(round_seed)); config.initial_step = static_cast(options.step); config.minimum_step = static_cast( std::max(1e-7, options.step * options.min_step_ratio)); config.cooling = static_cast(options.beta); config.initial_temperature = static_cast(options.temperature); config.minimum_temperature = static_cast(options.temperature * 0.04); config.jump_chance = static_cast(options.jump_chance); config.initial_state_jitter = static_cast(options.step * 0.10); config.stagnation_iterations = std::max(256, options.stagnation); config.degeneracy_weight = static_cast(options.degeneracy_weight); initial_pool_seed = static_cast(static_cast(round_seed)) ^ 0x6a09e667f3bcc909ULL; initial_pool_for_archive = make_hybrid_seed_pool( state, initial_pool_seed, std::min(128, effective_chain_count)); std::vector& initial_states = initial_pool_for_archive.states; 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(i)] = static_cast(state.best_x[i]); } initial_states.push_back(initial); initial_pool_for_archive.training_states.push_back(initial); initial_pool_for_archive.inputs.push_back( make_neural_input(state.best_x, state.best, false)); initial_pool_for_archive.metrics.push_back(state.best); initial_pool_for_archive.injected_metrics.push_back(state.best); initial_pool_for_archive.behavior_predictions.emplace_back(); initial_pool_for_archive.neural_guided.push_back(0); initial_pool_for_archive.plane_actions.push_back(-1); initial_pool_for_archive.move_actions.push_back(-1); initial_pool_for_archive.scale_actions.push_back(-1); initial_pool_for_archive.plane_propensities.push_back(1.0f); initial_pool_for_archive.move_propensities.push_back(1.0f); initial_pool_for_archive.scale_propensities.push_back(1.0f); initial_pool_for_archive.plane_sampling_probabilities.emplace_back(); initial_pool_for_archive.move_sampling_probabilities.emplace_back(); initial_pool_for_archive.scale_sampling_probabilities.emplace_back(); } 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(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; const std::array strategy_weights = choose_strategy_weights(state); state.last_round.strategy_weights = strategy_weights; run_config.baseline_weight = strategy_weights[0]; run_config.replica_exchange_weight = strategy_weights[1]; run_config.adaptive_move_weight = strategy_weights[2]; run_config.pbt_weight = strategy_weights[3]; run_config.injected_weight = strategy_weights[4]; 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( std::max(0.002, options.step * (depth ? 0.025 : 0.05))); run_config.injected_state_jitter = static_cast( std::max(0.002, options.step * (depth ? 0.018 : 0.04))); HybridSeedPool injected_pool = make_hybrid_seed_pool( state, static_cast(static_cast(round_seed)) + static_cast(batch + 1) * 0x9e3779b97f4a7c15ULL, 64); run_config.injected_states = injected_pool.states; state.last_round.neural_seed_count += static_cast( std::count( injected_pool.neural_guided.begin(), injected_pool.neural_guided.end(), static_cast(1))); 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; const int total_strategy_weight = std::accumulate( strategy_weights.begin(), strategy_weights.end(), 0); for (int chain = 0; chain < effective_chain_count; ++chain) { const int residue = chain % total_strategy_weight; int boundary = 0; for (int strategy = 0; strategy < 5; ++strategy) { boundary += strategy_weights[static_cast(strategy)]; if (residue < boundary) { state.last_round.strategy_evaluated[static_cast(strategy)] += static_cast(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 batch_verified; const bool batch_had_reference = state.has_state; const GlobalMetrics batch_reference = state.best; const VectorXd batch_reference_x = state.best_x; const szilassi::training::BatchConfig training_batch = make_cuda_training_batch_config( options, effective_chain_count, strategy_weights, run_config, gpu); if (!session_was_initialized && batch == 0) { state.last_round.training_seed_records += archive_hybrid_seed_pool( training_archive, options, state, static_cast( static_cast(round_seed)), training_batch, gpu.session_seed, gpu.session_batch_ordinal, new_trials, evaluated_states, initial_pool_for_archive, szilassi::training::SeedProposalPurpose::SessionInitial, initial_pool_seed, static_cast(effective_chain_count)); } const std::uint64_t injected_pool_seed = static_cast(static_cast(round_seed)) + static_cast(batch + 1) * 0x9e3779b97f4a7c15ULL; std::unordered_set verified_chain_ids; verified_chain_ids.reserve(gpu.shortlist.size() * 2U + 1U); for (const cuda_search::Candidate& candidate : gpu.shortlist) { verified_chain_ids.insert(candidate.chain_id); } for (const cuda_search::Candidate& candidate : gpu.population_sample) { if (training_archive == nullptr || training_archive->writer == nullptr || !training_archive->writer->collection_enabled()) { break; } szilassi::training::Fp32Candidate record; record.context = make_training_context( *training_archive, options, state, static_cast( static_cast(round_seed)), training_batch, gpu.session_seed, gpu.session_batch_ordinal, new_trials, evaluated_states); record.state = training_plane_state(candidate.state); record.approximate_metrics = training_cuda_metrics(candidate); record.candidate_hash = training_plane_hash(record.state); record.chain_initial_rng_state = cuda_chain_initial_rng_state( gpu.session_seed, candidate.producing_chain_id); record.rollout_rng_state = candidate.rollout_rng_state; record.rollout_rng_spare_normal = candidate.rollout_rng_spare_normal; record.rollout_rng_has_spare_normal = candidate.rollout_rng_has_spare_normal; record.chain_iterations = candidate.iterations; record.rollout_start_iterations = candidate.rollout_start_iterations; record.origin_batch_ordinal = candidate.rollout_start_batch_ordinal; record.rollout_start_batch_ordinal = candidate.rollout_start_batch_ordinal; record.parent_batch_ordinal = candidate.parent_batch_ordinal; record.best_found_iteration = candidate.best_found_iteration; record.best_found_batch_ordinal = candidate.best_found_batch_ordinal; record.strategy = static_cast(candidate.strategy); record.chain_index = static_cast(candidate.chain_id); record.producing_chain_index = static_cast( candidate.producing_chain_id); record.injected_seed_index = candidate.injected_seed_index >= 0 ? static_cast(candidate.injected_seed_index) : 0xffffffffU; record.pool_size = static_cast( std::max(0, candidate.origin_pool_size)); record.replica_count = static_cast( std::max(0, candidate.origin_seed_replica_count)); record.origin_kind = static_cast(candidate.origin); record.origin_seed_index = candidate.origin_seed_index >= 0 ? static_cast(candidate.origin_seed_index) : 0xffffffffU; record.parent_chain_index = candidate.parent_chain_id != UINT64_MAX ? static_cast(candidate.parent_chain_id) : 0xffffffffU; record.sampling_numerator = candidate.population_sample_numerator; record.sampling_denominator = candidate.population_sample_denominator; record.sampling_probability = candidate.population_inclusion_probability; record.flags |= szilassi::training::Fp32UnbiasedPopulationSample; if (verified_chain_ids.find(candidate.chain_id) != verified_chain_ids.end()) { record.flags |= szilassi::training::Fp32SelectedForCpuVerification; } if (candidate.origin == cuda_search::OriginKind::InjectedPool && candidate.rollout_start_batch_ordinal == gpu.session_batch_ordinal && candidate.injected_seed_index >= 0 && static_cast(candidate.injected_seed_index) < injected_pool.neural_guided.size() && injected_pool.neural_guided[static_cast( candidate.injected_seed_index)] != 0) { record.flags |= szilassi::training::Fp32NeuralGuided; } if (append_training_record(training_archive, record)) { state.last_round.training_fp32_records += 1; } } std::array batch_strategy_rewards{}; std::array batch_strategy_observations{}; struct NeuralVerifiedRollout { bool present = false; GlobalMetrics metrics; cuda_search::PlaneState state{}; szilassi::training::ApproximateMetrics cuda_metrics{}; std::uint64_t chain_id = 0; std::uint64_t iterations = 0; std::uint64_t rollout_rng_state = 0; float rollout_rng_spare_normal = 0.0f; std::uint32_t rollout_rng_has_spare_normal = 0; std::uint64_t rollout_start_iterations = 0; std::uint64_t origin_batch_ordinal = 0; std::uint32_t pool_size = 0; std::uint32_t replica_count = 0; bool entered_archive = false; bool improved_best = false; }; std::vector neural_rollouts( injected_pool.states.size()); const auto make_cuda_verified_record = [&] ( const cuda_search::Candidate& candidate, std::size_t shortlist_rank, std::size_t strategy) { szilassi::training::VerifiedCandidate record; record.context = make_training_context( *training_archive, options, state, static_cast( static_cast(round_seed)), training_batch, gpu.session_seed, gpu.session_batch_ordinal, new_trials, evaluated_states); record.backend = szilassi::training::BackendSource::CudaShortlist; record.strategy = static_cast(strategy); record.chain_index = static_cast(candidate.chain_id); record.injected_seed_index = candidate.injected_seed_index >= 0 ? static_cast(candidate.injected_seed_index) : 0xffffffffU; record.shortlist_rank = static_cast(shortlist_rank); record.chain_seed = cuda_chain_initial_rng_state( gpu.session_seed, candidate.producing_chain_id); record.chain_iterations = candidate.iterations; record.rollout_rng_state = candidate.rollout_rng_state; record.rollout_rng_spare_normal = candidate.rollout_rng_spare_normal; record.rollout_rng_has_spare_normal = candidate.rollout_rng_has_spare_normal; record.rollout_start_iterations = candidate.rollout_start_iterations; record.origin_batch_ordinal = candidate.rollout_start_batch_ordinal; record.rollout_start_batch_ordinal = candidate.rollout_start_batch_ordinal; record.parent_batch_ordinal = candidate.parent_batch_ordinal; record.best_found_iteration = candidate.best_found_iteration; record.best_found_batch_ordinal = candidate.best_found_batch_ordinal; record.origin_kind = static_cast(candidate.origin); record.origin_seed_index = candidate.origin_seed_index >= 0 ? static_cast(candidate.origin_seed_index) : 0xffffffffU; record.origin_pool_size = static_cast( std::max(0, candidate.origin_pool_size)); record.origin_seed_replica_count = static_cast( std::max(0, candidate.origin_seed_replica_count)); record.parent_chain_index = candidate.parent_chain_id != UINT64_MAX ? static_cast(candidate.parent_chain_id) : 0xffffffffU; record.producing_chain_index = static_cast( candidate.producing_chain_id); if (batch_had_reference) { record.reference_state = training_plane_state(batch_reference_x); record.reference_metrics = training_metrics(batch_reference); record.reference_hash = training_plane_hash(record.reference_state); record.flags |= szilassi::training::CandidateHasReference; } record.candidate_state = training_plane_state(candidate.state); record.candidate_hash = training_plane_hash(record.candidate_state); record.cuda_metrics = training_cuda_metrics(candidate); record.flags |= szilassi::training::CandidateSelectedForCpuVerification; if (candidate.origin == cuda_search::OriginKind::InjectedPool && candidate.rollout_start_batch_ordinal == gpu.session_batch_ordinal && candidate.injected_seed_index >= 0 && static_cast(candidate.injected_seed_index) < injected_pool.neural_guided.size() && injected_pool.neural_guided[static_cast( candidate.injected_seed_index)] != 0) { record.flags |= szilassi::training::CandidateNeuralGuided; } return record; }; for (std::size_t shortlist_rank = 0; shortlist_rank < gpu.shortlist.size(); ++shortlist_rank) { const cuda_search::Candidate& candidate = gpu.shortlist[shortlist_rank]; const std::size_t strategy = std::min( static_cast(candidate.strategy), state.last_round.strategy_verified.size() - 1); VectorXd x(GLOBAL_PLANE_VALUE_COUNT); for (int i = 0; i < GLOBAL_PLANE_VALUE_COUNT; ++i) { x[i] = static_cast( candidate.state.values[static_cast(i)]); } GlobalMetrics verified; bool verification_failed = !valid_plane_state(x); if (!verification_failed) { verified = evaluate_global_state(x, true, scratch); verification_failed = !std::isfinite(verified.energy); } if (verification_failed) { if (training_archive != nullptr && training_archive->writer != nullptr && training_archive->writer->collection_enabled()) { szilassi::training::VerifiedCandidate record = make_cuda_verified_record( candidate, shortlist_rank, strategy); record.flags |= szilassi::training::CandidateVerificationFailed; if (append_training_record(training_archive, record)) { state.last_round.training_verified_records += 1; } } continue; } state.last_round.verified_candidates += 1; state.last_round.strategy_verified[strategy] += 1; if (strategy < 5 && candidate.best_found_batch_ordinal == gpu.session_batch_ordinal) { batch_strategy_observations[strategy] += 1; double reward = batch_had_reference ? 0.0 : 1.0; if (batch_had_reference) { const int defect_gain = global_defects(batch_reference) - global_defects(verified); if (defect_gain > 0) { reward = static_cast(std::min(defect_gain, 4)); } else if (defect_gain == 0) { const int balance_gain = std::max(batch_reference.crossings, batch_reference.intersections) - std::max(verified.crossings, verified.intersections); if (balance_gain > 0) { reward = 0.20 * static_cast(std::min(balance_gain, 2)); } else if (better_global_metrics(verified, batch_reference)) { reward = 0.02; } } } batch_strategy_rewards[strategy] = std::max( batch_strategy_rewards[strategy], reward); } const bool archive_improved = add_global_archive_elite( state, verified, x, {}, static_cast(static_cast(round_seed)), candidate.iterations, true); if (archive_improved) { state.last_round.archive_improvements += 1; state.last_round.strategy_archive_improvements[strategy] += 1; } const bool improved_global = !state.has_state || better_global_metrics(verified, state.best); batch_verified.push_back(GlobalVerifiedCandidate{ verified, x, static_cast(strategy), candidate.chain_id}); if (training_archive != nullptr && training_archive->writer != nullptr && training_archive->writer->collection_enabled()) { szilassi::training::VerifiedCandidate record = make_cuda_verified_record( candidate, shortlist_rank, strategy); record.verified_metrics = training_metrics(verified); if (improved_global) { record.flags |= szilassi::training::CandidateImprovedBest; } if (archive_improved) { record.flags |= szilassi::training::CandidateEnteredArchive; } if (append_training_record(training_archive, record)) { state.last_round.training_verified_records += 1; } } if (candidate.origin == cuda_search::OriginKind::InjectedPool && candidate.rollout_start_batch_ordinal == gpu.session_batch_ordinal && candidate.origin_pool_size == static_cast(neural_rollouts.size()) && candidate.injected_seed_index >= 0 && static_cast(candidate.injected_seed_index) < neural_rollouts.size()) { const std::size_t seed_index = static_cast( candidate.injected_seed_index); NeuralVerifiedRollout& rollout = neural_rollouts[seed_index]; if (!rollout.present || better_global_metrics(verified, rollout.metrics)) { rollout.present = true; rollout.metrics = verified; rollout.state = candidate.state; rollout.cuda_metrics = training_cuda_metrics(candidate); rollout.chain_id = candidate.producing_chain_id; rollout.iterations = candidate.best_found_iteration >= candidate.rollout_start_iterations ? candidate.best_found_iteration - candidate.rollout_start_iterations : 0; rollout.rollout_rng_state = candidate.rollout_rng_state; rollout.rollout_rng_spare_normal = candidate.rollout_rng_spare_normal; rollout.rollout_rng_has_spare_normal = candidate.rollout_rng_has_spare_normal; rollout.rollout_start_iterations = candidate.rollout_start_iterations; rollout.origin_batch_ordinal = candidate.rollout_start_batch_ordinal; rollout.pool_size = static_cast( std::max(0, candidate.origin_pool_size)); rollout.replica_count = static_cast( std::max(0, candidate.origin_seed_replica_count)); rollout.entered_archive = archive_improved; rollout.improved_best = improved_global; } } if (improved_global) { state.best = verified; state.best_x = x; state.has_state = true; state.last_round.strategy_global_improvements[strategy] += 1; } } std::vector pending_neural_samples; pending_neural_samples.reserve(neural_rollouts.size()); std::vector proposal_only_mask( neural_rollouts.size(), static_cast(1)); for (std::size_t seed_index = 0; seed_index < neural_rollouts.size(); ++seed_index) { const NeuralVerifiedRollout& rollout = neural_rollouts[seed_index]; if (!rollout.present || seed_index >= injected_pool.states.size() || seed_index >= injected_pool.training_states.size() || seed_index >= injected_pool.inputs.size() || seed_index >= injected_pool.metrics.size() || seed_index >= injected_pool.injected_metrics.size() || seed_index >= injected_pool.behavior_predictions.size() || seed_index >= injected_pool.neural_guided.size() || seed_index >= injected_pool.plane_actions.size() || seed_index >= injected_pool.move_actions.size() || seed_index >= injected_pool.scale_actions.size() || seed_index >= injected_pool.plane_propensities.size() || seed_index >= injected_pool.move_propensities.size() || seed_index >= injected_pool.scale_propensities.size() || seed_index >= injected_pool.plane_sampling_probabilities.size() || seed_index >= injected_pool.move_sampling_probabilities.size() || seed_index >= injected_pool.scale_sampling_probabilities.size()) { continue; } bool trajectory_archived = false; if (training_archive != nullptr && training_archive->writer != nullptr && training_archive->writer->collection_enabled()) { szilassi::training::InjectedTrajectory record; record.context = make_training_context( *training_archive, options, state, static_cast( static_cast(round_seed)), training_batch, gpu.session_seed, gpu.session_batch_ordinal, new_trials, evaluated_states); record.anchor_state = training_plane_state( injected_pool.training_states[seed_index]); record.anchor_metrics = training_metrics(injected_pool.metrics[seed_index]); record.injected_state = training_plane_state( injected_pool.states[seed_index]); record.injected_metrics = training_metrics( injected_pool.injected_metrics[seed_index]); record.result_state = training_plane_state(rollout.state); record.cuda_result_metrics = rollout.cuda_metrics; record.result_metrics = training_metrics(rollout.metrics); record.anchor_hash = training_plane_hash(record.anchor_state); record.injected_hash = training_plane_hash(record.injected_state); record.result_hash = training_plane_hash(record.result_state); record.rollout_seed = injected_pool_seed; record.rollout_rng_state = rollout.rollout_rng_state; record.rollout_rng_spare_normal = rollout.rollout_rng_spare_normal; record.rollout_rng_has_spare_normal = rollout.rollout_rng_has_spare_normal; record.rollout_iterations = rollout.iterations; record.start_iterations = rollout.rollout_start_iterations; record.origin_batch_ordinal = rollout.origin_batch_ordinal; record.injected_seed_index = static_cast(seed_index); record.result_chain_index = static_cast(rollout.chain_id); record.pool_size = rollout.pool_size; record.replica_count = rollout.replica_count; // Best CPU-verified descendant among all chains assigned to // this proposal in the batch. record.selection_rule = 1; record.plane_action = injected_pool.plane_actions[seed_index]; record.move_action = injected_pool.move_actions[seed_index]; record.scale_action = injected_pool.scale_actions[seed_index]; record.plane_propensity = injected_pool.plane_propensities[seed_index]; record.move_propensity = injected_pool.move_propensities[seed_index]; record.scale_propensity = injected_pool.scale_propensities[seed_index]; record.plane_sampling_distribution = injected_pool.plane_sampling_probabilities[seed_index]; record.move_sampling_distribution = injected_pool.move_sampling_probabilities[seed_index]; record.scale_sampling_distribution = injected_pool.scale_sampling_probabilities[seed_index]; record.behavior = training_behavior_prediction( injected_pool.behavior_predictions[seed_index]); record.neural_input = injected_pool.inputs[seed_index]; record.flags |= szilassi::training::TrajectoryHasNeuralInput; if (injected_pool.neural_guided[seed_index] != 0) { record.flags |= szilassi::training::TrajectoryNeuralGuided; } if (injected_pool.behavior_predictions[seed_index].finite) { record.flags |= szilassi::training::TrajectoryPredictionSupplied; } if (better_global_metrics( rollout.metrics, injected_pool.metrics[seed_index])) { record.flags |= szilassi::training::TrajectoryImproved; } if (rollout.entered_archive) { record.flags |= szilassi::training::TrajectoryEnteredArchive; } trajectory_archived = append_training_record(training_archive, record); if (trajectory_archived) { state.last_round.training_rollout_records += 1; proposal_only_mask[seed_index] = 0; } } if (trajectory_archived && state.neural != nullptr) { pending_neural_samples.push_back(make_neural_training_sample( injected_pool.inputs[seed_index], injected_pool.metrics[seed_index], rollout.metrics, injected_pool.plane_actions[seed_index], injected_pool.move_actions[seed_index], injected_pool.scale_actions[seed_index], injected_pool.plane_propensities[seed_index], injected_pool.move_propensities[seed_index], injected_pool.scale_propensities[seed_index], ++state.neural_sequence)); } } state.last_round.training_seed_records += archive_hybrid_seed_pool( training_archive, options, state, static_cast( static_cast(round_seed)), training_batch, gpu.session_seed, gpu.session_batch_ordinal, new_trials, evaluated_states, injected_pool, szilassi::training::SeedProposalPurpose::BatchInjection, injected_pool_seed, gpu.injected_chains, &proposal_only_mask); if (!checkpoint_training_archive(training_archive) || (training_archive != nullptr && training_archive->fatal_error)) { break; } const bool training_collection_active = training_archive != nullptr && training_archive->writer != nullptr && !training_archive->fatal_error && training_archive->writer->collection_enabled(); if (training_collection_active && state.neural != nullptr) { for (const auto& sample : pending_neural_samples) { if (state.neural->add_sample(sample)) { state.last_round.neural_samples_added += 1; state.neural_dirty = true; } } } update_diagonal_cem(state, batch_verified); update_strategy_portfolio( state, batch_strategy_rewards, batch_strategy_observations); if (training_collection_active && state.neural != nullptr && state.neural->replay_size() >= 128) { const szilassi::surrogate::TrainStats train_stats = state.neural->train(depth ? 2 : 1, 64); state.last_round.neural_training_steps += train_stats.completed_steps; state.neural_dirty = state.neural_dirty || train_stats.completed_steps != 0 || train_stats.rejected_updates != 0; } } if (completed_batches == 0) { return false; } if (depth && state.has_state && !state.last_round.backend_error && !stop_requested.load(std::memory_order_relaxed) && (training_archive == nullptr || !training_archive->fatal_error)) { state.last_round.spsa_attempted = true; const VectorXd refinement_start_x = state.best_x; const GlobalMetrics refinement_start_metrics = state.best; const std::uint64_t refinement_seed = static_cast(static_cast(round_seed)) ^ 0xbb67ae8584caa73bULL; VectorXd refined; GlobalMetrics refined_metrics; std::vector refinement_trace; const bool refinement_constructed = refine_with_spsa_adam( refinement_start_x, refinement_seed, 6, refined, refinement_trace); std::uint64_t refinement_evaluations = std::accumulate( refinement_trace.begin(), refinement_trace.end(), std::uint64_t{0}, [](std::uint64_t count, const szilassi::training::RefinementStep& step) { if ((step.flags & szilassi::training::RefinementStepPlusEvaluated) != 0) { ++count; } if ((step.flags & szilassi::training::RefinementStepMinusEvaluated) != 0) { ++count; } return count; }); bool refinement_verified = false; if (refinement_constructed) { ++refinement_evaluations; refinement_verified = round_trip_global_state_to_fp32(refined, refined_metrics); } evaluated_states += refinement_evaluations; if (refinement_verified) { const bool archive_improved = add_global_archive_elite( state, refined_metrics, refined, {}, static_cast(static_cast(round_seed)), 0, true); if (archive_improved) { state.last_round.archive_improvements += 1; } const bool improved_global = better_global_metrics( refined_metrics, state.best); if (training_archive != nullptr && training_archive->writer != nullptr && training_archive->writer->collection_enabled()) { szilassi::training::RefinementTrajectory record; szilassi::training::BatchConfig refinement_batch; refinement_batch.chain_count = 1; refinement_batch.iterations_per_chain = 6; refinement_batch.shortlist_size = 1; refinement_batch.baseline_chains = 1; refinement_batch.evaluated_states = refinement_evaluations; refinement_batch.proposal_scale = static_cast(options.step); refinement_batch.degeneracy_weight = static_cast(options.degeneracy_weight); refinement_batch.strategy_weights[0] = 1.0f; record.context = make_training_context( *training_archive, options, state, static_cast( static_cast(round_seed)), refinement_batch, 0, 0, new_trials, evaluated_states); record.method = szilassi::training::RefinementMethod::Spsa; record.start_state = training_plane_state(refinement_start_x); record.start_metrics = training_metrics(refinement_start_metrics); record.result_state = training_plane_state(refined); record.result_metrics = training_metrics(refined_metrics); record.start_hash = training_plane_hash(record.start_state); record.result_hash = training_plane_hash(record.result_state); record.seed = refinement_seed; record.iterations = 6; record.evaluated_states = refinement_evaluations; record.steps = refinement_trace; if (improved_global || archive_improved) { record.flags |= szilassi::training::RefinementAccepted; } if (improved_global) { record.flags |= szilassi::training::RefinementImprovedBest; } if (archive_improved) { record.flags |= szilassi::training::RefinementEnteredArchive; } if (append_training_record(training_archive, record)) { state.last_round.training_refinement_records += 1; } } if (improved_global) { state.best = refined_metrics; state.best_x = refined; state.last_round.spsa_accepted = true; } else if (archive_improved) { state.last_round.spsa_accepted = true; } } else if (training_archive != nullptr && training_archive->writer != nullptr && training_archive->writer->collection_enabled()) { szilassi::training::RefinementTrajectory record; szilassi::training::BatchConfig refinement_batch; refinement_batch.chain_count = 1; refinement_batch.iterations_per_chain = 6; refinement_batch.shortlist_size = 1; refinement_batch.baseline_chains = 1; refinement_batch.evaluated_states = refinement_evaluations; refinement_batch.proposal_scale = static_cast(options.step); refinement_batch.degeneracy_weight = static_cast(options.degeneracy_weight); refinement_batch.strategy_weights[0] = 1.0f; record.context = make_training_context( *training_archive, options, state, static_cast( static_cast(round_seed)), refinement_batch, 0, 0, new_trials, evaluated_states); record.method = szilassi::training::RefinementMethod::Spsa; record.start_state = training_plane_state(refinement_start_x); record.start_metrics = training_metrics(refinement_start_metrics); // The exact intermediate double states remain in steps. Reusing // the start as the summary result avoids inventing a valid FP32 // endpoint when final verification failed. record.result_state = record.start_state; record.result_metrics = record.start_metrics; record.start_hash = training_plane_hash(record.start_state); record.result_hash = record.start_hash; record.seed = refinement_seed; record.iterations = 6; record.evaluated_states = refinement_evaluations; record.steps = refinement_trace; record.flags |= szilassi::training::RefinementVerificationFailed; if (append_training_record(training_archive, record)) { state.last_round.training_refinement_records += 1; } } } checkpoint_training_archive(training_archive); 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( 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& stop_requested, int round_seed, int effective_cuda_chains, cuda_search::BatchSession* cuda_session, TrainingArchiveSession* training_archive ) { 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, training_archive); } 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 next_trial{0}; std::atomic completed{0}; std::atomic completed_iterations{0}; std::mutex best_mutex; std::vector workers; workers.reserve(worker_count); struct CpuTrialArchiveOutcome { bool completed = false; bool had_reference = false; int seed = 0; GlobalMetrics reference_metrics; VectorXd reference_x; GlobalTrialResult result; }; std::vector archive_outcomes( static_cast(trial_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; GlobalMetrics start_metrics; bool has_start = false; { std::lock_guard lock(best_mutex); has_start = shared_has_state; if (has_start) { start_x = shared_best_x; start_metrics = shared_best; } } const bool start_fresh = !has_start || (depth ? trial % 4 == 0 : trial % 8 != 0); const std::uint64_t seed_value = static_cast(static_cast(round_seed)) + static_cast(trial + 1) * 1000003ULL + static_cast(state.topology + 1) * 104729ULL; const int seed = static_cast(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{}; } CpuTrialArchiveOutcome& archive_outcome = archive_outcomes[static_cast(trial)]; archive_outcome.had_reference = has_start && !start_fresh; archive_outcome.seed = seed; archive_outcome.reference_metrics = start_metrics; archive_outcome.reference_x = start_x; archive_outcome.result = result; archive_outcome.completed = true; { std::lock_guard 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( result.evaluated_states, 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(); } const std::uint64_t round_trials = static_cast(completed.load(std::memory_order_relaxed)); const std::uint64_t round_iterations = completed_iterations.load(std::memory_order_relaxed); szilassi::training::BatchConfig training_batch = make_cpu_training_batch_config(options, trial_count, iterations); training_batch.evaluated_states = completed_iterations.load(std::memory_order_relaxed); for (std::size_t trial = 0; trial < archive_outcomes.size(); ++trial) { const CpuTrialArchiveOutcome& outcome = archive_outcomes[trial]; if (!outcome.completed || !std::isfinite(outcome.result.best.energy) || outcome.result.best_x.size() != GLOBAL_PLANE_VALUE_COUNT) { continue; } const bool archive_improved = add_global_archive_elite( state, outcome.result.best, outcome.result.best_x, {}, static_cast(static_cast(outcome.seed)), static_cast(std::max(0, outcome.result.iterations)), true); if (archive_improved) { state.last_round.archive_improvements += 1; state.last_round.strategy_archive_improvements[0] += 1; } if (training_archive == nullptr || training_archive->writer == nullptr || !training_archive->writer->collection_enabled()) { continue; } szilassi::training::VerifiedCandidate record; record.context = make_training_context( *training_archive, options, state, static_cast(static_cast(round_seed)), training_batch, 0, 0, round_trials, round_iterations); record.backend = szilassi::training::BackendSource::CpuTrialBest; record.strategy = 0; record.chain_index = static_cast(trial); record.shortlist_rank = static_cast(trial); record.chain_seed = static_cast( static_cast(outcome.seed)); record.chain_iterations = static_cast( std::max(0, outcome.result.iterations)); if (outcome.had_reference) { record.reference_state = training_plane_state(outcome.reference_x); record.reference_metrics = training_metrics(outcome.reference_metrics); record.reference_hash = training_plane_hash(record.reference_state); record.flags |= szilassi::training::CandidateHasReference; } record.candidate_state = training_plane_state(outcome.result.best_x); record.candidate_hash = training_plane_hash(record.candidate_state); record.verified_metrics = training_metrics(outcome.result.best); const bool improved_relevant_best = outcome.had_reference ? better_global_metrics( outcome.result.best, outcome.reference_metrics) : (!previous_has_state || better_global_metrics( outcome.result.best, previous_best)); if (improved_relevant_best) { record.flags |= szilassi::training::CandidateImprovedBest; } if (archive_improved) { record.flags |= szilassi::training::CandidateEnteredArchive; } if (append_training_record(training_archive, record)) { state.last_round.training_verified_records += 1; } } checkpoint_training_archive(training_archive); state.visits += 1; state.run_visits += 1; 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(static_cast(round_seed)), state.iterations, true)) { state.last_round.archive_improvements += 1; state.last_round.strategy_archive_improvements[0] += 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( 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(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 stop_requested{false}; std::atomic 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 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 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( 8, static_cast(std::floor( std::log2(static_cast(states[topology].visits) + 1.0)))); if (options.use_cuda) { initialize_neural_model(options, states[topology]); } } std::uint64_t existing_neural_bytes = 0; std::array existing_final_neural_snapshot_bytes{}; std::string training_setup_error; if (!scan_neural_cache_bytes( root, existing_neural_bytes, existing_final_neural_snapshot_bytes, training_setup_error)) { std::cerr << "Cannot account existing neural cache: " << training_setup_error << std::endl; stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 2; } const std::uint64_t maximum_neural_snapshot = szilassi::surrogate::maximum_snapshot_bytes( expected_neural_config(0), true); if (maximum_neural_snapshot > NEURAL_MODEL_BUDGET_PER_TOPOLOGY) { std::cerr << "Configured neural snapshot no longer fits its exact cache budget; " << "revise the program before collecting more training data." << std::endl; stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 2; } // Count every byte that is already present, then reserve the missing // capacity of each canonical final-model slot independently. Taking the // maximum of aggregate existing/reserved bytes would undercount orphan or // foreign files when some topology snapshots have not been created yet. std::uint64_t accounted_neural_files = existing_neural_bytes; for (const std::uint64_t existing_snapshot_bytes : existing_final_neural_snapshot_bytes) { if (existing_snapshot_bytes >= NEURAL_MODEL_BUDGET_PER_TOPOLOGY) { continue; } const std::uint64_t missing_slot_bytes = NEURAL_MODEL_BUDGET_PER_TOPOLOGY - existing_snapshot_bytes; if (missing_slot_bytes > std::numeric_limits::max() - accounted_neural_files) { std::cerr << "Neural cache byte-count overflow." << std::endl; stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 2; } accounted_neural_files += missing_slot_bytes; } if (maximum_neural_snapshot > std::numeric_limits::max() - accounted_neural_files) { std::cerr << "Neural cache byte-count overflow." << std::endl; stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 2; } // Atomic save temporarily keeps the previous final snapshot and one full // replacement. Existing orphan/foreign files remain separately accounted. const std::uint64_t neural_model_reserve = accounted_neural_files + maximum_neural_snapshot; szilassi::training::RecoveryReport training_recovery; const std::uint64_t training_archive_allocation = neural_model_reserve < TRAINING_CACHE_LIMIT_BYTES ? TRAINING_CACHE_LIMIT_BYTES - neural_model_reserve : 0; if (!recover_training_archive_tree( root / "runs", training_archive_allocation, training_recovery, training_setup_error)) { std::cerr << "Cannot recover interrupted training archive: " << training_setup_error << std::endl; stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 2; } if (training_recovery.shards_sealed != 0 || training_recovery.records_recovered != 0 || training_recovery.torn_bytes_discarded != 0) { std::cout << "Recovered " << training_recovery.records_recovered << " training records from " << training_recovery.shards_sealed << " interrupted run WAL(s); discarded " << training_recovery.torn_bytes_discarded << " torn tail byte(s)." << std::endl; } if (training_recovery.wals_deferred_at_limit != 0) { std::cout << "Preserved " << training_recovery.wals_deferred_at_limit << " committed training WAL(s) because the exact cache allocation " "has no room for their footer; collection will remain frozen." << std::endl; } std::uint64_t existing_training_bytes = 0; if (!scan_training_archive_tree_bytes( root / "runs", existing_training_bytes, training_setup_error)) { std::cerr << "Cannot account existing training archive: " << training_setup_error << std::endl; stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 2; } if (existing_training_bytes > std::numeric_limits::max() - neural_model_reserve) { std::cerr << "Training cache byte count overflow." << std::endl; stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 2; } const std::uint64_t initial_training_cache_bytes = existing_training_bytes + neural_model_reserve; const std::filesystem::path run_directory = root / "runs" / run_identity.run_id; std::filesystem::create_directories(run_directory); TrainingArchiveSession training_archive; training_archive.run_id = training_run_id_from_string(run_identity.run_id); training_archive.run_directory = run_directory; szilassi::training::WriterConfig training_writer_config; training_writer_config.shard_directory = run_directory / "training"; training_writer_config.global_cap_bytes = TRAINING_CACHE_LIMIT_BYTES; training_writer_config.initial_existing_bytes = initial_training_cache_bytes; training_writer_config.target_payload_bytes = 64U * 1024U * 1024U; training_archive.writer = std::make_unique( training_writer_config); poll_training_archive_notices(training_archive); if (training_archive.fatal_error) { stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 2; } const int effective_cuda_session_cache = options.use_cuda ? std::min( CUDA_SESSION_CACHE_LIMIT, static_cast(active_topologies.size())) : 0; if (!write_run_manifest( run_directory, options, run_identity, effective_cuda_chains, effective_cuda_session_cache, training_archive.writer->committed_bytes(), neural_model_reserve, training_archive.writer->collection_enabled(), training_recovery)) { 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" << "\tbest_crossing_loss\tbest_intersection_loss" << "\tneural_replay\tneural_positive\tneural_samples\tneural_train_steps\tneural_seeds" << "\ttraining_verified\ttraining_fp32\ttraining_seeds" << "\ttraining_rollouts\ttraining_refinements" << "\ttraining_cache_accounted_bytes\ttraining_wal_bytes" << "\ttraining_collection_enabled\ttraining_limit_reached" << "\tweight_baseline\tweight_replica\tweight_adaptive\tweight_pbt\tweight_injected" << "\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(); if (!export_legacy_neural_replay( options, active_topologies, states, run_identity.run_id, training_archive, stop_requested)) { std::cerr << "Legacy neural replay export failed; stopping safely." << std::endl; seal_training_archive(&training_archive); stop_requested.store(true, std::memory_order_relaxed); watcher_done.store(true, std::memory_order_relaxed); watcher.join(); return 4; } if (training_archive.legacy_records != 0) { std::cout << "Archived " << training_archive.legacy_records << " retained neural replay samples for long-term training." << std::endl; } 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; if (options.use_cuda) { std::cout << "Guidance : smooth-I repair + per-topology online neural ensemble" << std::endl; } std::cout << "Training data: CPU-verified states, sampled chain bests and injected trajectories, " << training_archive.writer->committed_bytes() << "/" << TRAINING_CACHE_LIMIT_BYTES << " cache bytes accounted" << (training_archive.writer->collection_enabled() ? "" : " (collection and online learning frozen)") << 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 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) { const bool checkpoint_changed = state.run_visits > state.last_checkpointed_run_visits; const bool archive_changed = !state.archive_dirty.empty(); const bool neural_changed = state.neural_dirty; if (state.run_visits == 0 && !archive_changed && !neural_changed) { return true; } if (!checkpoint_changed && !archive_changed && !neural_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); const bool neural_saved = persist_neural_model(options, state); if (checkpoint_saved && archive_saved && neural_saved && state.archive_dirty.empty() && !state.neural_dirty) { 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::infinity(); double current_worst_severity = -std::numeric_limits::infinity(); int current_best_defects = std::numeric_limits::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(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( 1 + (static_cast(static_cast(options.seed)) + static_cast(state.visits + 1) * 15485863ULL + static_cast(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, &training_archive); if (training_archive.fatal_error) { fatal_search_error = true; stop_requested.store(true, std::memory_order_relaxed); std::cerr << "Training archive failed; stopping search safely." << std::endl; } 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(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( 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() << (state.neural != nullptr ? ", neural replay " + std::to_string(state.neural->replay_size()) : std::string{}) << ", 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) << "\t" << state.best.crossing_loss << "\t" << state.best.intersection_loss << "\t" << (state.neural != nullptr ? state.neural->replay_size() : 0) << "\t" << (state.neural != nullptr ? state.neural->positive_replay_size() : 0) << "\t" << state.last_round.neural_samples_added << "\t" << state.last_round.neural_training_steps << "\t" << state.last_round.neural_seed_count << "\t" << state.last_round.training_verified_records << "\t" << state.last_round.training_fp32_records << "\t" << state.last_round.training_seed_records << "\t" << state.last_round.training_rollout_records << "\t" << state.last_round.training_refinement_records << "\t" << training_archive.writer->committed_bytes() << "\t" << training_archive.writer->pending_bytes() << "\t" << (training_archive.writer->collection_enabled() ? 1 : 0) << "\t" << (training_archive.writer->limit_reached() ? 1 : 0); for (int strategy = 0; strategy < 5; ++strategy) { metrics_log << "\t" << state.last_round.strategy_weights[strategy]; } 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 breadth_order = active_topologies; const int active_count = static_cast(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 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(i + 1), static_cast(missing.size())); } int cycle = 0; while (!stop_requested.load(std::memory_order_relaxed)) { std::vector 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::infinity(); double worst_severity = -std::numeric_limits::infinity(); int best_defects = std::numeric_limits::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 bandit_scores; for (int topology : active_topologies) { bandit_scores[topology] = topology_bandit_score( states[topology], total_scheduler_pulls, static_cast(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(options.seed + cycle * 7919)); std::shuffle(ranked.begin(), ranked.end(), shuffle_rng); } ranked.resize(std::min(keep, static_cast(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(i + 1), static_cast(ranked.size())); } cycle += 1; } watcher_done.store(true, std::memory_order_relaxed); watcher.join(); bool training_persistence_succeeded = seal_training_archive(&training_archive); if (!training_persistence_succeeded) { training_persistence_succeeded = seal_training_archive(&training_archive); } bool topology_persistence_succeeded = sweep_due_topologies( -1, true, szilassi::checkpoint::CheckpointReason::Stop); if (!topology_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. topology_persistence_succeeded = sweep_due_topologies( -1, true, szilassi::checkpoint::CheckpointReason::Stop); } const bool persistence_succeeded = training_persistence_succeeded && topology_persistence_succeeded; 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 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 normal(0.0, 1.0); std::uniform_real_distribution uniform(0.0, 1.0); std::uniform_int_distribution 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(std::ceil(std::min( minutes * 60.0, static_cast(std::numeric_limits::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; }