pdp_lns 0.1.1

use crate::crossover;
use crate::ges;
use crate::infeasible_tabu;
use crate::initial_construction;
use crate::instance::Instance;
use crate::lns;
use crate::local_search::{self, LsStats};
use crate::set_covering::RoutePool;
use crate::shaking;
use crate::solution::{PenaltyState, RouteInfo, Solution};
use crate::vlsn;
use rand::prelude::*;
use rand::rngs::SmallRng;
use std::collections::HashSet;
use std::thread::ScopedJoinHandle;
use std::time::{Duration, Instant};

/// Solve SC only after accumulating this many new routes.
const SC_MIN_NEW_ROUTES: usize = 48;
/// Force an SC attempt at least every N main-loop iterations.
const SC_MAX_ITERS_BETWEEN_SOLVES: usize = 10;
/// Also try SC under deeper stagnation (even if batch size is not reached).
const SC_STAGNATION_TRIGGER: u32 = 10;
/// Trigger tabu-search kick on prolonged stagnation.
const TABU_STAGNATION_TRIGGER: u32 = 7;
/// Retry interval for tabu kicks during stagnation.
const TABU_RETRY_GAP_ITERS: usize = 4;
/// Trigger LCS-SREX crossover during stagnation.
const CREX_STAGNATION_TRIGGER: u32 = 8;
/// Retry interval for crossover kicks during stagnation.
const CREX_RETRY_GAP_ITERS: usize = 8;
/// Skip CREX too close to deadline (it needs room for LS cleanup).
const CREX_MIN_REMAINING: Duration = Duration::from_secs(30);
/// Only use incumbent as 2nd parent when it is close enough to best.
const CREX_PARENT_MAX_DIST_RATIO: f64 = 1.03;
/// Perturbation fraction in regular stagnation.
const PERTURB_FRAC_LO: f64 = 0.2;
/// Perturbation fraction in deeper stagnation.
const PERTURB_FRAC_HI: f64 = 0.4;
/// Stagnation threshold to switch perturbation intensity.
const PERTURB_STAG_TRIGGER: u32 = 5;
/// Minimum number of perturbation moves.
const PERTURB_MIN_MOVES: usize = 20;
/// Tabu slice in mild stagnation.
const TABU_SLICE_MS_LO: u64 = 900;
/// Tabu slice in medium stagnation.
const TABU_SLICE_MS_MID: u64 = 1300;
/// Tabu slice in deep stagnation.
const TABU_SLICE_MS_HI: u64 = 1800;
/// Deep stagnation threshold for elite-pool SREX perturbation.
const SREX_STAGNATION_TRIGGER: u32 = 14;
/// Retry interval for elite-pool SREX perturbation.
const SREX_RETRY_GAP_ITERS: usize = 8;
/// Need at least this much time to run SREX perturbation + LS.
const SREX_MIN_REMAINING_SECS: u64 = 30;
/// Elite pool size target (requested 5-10).
const ELITE_POOL_SIZE: usize = 8;
/// Minimum diversity threshold for elite-pool insertion (broken-pairs distance).
const ELITE_MIN_BPD: f64 = 0.08;
/// SREX route-swap lower bound.
const SREX_SWAP_ROUTES_MIN: usize = 1;
/// SREX route-swap upper bound.
const SREX_SWAP_ROUTES_MAX: usize = 3;

#[derive(Clone, Copy)]
struct HyperParams {
    sc_min_new_routes: usize,
    sc_max_iters_between_solves: usize,
    sc_stagnation_trigger: u32,
    tabu_stagnation_trigger: u32,
    tabu_retry_gap_iters: usize,
    tabu_slice_ms_lo: u64,
    tabu_slice_ms_mid: u64,
    tabu_slice_ms_hi: u64,
    crex_stagnation_trigger: u32,
    crex_retry_gap_iters: usize,
    crex_min_remaining_secs: u64,
    crex_parent_max_dist_ratio: f64,
    perturb_frac_lo: f64,
    perturb_frac_hi: f64,
    perturb_stag_trigger: u32,
    perturb_min_moves: usize,
    srex_stagnation_trigger: u32,
    srex_retry_gap_iters: usize,
    srex_min_remaining_secs: u64,
    srex_swap_routes_min: usize,
    srex_swap_routes_max: usize,
    elite_pool_size: usize,
    elite_min_bpd: f64,
    op_log: bool,
}

#[derive(Clone)]
struct EliteEntry {
    sol: Solution,
    dist: f64,
}

fn env_parse_usize(name: &str, default: usize) -> usize {
    std::env::var(name)
        .ok()
        .and_then(|s| s.parse::<usize>().ok())
        .unwrap_or(default)
}

fn env_parse_u32(name: &str, default: u32) -> u32 {
    std::env::var(name)
        .ok()
        .and_then(|s| s.parse::<u32>().ok())
        .unwrap_or(default)
}

fn env_parse_u64(name: &str, default: u64) -> u64 {
    std::env::var(name)
        .ok()
        .and_then(|s| s.parse::<u64>().ok())
        .unwrap_or(default)
}

fn env_parse_f64(name: &str, default: f64) -> f64 {
    std::env::var(name)
        .ok()
        .and_then(|s| s.parse::<f64>().ok())
        .unwrap_or(default)
}

fn env_parse_bool(name: &str, default: bool) -> bool {
    match std::env::var(name) {
        Ok(v) => {
            let s = v.trim().to_ascii_lowercase();
            !(s == "0" || s == "false" || s == "no" || s == "off")
        }
        Err(_) => default,
    }
}

fn load_hyper_params() -> HyperParams {
    let sc_min_new_routes = env_parse_usize("HP_SC_MIN_NEW_ROUTES", SC_MIN_NEW_ROUTES).max(1);
    let sc_max_iters_between_solves = env_parse_usize(
        "HP_SC_MAX_ITERS_BETWEEN_SOLVES",
        SC_MAX_ITERS_BETWEEN_SOLVES,
    )
    .max(1);
    let sc_stagnation_trigger =
        env_parse_u32("HP_SC_STAGNATION_TRIGGER", SC_STAGNATION_TRIGGER).max(1);

    let tabu_stagnation_trigger =
        env_parse_u32("HP_TABU_STAGNATION_TRIGGER", TABU_STAGNATION_TRIGGER).max(1);
    let tabu_retry_gap_iters =
        env_parse_usize("HP_TABU_RETRY_GAP_ITERS", TABU_RETRY_GAP_ITERS).max(1);
    let tabu_slice_ms_lo = env_parse_u64("HP_TABU_SLICE_MS_LO", TABU_SLICE_MS_LO).max(50);
    let tabu_slice_ms_mid = env_parse_u64("HP_TABU_SLICE_MS_MID", TABU_SLICE_MS_MID).max(50);
    let tabu_slice_ms_hi = env_parse_u64("HP_TABU_SLICE_MS_HI", TABU_SLICE_MS_HI).max(50);

    let crex_stagnation_trigger =
        env_parse_u32("HP_CREX_STAGNATION_TRIGGER", CREX_STAGNATION_TRIGGER).max(1);
    let crex_retry_gap_iters =
        env_parse_usize("HP_CREX_RETRY_GAP_ITERS", CREX_RETRY_GAP_ITERS).max(1);
    let crex_min_remaining_secs =
        env_parse_u64("HP_CREX_MIN_REMAINING_SECS", CREX_MIN_REMAINING.as_secs());
    let crex_parent_max_dist_ratio =
        env_parse_f64("HP_CREX_PARENT_MAX_DIST_RATIO", CREX_PARENT_MAX_DIST_RATIO).clamp(1.0, 1.2);

    let perturb_frac_lo = env_parse_f64("HP_PERTURB_FRAC_LO", PERTURB_FRAC_LO).clamp(0.0, 1.0);
    let perturb_frac_hi = env_parse_f64("HP_PERTURB_FRAC_HI", PERTURB_FRAC_HI).clamp(0.0, 1.0);
    let perturb_stag_trigger =
        env_parse_u32("HP_PERTURB_STAG_TRIGGER", PERTURB_STAG_TRIGGER).max(1);
    let perturb_min_moves = env_parse_usize("HP_PERTURB_MIN_MOVES", PERTURB_MIN_MOVES).max(1);
    let srex_stagnation_trigger =
        env_parse_u32("HP_SREX_STAGNATION_TRIGGER", SREX_STAGNATION_TRIGGER).max(1);
    let srex_retry_gap_iters =
        env_parse_usize("HP_SREX_RETRY_GAP_ITERS", SREX_RETRY_GAP_ITERS).max(1);
    let srex_min_remaining_secs =
        env_parse_u64("HP_SREX_MIN_REMAINING_SECS", SREX_MIN_REMAINING_SECS).max(1);
    let srex_swap_routes_min =
        env_parse_usize("HP_SREX_SWAP_ROUTES_MIN", SREX_SWAP_ROUTES_MIN).max(1);
    let srex_swap_routes_max =
        env_parse_usize("HP_SREX_SWAP_ROUTES_MAX", SREX_SWAP_ROUTES_MAX).max(srex_swap_routes_min);
    let elite_pool_size = env_parse_usize("HP_ELITE_POOL_SIZE", ELITE_POOL_SIZE).clamp(5, 10);
    let elite_min_bpd = env_parse_f64("HP_ELITE_MIN_BPD", ELITE_MIN_BPD).clamp(0.0, 1.0);
    let op_log = env_parse_bool("HP_OP_LOG", true);

    HyperParams {
        sc_min_new_routes,
        sc_max_iters_between_solves,
        sc_stagnation_trigger,
        tabu_stagnation_trigger,
        tabu_retry_gap_iters,
        tabu_slice_ms_lo,
        tabu_slice_ms_mid,
        tabu_slice_ms_hi,
        crex_stagnation_trigger,
        crex_retry_gap_iters,
        crex_min_remaining_secs,
        crex_parent_max_dist_ratio,
        perturb_frac_lo,
        perturb_frac_hi,
        perturb_stag_trigger,
        perturb_min_moves,
        srex_stagnation_trigger,
        srex_retry_gap_iters,
        srex_min_remaining_secs,
        srex_swap_routes_min,
        srex_swap_routes_max,
        elite_pool_size,
        elite_min_bpd,
        op_log,
    }
}

/// Constructive heuristic used for initial solution generation.
///
/// With multi-start initialization (the default), all methods run in parallel
/// and this enum is only relevant when providing an explicit initial solution.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum InitialSolutionMethod {
    Legacy,
    LuDessouky2006,
    RopkePisinger2006,
    Cluster2004,
    Hosny2012,
}

impl InitialSolutionMethod {
    #[allow(clippy::should_implement_trait)]
    pub fn from_str(raw: &str) -> Option<Self> {
        match raw.to_ascii_lowercase().as_str() {
            "legacy" | "old" | "empty" => Some(Self::Legacy),
            "lu2006" | "lu_dessouky_2006" | "lu-dessouky-2006" | "ludessouky" => {
                Some(Self::LuDessouky2006)
            }
            "ropke2006" | "ropke_pisinger_2006" | "ropke-pisinger-2006" | "ropke" => {
                Some(Self::RopkePisinger2006)
            }
            "cluster2004" | "cluster_2004" | "cluster" | "pankratz" => Some(Self::Cluster2004),
            "hosny2012" | "initial2012" | "best2012" | "pbq2012" | "hosny" => Some(Self::Hosny2012),
            _ => None,
        }
    }

    pub fn as_str(self) -> &'static str {
        match self {
            Self::Legacy => "legacy",
            Self::LuDessouky2006 => "lu2006",
            Self::RopkePisinger2006 => "ropke2006",
            Self::Cluster2004 => "cluster2004",
            Self::Hosny2012 => "hosny2012",
        }
    }
}

#[inline(always)]
fn is_better_than_with_cached_best(
    candidate: &Solution,
    best: &Solution,
    best_dist: f64,
    inst: &Instance,
) -> bool {
    let u1 = candidate.unassigned.len();
    let u2 = best.unassigned.len();
    if u1 != u2 {
        return u1 < u2;
    }

    let v1 = candidate.num_vehicles();
    let v2 = best.num_vehicles();
    if v1 != v2 {
        return v1 < v2;
    }

    candidate.total_distance(inst) < best_dist
}

#[inline]
fn is_better_lex(candidate: &Solution, cand_dist: f64, other: &Solution, other_dist: f64) -> bool {
    let u1 = candidate.unassigned.len();
    let u2 = other.unassigned.len();
    if u1 != u2 {
        return u1 < u2;
    }
    let v1 = candidate.num_vehicles();
    let v2 = other.num_vehicles();
    if v1 != v2 {
        return v1 < v2;
    }
    cand_dist < other_dist
}

fn route_arc_set(sol: &Solution, node_base: usize) -> HashSet<u64> {
    let mut arcs = HashSet::new();
    for route in &sol.routes {
        for w in route.windows(2) {
            let a = w[0] as u64;
            let b = w[1] as u64;
            arcs.insert(a * node_base as u64 + b);
        }
    }
    arcs
}

/// Broken-pairs diversity proxy using directed route arcs.
/// 0.0 means very similar, 1.0 means highly different.
fn broken_pairs_distance(a: &Solution, b: &Solution, node_base: usize) -> f64 {
    let aa = route_arc_set(a, node_base);
    let bb = route_arc_set(b, node_base);
    if aa.is_empty() && bb.is_empty() {
        return 0.0;
    }
    let common = aa.intersection(&bb).count() as f64;
    let denom = aa.len().max(bb.len()) as f64;
    1.0 - (common / denom)
}

fn try_insert_elite(
    elite_pool: &mut Vec<EliteEntry>,
    cand: &Solution,
    cand_dist: f64,
    node_base: usize,
    hp: &HyperParams,
) {
    if elite_pool.is_empty() {
        elite_pool.push(EliteEntry {
            sol: cand.clone(),
            dist: cand_dist,
        });
        return;
    }

    let mut closest_idx = 0usize;
    let mut closest_bpd = f64::MAX;
    for (i, e) in elite_pool.iter().enumerate() {
        let bpd = broken_pairs_distance(cand, &e.sol, node_base);
        if bpd < closest_bpd {
            closest_bpd = bpd;
            closest_idx = i;
        }
    }

    if closest_bpd < hp.elite_min_bpd {
        let cur = &elite_pool[closest_idx];
        if is_better_lex(cand, cand_dist, &cur.sol, cur.dist) {
            elite_pool[closest_idx] = EliteEntry {
                sol: cand.clone(),
                dist: cand_dist,
            };
        }
        return;
    }

    if elite_pool.len() < hp.elite_pool_size {
        elite_pool.push(EliteEntry {
            sol: cand.clone(),
            dist: cand_dist,
        });
        return;
    }

    let mut worst_idx = 0usize;
    for i in 1..elite_pool.len() {
        let w = &elite_pool[worst_idx];
        let c = &elite_pool[i];
        if is_better_lex(&w.sol, w.dist, &c.sol, c.dist) {
            worst_idx = i;
        }
    }

    let worst = &elite_pool[worst_idx];
    if is_better_lex(cand, cand_dist, &worst.sol, worst.dist) {
        elite_pool[worst_idx] = EliteEntry {
            sol: cand.clone(),
            dist: cand_dist,
        };
    }
}

/// Main algorithm: LS + AGES + LNS (Fig. 1 and Section 3 of the paper).
/// `time_limit` is the total allowed computation time.
pub fn solve(inst: &Instance, time_limit: Duration, seed: u64) -> Solution {
    solve_with_initial_method(
        inst,
        time_limit,
        seed,
        None,
        InitialSolutionMethod::Legacy,
        None,
    )
}

/// Same as `solve`, but allows passing an optional initial solution.
/// The initial solution should be validated beforehand (e.g. via `Solution::from_routes`).
pub fn solve_with_initial(
    inst: &Instance,
    time_limit: Duration,
    seed: u64,
    initial: Option<Solution>,
) -> Solution {
    solve_with_initial_method(
        inst,
        time_limit,
        seed,
        initial,
        InitialSolutionMethod::Legacy,
        None,
    )
}

/// Same as `solve_with_initial`, but allows selecting the initial construction method.
///
/// If `process_start` is provided, the deadline is computed from that instant
/// so that preprocessing time (LP, arc filtering) counts toward `time_limit`.
pub fn solve_with_initial_method(
    inst: &Instance,
    time_limit: Duration,
    seed: u64,
    initial: Option<Solution>,
    _initial_method: InitialSolutionMethod,
    process_start: Option<Instant>,
) -> Solution {
    let mut rng = SmallRng::seed_from_u64(seed);
    let hp = load_hyper_params();
    let ges_stats_log = env_parse_bool("HP_GES_STATS_LOG", false);
    let short_run = time_limit <= Duration::from_secs(180);
    let eff_crex_stagnation_trigger = if short_run {
        hp.crex_stagnation_trigger.min(3)
    } else {
        hp.crex_stagnation_trigger
    };
    let eff_crex_min_remaining_secs = if short_run {
        hp.crex_min_remaining_secs.min(5)
    } else {
        hp.crex_min_remaining_secs
    };
    let eff_srex_stagnation_trigger = if short_run {
        hp.srex_stagnation_trigger.min(4)
    } else {
        hp.srex_stagnation_trigger
    };
    let eff_srex_min_remaining_secs = if short_run {
        hp.srex_min_remaining_secs.min(5)
    } else {
        hp.srex_min_remaining_secs
    };
    let eff_tabu_stagnation_trigger = if short_run {
        hp.tabu_stagnation_trigger.min(3)
    } else {
        hp.tabu_stagnation_trigger
    };

    // Use process_start (if given) so that preprocessing + init time counts
    // toward the total budget, preventing large instances from overshooting.
    let start = process_start.unwrap_or_else(Instant::now);
    let deadline = start + time_limit;

    // Step 1: Create or load initial solution(s) and apply local search.
    // When no explicit initial solution is provided, run all 4 constructive
    // heuristics in parallel (multi-start), polish each with LS, then pick the
    // best and seed pools with all of them.
    let disabled_penalty = PenaltyState::disabled();

    let all_solutions: Vec<(String, Solution)> = if let Some(initial_sol) = initial {
        eprintln!(
            "  Initial: {} vehicles, {} unassigned pairs (provided)",
            initial_sol.num_vehicles(),
            initial_sol.unassigned.len()
        );
        vec![("provided".to_string(), initial_sol)]
    } else {
        let results: Vec<(String, Solution)> = std::thread::scope(|s| {
            let h1 = s.spawn(|| {
                let mut sol = initial_construction::construct_initial_solution(inst);
                local_search::local_search(inst, &mut sol, &disabled_penalty, None);
                ("Lu&Dessouky 2006".to_string(), sol)
            });
            let h2 = s.spawn(|| {
                let mut sol = initial_construction::construct_initial_solution_ropke(inst);
                local_search::local_search(inst, &mut sol, &disabled_penalty, None);
                ("Ropke&Pisinger 2006".to_string(), sol)
            });
            let h3 = s.spawn(|| {
                let mut sol = initial_construction::construct_initial_solution_cluster(inst, seed);
                local_search::local_search(inst, &mut sol, &disabled_penalty, None);
                ("Cluster 2004".to_string(), sol)
            });
            let h4 = s.spawn(|| {
                let mut sol = initial_construction::construct_initial_solution_hosny2012(inst);
                local_search::local_search(inst, &mut sol, &disabled_penalty, None);
                ("Hosny 2012".to_string(), sol)
            });
            vec![
                h1.join().unwrap(),
                h2.join().unwrap(),
                h3.join().unwrap(),
                h4.join().unwrap(),
            ]
        });
        for (name, sol) in &results {
            eprintln!(
                "  Initial ({}): {} vehicles, dist={:.2}",
                name,
                sol.num_vehicles(),
                sol.total_distance(inst)
            );
        }
        results
    };

    // Pick the best solution (lexicographic: fewest vehicles, then shortest distance).
    let mut best_idx = 0;
    for (i, (_, sol)) in all_solutions.iter().enumerate().skip(1) {
        let (_, ref best_sol) = all_solutions[best_idx];
        let sol_nv = sol.num_vehicles();
        let best_nv = best_sol.num_vehicles();
        if sol_nv < best_nv
            || (sol_nv == best_nv && sol.total_distance(inst) < best_sol.total_distance(inst))
        {
            best_idx = i;
        }
    }
    let (ref best_name, _) = all_solutions[best_idx];
    eprintln!("  Best initial: {best_name} (index {best_idx})");

    // Destructure: sol is the working copy, best tracks the overall best.
    let mut sol = all_solutions[best_idx].1.clone();

    if hp.op_log {
        eprintln!(
            "  Operator profile: short_run={short_run} crex(stag>={eff_crex_stagnation_trigger}, rem>={eff_crex_min_remaining_secs}s) srex(stag>={eff_srex_stagnation_trigger}, rem>={eff_srex_min_remaining_secs}s) tabu(stag>={eff_tabu_stagnation_trigger})"
        );
    }

    // Initialize adaptive penalty for capacity relaxation.
    // penalty_cap scaled relative to distance/demand ratio so penalties are
    // meaningful across different instance scales.
    let max_dist = inst.dist_flat.iter().copied().fold(0.0f64, f64::max);
    let max_demand = f64::from(inst.demand.iter().map(|&d| d.abs()).max().unwrap_or(1));
    let initial_penalty_cap = (max_dist / max_demand).clamp(0.1, 1000.0);
    let initial_penalty_tw = initial_penalty_cap; // same scale (distance ≡ time)
    let mut penalty_state = PenaltyState::new(initial_penalty_cap, initial_penalty_tw);

    let mut best = sol.clone();
    let mut best_dist = best.total_distance(inst);

    // Route pool for set covering MIP (Goeke 2019, Section 3.6)
    let mut route_pool = RoutePool::new(inst);

    let mut iteration: usize = 0;
    let mut info_buf = RouteInfo::new();
    // Whole-run GES instrumentation (aggregated across outer iterations).
    let mut ges_total_stats = ges::GesRunStats::default();
    let mut ges_total_runs: u64 = 0;

    // Adaptive GES time: full budget on first attempt, reduced after consecutive failures
    let mut ges_failures_at_current_vehicles: u32 = 0;
    let mut last_ges_vehicle_target: usize = best.num_vehicles();

    // Number of parallel GES threads (override with GES_THREADS env var)
    let num_threads = std::env::var("GES_THREADS")
        .ok()
        .and_then(|s| s.parse::<usize>().ok())
        .unwrap_or(4)
        .max(1);
    // Runtime switch for CREX: set CREX_ENABLE=0 to disable crossover.
    let crex_enabled = std::env::var("CREX_ENABLE")
        .map(|v| v != "0")
        .unwrap_or(true);

    // Stagnation detection: count consecutive iterations without best improvement
    let mut no_improve_count: u32 = 0;
    // Diversification stagnation: counts iterations without main-thread improvement.
    // Unlike no_improve_count, this is NOT reset by worker-thread-only improvements,
    // ensuring Shake/CREX/SREX/Tabu fire even when worker threads keep finding better
    // solutions. Without this, diversification operators never activate in distance mode.
    let mut diversify_stag: u32 = 0;
    // Track consecutive iterations at the same vehicle count for distance-mode tuning.
    let mut stable_vehicle_iters: u32 = 0;
    let mut last_vehicle_count: usize = best.num_vehicles();
    // VLSN: only try once per best solution (no point rebuilding identical graph)
    let mut vlsn_tried_at_dist: f64 = f64::MAX;
    // Tabu: retry periodically while stagnation persists.
    let mut last_tabu_iter: usize = 0;
    // Crossover: retry periodically while stagnation persists.
    let mut last_crex_iter: usize = 0;
    let mut last_srex_iter: usize = 0;
    let mut last_sc_iter: usize = 0;
    let node_base = inst.n + 1;

    // Seed elite pool and route pool with ALL initial solutions.
    let mut elite_pool: Vec<EliteEntry> = Vec::new();
    for (_, init_sol) in &all_solutions {
        let d = init_sol.total_distance(inst);
        route_pool.add_solution(inst, init_sol);
        try_insert_elite(&mut elite_pool, init_sol, d, node_base, &hp);
    }

    // Persistent acceptance state: survives across lns() calls so that
    // LAHC/SA/RRT history accumulates and adaptive weights evolve.
    let mut accept_state = lns::AcceptanceState::new(best.cost(inst));
    let mut alns_state = lns::AlnsState::new();
    let mut arc_history = lns::ArcHistory::new(inst.n);
    arc_history.add_solution(&best);

    // Step 3: Main loop (while time remaining)
    while Instant::now() < deadline {
        iteration += 1;
        let mut ls_stats = LsStats::default();

        // Step 5: Try to reduce vehicles in best solution using GES
        // Adaptive time budget: 1/12 on first attempt at this vehicle count,
        // reduced after consecutive failures to leave more time for LS+LNS.
        let current_vehicles = best.num_vehicles();
        if current_vehicles != last_ges_vehicle_target {
            last_ges_vehicle_target = current_vehicles;
            ges_failures_at_current_vehicles = 0;
        }

        let ges_start = Instant::now();
        let remaining = deadline - ges_start;
        let large_instance = inst.m >= 200;
        // During stagnation, boost GES budget: reset the failure-based decay
        // so GES gets a fresh generous allocation.  GES runs in parallel with
        // LS+LNS, so the extra wall-time is essentially free.
        let effective_ges_failures = if no_improve_count >= 8 {
            // Deep stagnation: give GES first-attempt budget every time
            0
        } else if no_improve_count >= 4 {
            // Moderate stagnation: cap decay at 1 failure level
            ges_failures_at_current_vehicles.min(1)
        } else {
            ges_failures_at_current_vehicles
        };
        let ges_fraction = if large_instance {
            match effective_ges_failures {
                0 => 6,
                1 => 18,
                2 => 54,
                _ => 135,
            }
        } else {
            match effective_ges_failures {
                0 => 3,
                1 => 10,
                2 => 30,
                _ => 75,
            }
        };
        let ges_budget = remaining / ges_fraction;
        // In concurrent mode, GES runs in parallel with main-thread LS+LNS (~2-3s).
        // A minimum GES budget of 2s is essentially free since LS+LNS takes at
        // least that long. This gives GES more search time on hard reductions.
        let ges_budget = if large_instance {
            // Reduce GES floor during stagnation: faster iterations = more
            // perturbation cycles = more diverse starting points for GES.
            let floor = if ges_failures_at_current_vehicles >= 3 {
                5
            } else {
                8
            };
            ges_budget.max(Duration::from_secs(floor)).min(remaining)
        } else {
            ges_budget.max(Duration::from_secs(3)).min(remaining)
        };
        let ges_deadline = (ges_start + ges_budget).min(deadline);

        // When GES has failed many times, those threads waste compute. Switch them
        // to independent LNS searches for productive parallel distance exploration.
        // Switch to LNS threads early: after just 2 GES failures, vehicle reduction
        // is very unlikely. Spending extra time on LNS threads gives more parallel
        // distance search (~40s extra). Thread 0 still runs GES as a safety net.
        let use_lns_threads = ges_failures_at_current_vehicles >= 2;
        let use_starvation = ges::ges_use_starvation();

        // Compute LP dual values for SC-guided GES route selection
        let sc_duals: Option<Vec<f64>> = if ges::ges_use_sc_duals() && !route_pool.is_empty() {
            route_pool.solve_lp_duals(inst)
        } else {
            None
        };

        // Use more threads for LNS mode (parallel distance exploration benefits
        // from breadth), fewer for GES (diminishing returns on vehicle reduction).
        let effective_threads = if use_lns_threads {
            num_threads.max(8)
        } else {
            num_threads
        };

        // Generate independent seeds from main rng
        let seeds: Vec<u64> = (0..effective_threads).map(|_| rng.random()).collect();

        // Prepare speculative solution for the "GES fails" path.
        // On improvement or vehicle mismatch, restart from best.
        // During stagnation at the same vehicle count, keep current sol
        // so that LAHC acceptance inside LNS accumulates exploration
        // instead of resetting to the same local optimum every iteration.
        if no_improve_count == 0 || sol.num_vehicles() > best.num_vehicles() {
            sol.clone_from(&best);
            accept_state.reset();
            alns_state.reset();
        } else if hp.op_log {
            eprintln!(
                "  [iter {}] Continuing from current sol (veh={}, dist={:.2}, best_dist={:.2})",
                iteration,
                sol.num_vehicles(),
                sol.total_distance(inst),
                best_dist
            );
        }
        // Adaptive perturbation: heavier early in the search for diversification,
        // lighter near the end for intensification (survey §5.1.1).
        let elapsed_frac = start.elapsed().as_secs_f64() / time_limit.as_secs_f64();
        let time_decay = 1.0 - 0.5 * elapsed_frac; // 1.0 at start, 0.5 at end
        let base_frac = if no_improve_count >= hp.perturb_stag_trigger {
            hp.perturb_frac_hi
        } else {
            hp.perturb_frac_lo
        };
        let perturb_fraction = base_frac * time_decay;
        #[allow(
            clippy::cast_possible_truncation,
            clippy::cast_sign_loss,
            clippy::cast_precision_loss
        )]
        let num_perturb_moves =
            ((inst.m as f64 * perturb_fraction).ceil() as usize).max(hp.perturb_min_moves);
        ges::perturb_solution(inst, &mut sol, num_perturb_moves, &mut info_buf, &mut rng);

        // Run GES threads and main-thread LS+LNS concurrently.
        // GES threads work on clones of best; main thread works on perturbed sol.
        // On a multi-core machine (32 cores), this overlaps GES wall time with
        // LS+LNS work, saving 0.5-2s per iteration during vehicle reduction.
        let mut best_improved_this_iter = false;
        let mut main_improved_this_iter = false;

        let (best_ges, ges_stats_agg, ges_runs): (Option<Solution>, ges::GesRunStats, usize) =
            std::thread::scope(|s| {
                struct ThreadOutcome {
                    candidate: Option<Solution>,
                    ges_stats: Option<ges::GesRunStats>,
                }

                // Spawn worker threads: GES (vehicle reduction) or LNS (distance exploration).
                // IMPORTANT: Thread #0 ALWAYS runs GES regardless of use_lns_threads.
                // Do not disable this — vehicle reduction opportunities can appear at
                // any point during distance optimization and are worth the 1-thread cost.
                // Removing the dedicated GES thread means we can never escape a sub-optimal
                // vehicle count once we enter LNS mode.
                let handles: Vec<_> = seeds
                    .iter()
                    .enumerate()
                    .map(|(thread_idx, &seed)| {
                        // Thread 0: always GES. Threads 1..N: GES or LNS based on mode.
                        let lns_mode = use_lns_threads && thread_idx > 0;

                        // Diversify starting point:
                        // - LNS threads: half start from elite pool
                        // - GES thread: after 2+ failures, try elite pool solutions
                        //   (different route structures may enable different removals)
                        let use_elite = if lns_mode {
                            thread_idx >= effective_threads / 2
                        } else {
                            ges_failures_at_current_vehicles >= 2
                        };
                        let mut candidate = if use_elite && elite_pool.len() > 1 {
                            let best_veh = best.num_vehicles();
                            let elite_idx = (seed as usize) % elite_pool.len();
                            let entry = &elite_pool[elite_idx];
                            if entry.sol.num_vehicles() == best_veh
                                && (entry.dist - best_dist).abs() > 1e-6
                            {
                                entry.sol.clone()
                            } else {
                                best.clone()
                            }
                        } else {
                            best.clone()
                        };
                        let thread_budget = ges_budget;
                        let thread_deadline = ges_deadline;
                        let ges_failures_copy = ges_failures_at_current_vehicles;
                        let sc_duals_ref = sc_duals.as_deref();
                        s.spawn(move || {
                            let mut thread_rng = SmallRng::seed_from_u64(seed);
                            if lns_mode {
                                // Restart loop: LNS stops after MIN_NO_IMPROVE=500
                                // (~0.5s), wasting the remaining ~7s budget. Restart
                                // with fresh perturbation to get multiple independent
                                // searches per thread (like the main-thread loop).
                                let mut info_buf = RouteInfo::new();
                                let disabled = PenaltyState::disabled();
                                #[allow(
                                    clippy::cast_possible_truncation,
                                    clippy::cast_sign_loss,
                                    clippy::cast_precision_loss
                                )]
                                let n_moves = ((inst.m as f64 * 0.25).ceil() as usize).max(15);
                                ges::perturb_solution(
                                    inst,
                                    &mut candidate,
                                    n_moves,
                                    &mut info_buf,
                                    &mut thread_rng,
                                );
                                local_search::local_search(inst, &mut candidate, &disabled, None);
                                let mut thread_best = candidate.clone();
                                lns::lns(
                                    inst,
                                    &mut candidate,
                                    &mut thread_best,
                                    thread_budget,
                                    &mut thread_rng,
                                    thread_deadline,
                                    None,
                                    None,
                                    None,
                                    None,
                                    None,
                                );
                                // Restart: perturb from thread_best and search again
                                while Instant::now() < thread_deadline {
                                    candidate.clone_from(&thread_best);
                                    ges::perturb_solution(
                                        inst,
                                        &mut candidate,
                                        n_moves,
                                        &mut info_buf,
                                        &mut thread_rng,
                                    );
                                    local_search::local_search(
                                        inst,
                                        &mut candidate,
                                        &disabled,
                                        None,
                                    );
                                    lns::lns(
                                        inst,
                                        &mut candidate,
                                        &mut thread_best,
                                        thread_budget,
                                        &mut thread_rng,
                                        thread_deadline,
                                        None,
                                        None,
                                        None,
                                        None,
                                        None,
                                    );
                                }
                                if thread_best.is_feasible(inst) {
                                    ThreadOutcome {
                                        candidate: Some(thread_best),
                                        ges_stats: None,
                                    }
                                } else {
                                    ThreadOutcome {
                                        candidate: None,
                                        ges_stats: None,
                                    }
                                }
                            } else {
                                // Pre-GES starvation: evacuate easy pairs from victim route
                                if use_starvation
                                    && ges_failures_copy >= 2
                                    && let Some(starve_budget) =
                                        thread_deadline.checked_duration_since(Instant::now())
                                {
                                    let starve_deadline = Instant::now() + starve_budget / 10;
                                    let evacuated = ges::starve_route(
                                        inst,
                                        &mut candidate,
                                        &mut thread_rng,
                                        starve_deadline,
                                        sc_duals_ref,
                                    );
                                    if evacuated > 0 {
                                        local_search::local_search(
                                            inst,
                                            &mut candidate,
                                            &PenaltyState::disabled(),
                                            None,
                                        );
                                    }
                                }
                                // Normal GES: try to reduce vehicles
                                let mut ges_stats = ges::GesRunStats::default();
                                let ok = ges::guided_ejection_search(
                                    inst,
                                    &mut candidate,
                                    &mut thread_rng,
                                    thread_deadline,
                                    Some(&mut ges_stats),
                                    sc_duals_ref,
                                );
                                if ok {
                                    local_search::local_search(
                                        inst,
                                        &mut candidate,
                                        &PenaltyState::disabled(),
                                        None,
                                    );
                                    ThreadOutcome {
                                        candidate: Some(candidate),
                                        ges_stats: Some(ges_stats),
                                    }
                                } else {
                                    ThreadOutcome {
                                        candidate: None,
                                        ges_stats: Some(ges_stats),
                                    }
                                }
                            }
                        })
                    })
                    .collect();

                // Main thread: run LS + LNS in a loop while GES threads work.
                // When LNS stagnates, perturb and restart to keep doing useful work
                // instead of idling while waiting for GES threads to finish.
                let main_deadline = ges_deadline.min(deadline);
                let main_start = Instant::now();
                let mut lns_iters_total: usize = 0;
                let mut lns_rounds: usize = 0;
                let mut perturb_info_buf = RouteInfo::new();
                let dist_before_main = sol.total_distance(inst);

                // First round: full LS
                local_search::local_search_full(
                    inst,
                    &mut sol,
                    &disabled_penalty,
                    if hp.op_log { Some(&mut ls_stats) } else { None },
                );
                if sol.is_feasible(inst)
                    && is_better_than_with_cached_best(&sol, &best, best_dist, inst)
                {
                    best.clone_from(&sol);
                    best_improved_this_iter = true;
                    main_improved_this_iter = true;
                }

                if !best.is_feasible(inst) {
                    penalty_state.penalty_tw *= 2.0;
                }

                loop {
                    // Exit when all GES threads have finished — no point in idling,
                    // but also no point in blocking the next outer iteration.
                    let ges_threads_done = handles.iter().all(ScopedJoinHandle::is_finished);
                    if ges_threads_done && lns_rounds > 0 {
                        break;
                    }
                    if Instant::now() >= main_deadline {
                        break;
                    }

                    // LNS round with short budget — stagnation exit will trigger naturally
                    let round_budget = Duration::from_secs(3);
                    accept_state.select_criterion(&mut rng);
                    let iters = lns::lns(
                        inst,
                        &mut sol,
                        &mut best,
                        round_budget,
                        &mut rng,
                        main_deadline,
                        Some(&mut accept_state),
                        Some(&mut alns_state),
                        Some(&arc_history),
                        Some(&mut penalty_state),
                        Some(&mut route_pool),
                    );
                    lns_iters_total += iters;
                    lns_rounds += 1;
                    best_dist = best.total_distance(inst);

                    if sol.is_feasible(inst)
                        && is_better_than_with_cached_best(&sol, &best, best_dist, inst)
                    {
                        best.clone_from(&sol);
                        best_dist = best.total_distance(inst);
                        best_improved_this_iter = true;
                        main_improved_this_iter = true;
                    }

                    if handles.iter().all(ScopedJoinHandle::is_finished) {
                        break;
                    }
                    if Instant::now() >= main_deadline {
                        break;
                    }

                    // Perturb + quick LS to give next LNS round a fresh starting point
                    #[allow(
                        clippy::cast_possible_truncation,
                        clippy::cast_sign_loss,
                        clippy::cast_precision_loss
                    )]
                    let n_moves = ((inst.m as f64 * 0.15).ceil() as usize).max(10);
                    let snap = local_search::snapshot_routes(&sol.routes);
                    ges::perturb_solution(inst, &mut sol, n_moves, &mut perturb_info_buf, &mut rng);
                    let dirty = local_search::dirty_route_indices(&sol.routes, &snap);
                    local_search::local_search_dirty(
                        inst,
                        &mut sol,
                        &disabled_penalty,
                        if hp.op_log { Some(&mut ls_stats) } else { None },
                        &dirty,
                    );
                    if sol.is_feasible(inst)
                        && is_better_than_with_cached_best(&sol, &best, best_dist, inst)
                    {
                        best.clone_from(&sol);
                        best_improved_this_iter = true;
                        main_improved_this_iter = true;
                    }
                }

                if !best.is_feasible(inst) {
                    penalty_state.penalty_tw /= 2.0;
                }
                best_dist = best.total_distance(inst);

                // HGS-style repair: if LNS left sol infeasible, try pushing it
                // toward feasibility by re-running LS with 10x penalty (50% prob).
                if !sol.is_feasible(inst) && rng.random_bool(0.5) {
                    let repair_penalty = PenaltyState::new(
                        penalty_state.penalty_cap.max(1.0) * 10.0,
                        penalty_state.penalty_tw.max(1.0) * 10.0,
                    );
                    local_search::local_search(
                        inst,
                        &mut sol,
                        &repair_penalty,
                        if hp.op_log { Some(&mut ls_stats) } else { None },
                    );
                    if sol.is_feasible(inst)
                        && is_better_than_with_cached_best(&sol, &best, best_dist, inst)
                    {
                        best.clone_from(&sol);
                        best_dist = best.total_distance(inst);
                        best_improved_this_iter = true;
                        main_improved_this_iter = true;
                    }
                }

                let main_elapsed = main_start.elapsed();
                if hp.op_log {
                    eprintln!(
                        "  [iter {}] main-thread: {:.1}s, {} rounds, {} iters (dist {:.2}->{:.2}), best={:.2} (t={:.1}s)",
                        iteration,
                        main_elapsed.as_secs_f64(),
                        lns_rounds,
                        lns_iters_total,
                        dist_before_main,
                        sol.total_distance(inst),
                        best_dist,
                        start.elapsed().as_secs_f64(),
                    );
                }

                // Join worker threads: collect GES instrumentation + best candidate.
                let mut best_candidate: Option<(usize, f64, Solution)> = None;
                let mut stats_agg = ges::GesRunStats::default();
                let mut runs = 0usize;
                for h in handles {
                    let out = h.join().unwrap();
                    if let Some(st) = out.ges_stats {
                        stats_agg.merge_from(&st);
                        runs += 1;
                    }
                    if let Some(candidate) = out.candidate {
                        let nv = candidate.num_vehicles();
                        let dist = candidate.total_distance(inst);
                        match &best_candidate {
                            None => best_candidate = Some((nv, dist, candidate)),
                            Some((bnv, bdist, _)) => {
                                let better = nv < *bnv
                                    || (nv == *bnv
                                        && dist
                                            .partial_cmp(bdist)
                                            .is_some_and(std::cmp::Ordering::is_lt));
                                if better {
                                    best_candidate = Some((nv, dist, candidate));
                                }
                            }
                        }
                    }
                }
                let best_candidate = best_candidate.map(|(_, _, candidate)| candidate);
                (best_candidate, stats_agg, runs)
            });
        ges_total_runs += ges_runs as u64;
        ges_total_stats.merge_from(&ges_stats_agg);

        if ges_stats_log {
            eprintln!(
                "  [iter {}] GES stats: runs={}, inner_iters={}, route_attempts={}, direct_ok/fail={}/{}, checks={}, cache_hits={}, prefilter(tp/d/td)={}/{}/{}, conflict_rej={}, eject_ok/calls={}/{}, squeeze_ok/attempts={}/{}, perturbs={}, pick_easy/hard={}/{}",
                iteration,
                ges_runs,
                ges_stats_agg.inner_iterations,
                ges_stats_agg.route_attempts,
                ges_stats_agg.direct_successes,
                ges_stats_agg.direct_failures,
                ges_stats_agg.direct_route_checks,
                ges_stats_agg.direct_cache_hits,
                ges_stats_agg.direct_arc_tp_rejects,
                ges_stats_agg.direct_arc_d_rejects,
                ges_stats_agg.direct_arc_td_rejects,
                ges_stats_agg.direct_conflict_rejects,
                ges_stats_agg.ejection_successes,
                ges_stats_agg.ejection_calls,
                ges_stats_agg.squeeze_successes,
                ges_stats_agg.squeeze_attempts,
                ges_stats_agg.perturb_calls,
                ges_stats_agg.pickup_easiest_picks,
                ges_stats_agg.pickup_hardest_picks,
            );
        }
        // Check thread results
        if let Some(thread_sol) = best_ges {
            if !thread_sol.is_feasible(inst) {
                if !use_lns_threads {
                    eprintln!(
                        "  [iter {}] GES returned infeasible solution ({} vehicles), discarding",
                        iteration,
                        thread_sol.num_vehicles(),
                    );
                }
                ges_failures_at_current_vehicles += 1;
            } else if thread_sol.num_vehicles() < best.num_vehicles() {
                // Vehicle reduction — from the dedicated GES thread (#0) or
                // from an all-GES round. Accept unconditionally and reset state.
                // Feed pre-reduction routes into SC pool: these efficient routes
                // from higher vehicle counts can be recombined by SC to form
                // better solutions at lower vehicle counts.
                if best.is_feasible(inst) {
                    route_pool.add_solution(inst, &best);
                }
                best = thread_sol;
                best_dist = best.total_distance(inst);
                ges_failures_at_current_vehicles = 0;
                no_improve_count = 0;
                best_improved_this_iter = true;
                main_improved_this_iter = true;
                eprintln!(
                    "  [iter {}] GES reduced to {} vehicles, dist={:.2} (t={:.1}s, {}threads)",
                    iteration,
                    best.num_vehicles(),
                    best_dist,
                    start.elapsed().as_secs_f64(),
                    num_threads,
                );
                sol.clone_from(&best);
                accept_state.reset();
                alns_state.reset();
                try_insert_elite(&mut elite_pool, &best, best_dist, node_base, &hp);
                // Basic LS on post-GES result — quick route structure cleanup
                // without expensive intra-route operators. Saves time for the
                // next GES attempt which is higher priority.
                if Instant::now() < deadline {
                    local_search::local_search(
                        inst,
                        &mut sol,
                        &disabled_penalty,
                        if hp.op_log { Some(&mut ls_stats) } else { None },
                    );
                    if sol.is_feasible(inst)
                        && is_better_than_with_cached_best(&sol, &best, best_dist, inst)
                    {
                        best.clone_from(&sol);
                        best_dist = best.total_distance(inst);
                    }
                }
            } else if use_lns_threads {
                // LNS threads: accept if better than current best (distance improvement)
                let thread_dist = thread_sol.total_distance(inst);
                // Always feed LNS-thread results into SC pool for route diversity.
                // Thread results contain diverse route structures from independent
                // searches that the main thread wouldn't discover.
                if thread_sol.is_feasible(inst) {
                    route_pool.add_solution(inst, &thread_sol);
                }
                if is_better_than_with_cached_best(&thread_sol, &best, best_dist, inst) {
                    best = thread_sol;
                    best_dist = best.total_distance(inst);
                    best_improved_this_iter = true;
                    try_insert_elite(&mut elite_pool, &best, best_dist, node_base, &hp);
                    eprintln!(
                        "  [iter {}] LNS-thread improved: {} vehicles, dist={:.2} (t={:.1}s, {}threads)",
                        iteration,
                        best.num_vehicles(),
                        best_dist,
                        start.elapsed().as_secs_f64(),
                        effective_threads,
                    );
                } else {
                    // Even non-improving results enrich elite pool diversity
                    try_insert_elite(&mut elite_pool, &thread_sol, thread_dist, node_base, &hp);
                    if hp.op_log {
                        eprintln!(
                            "  [iter {}] LNS-threads no improvement (best thread: veh={}, dist={:.2})",
                            iteration,
                            thread_sol.num_vehicles(),
                            thread_dist,
                        );
                    }
                }
            } else {
                // All-GES mode, but no vehicle reduction found
                // (GES threads return None on failure, so this path is rare —
                //  it means a GES thread returned a same-vehicle-count solution)
                ges_failures_at_current_vehicles += 1;
            }
        } else {
            ges_failures_at_current_vehicles += 1;
        }

        if Instant::now() >= deadline {
            break;
        }

        // Distance mode: vehicle count has been stable for 3+ iterations.
        // Used to tune CREX/SREX/SC parameters for distance optimization.
        let in_distance_mode = stable_vehicle_iters >= 3;

        // Shaking operators: graduated response to stagnation
        // Level 1: inter-route 2-pair relocate with threshold acceptance.
        // Worsening threshold grows with stagnation: at stag=3 accept 0.5% worsening,
        // scaling up to 2.5% at deep stagnation. Allows structural perturbation
        // that LS can then recover from.
        if diversify_stag >= 3 && Instant::now() < deadline {
            #[allow(clippy::cast_precision_loss)]
            let theta = best_dist * 5e-3 * (f64::from(diversify_stag) / 3.0).min(5.0);
            if hp.op_log {
                eprintln!(
                    "  [iter {}] TRY Shake2pair (dstag={}, theta={:.1}, t={:.1}s)",
                    iteration,
                    diversify_stag,
                    theta,
                    start.elapsed().as_secs_f64()
                );
            }
            let shake_snap = local_search::snapshot_routes(&sol.routes);
            if shaking::shake_inter_route_2pair(inst, &mut sol, theta) {
                let shake_dirty = local_search::dirty_route_indices(&sol.routes, &shake_snap);
                local_search::local_search_full_dirty(
                    inst,
                    &mut sol,
                    &disabled_penalty,
                    if hp.op_log { Some(&mut ls_stats) } else { None },
                    &shake_dirty,
                );
                if sol.is_feasible(inst)
                    && is_better_than_with_cached_best(&sol, &best, best_dist, inst)
                {
                    best.clone_from(&sol);
                    best_dist = best.total_distance(inst);
                    best_improved_this_iter = true;
                    main_improved_this_iter = true;
                    no_improve_count = 0;
                    eprintln!(
                        "  [iter {}] Shake2pair improved: {} vehicles, dist={:.2} (t={:.1}s)",
                        iteration,
                        best.num_vehicles(),
                        best_dist,
                        start.elapsed().as_secs_f64()
                    );
                } else if hp.op_log {
                    eprintln!(
                        "  [iter {}] Shake2pair applied, no best improvement (veh={}, dist={:.2})",
                        iteration,
                        sol.num_vehicles(),
                        sol.total_distance(inst),
                    );
                }
            } else if hp.op_log {
                eprintln!("  [iter {iteration}] Shake2pair no feasible move");
            }
        }

        // LCS-SREX crossover kick: combine best with current incumbent under stagnation.
        // In distance mode, activate earlier (stag>=4) with smaller retry gap (4 iters).
        let dm_crex_trigger = if in_distance_mode {
            eff_crex_stagnation_trigger.min(4)
        } else {
            eff_crex_stagnation_trigger
        };
        let dm_crex_gap = if in_distance_mode {
            hp.crex_retry_gap_iters.min(4)
        } else {
            hp.crex_retry_gap_iters
        };
        let remaining_before_crex = deadline.saturating_duration_since(Instant::now());
        if diversify_stag >= dm_crex_trigger
            && iteration.saturating_sub(last_crex_iter) >= dm_crex_gap
            && remaining_before_crex >= Duration::from_secs(eff_crex_min_remaining_secs)
            && crex_enabled
        {
            if hp.op_log {
                eprintln!(
                    "  [iter {}] TRY CREX (dstag={}, rem={:.1}s)",
                    iteration,
                    diversify_stag,
                    remaining_before_crex.as_secs_f64()
                );
            }
            // Build parent_b: use sol if close enough, otherwise perturb best.
            let sol_vehicles = sol.num_vehicles();
            let sol_dist = sol.total_distance(inst);
            let sol_ok = sol_vehicles <= best.num_vehicles() + 1
                && sol_dist <= best_dist * hp.crex_parent_max_dist_ratio;

            let parent_b = if sol_ok {
                sol.clone()
            } else {
                // sol drifted too far (e.g. after ShakeDestroy); create a
                // diverse parent from best via perturbation + LNS.
                let mut pb = best.clone();
                #[allow(
                    clippy::cast_possible_truncation,
                    clippy::cast_sign_loss,
                    clippy::cast_precision_loss
                )]
                let n_moves = ((inst.m as f64 * hp.perturb_frac_hi).ceil() as usize)
                    .max(hp.perturb_min_moves);
                ges::perturb_solution(inst, &mut pb, n_moves, &mut info_buf, &mut rng);
                local_search::local_search(inst, &mut pb, &disabled_penalty, None);
                lns::lns(
                    inst,
                    &mut pb,
                    &mut best,
                    Duration::from_millis(200),
                    &mut rng,
                    deadline,
                    None,
                    None,
                    None,
                    None,
                    None,
                );
                if hp.op_log {
                    eprintln!(
                        "  [iter {}] CREX parent_b generated from best+perturb (veh={}, dist={:.2})",
                        iteration,
                        pb.num_vehicles(),
                        pb.total_distance(inst)
                    );
                }
                pb
            };

            last_crex_iter = iteration;
            let crex_seed = seed
                ^ ((iteration as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15))
                ^ (u64::from(diversify_stag) << 32);
            let mut crex_rng = SmallRng::seed_from_u64(crex_seed);
            // EAX with diverse elite parent; fall back to LCS-SREX with standard parent.
            let eax_child = if !elite_pool.is_empty() {
                let mut eidx = crex_rng.random_range(0..elite_pool.len());
                if elite_pool.len() > 1 && elite_pool[eidx].dist == best_dist {
                    eidx = (eidx + 1) % elite_pool.len();
                }
                crossover::eax_crossover(inst, &best, &elite_pool[eidx].sol, &mut crex_rng)
            } else {
                crossover::eax_crossover(inst, &best, &parent_b, &mut crex_rng)
            };
            let crossover_child = eax_child
                .or_else(|| crossover::lcs_srex_crossover(inst, &best, &parent_b, &mut crex_rng));
            if let Some(mut child) = crossover_child {
                // Cheap filter first; only spend full LS if child is already promising.
                local_search::local_search(inst, &mut child, &disabled_penalty, None);
                if is_better_than_with_cached_best(&child, &best, best_dist, inst) {
                    local_search::local_search_full(inst, &mut child, &disabled_penalty, None);
                    if child.is_feasible(inst)
                        && is_better_than_with_cached_best(&child, &best, best_dist, inst)
                    {
                        best.clone_from(&child);
                        best_dist = best.total_distance(inst);
                        sol.clone_from(&best);
                        accept_state.reset();
                        alns_state.reset();
                        best_improved_this_iter = true;
                        main_improved_this_iter = true;
                        no_improve_count = 0;
                        eprintln!(
                            "  [iter {}] CREX improved: {} vehicles, dist={:.2} (t={:.1}s)",
                            iteration,
                            best.num_vehicles(),
                            best_dist,
                            start.elapsed().as_secs_f64()
                        );
                        try_insert_elite(&mut elite_pool, &best, best_dist, node_base, &hp);
                    } else if hp.op_log {
                        eprintln!(
                            "  [iter {}] CREX no best improvement after full LS (veh={}, dist={:.2})",
                            iteration,
                            child.num_vehicles(),
                            child.total_distance(inst),
                        );
                    }
                } else if hp.op_log {
                    eprintln!(
                        "  [iter {}] CREX rejected after quick LS (veh={}, dist={:.2})",
                        iteration,
                        child.num_vehicles(),
                        child.total_distance(inst),
                    );
                }
            } else if hp.op_log {
                eprintln!("  [iter {iteration}] CREX produced no feasible child");
            }
        }

        // Deep-stagnation perturbation: SREX/LCS-SREX with elite pool.
        // In distance mode, activate earlier (stag>=8) with smaller retry gap (4 iters).
        let dm_srex_trigger = if in_distance_mode {
            eff_srex_stagnation_trigger.min(8)
        } else {
            eff_srex_stagnation_trigger
        };
        let dm_srex_gap = if in_distance_mode {
            hp.srex_retry_gap_iters.min(4)
        } else {
            hp.srex_retry_gap_iters
        };
        if diversify_stag >= dm_srex_trigger
            && iteration.saturating_sub(last_srex_iter) >= dm_srex_gap
            && deadline.saturating_duration_since(Instant::now())
                >= Duration::from_secs(eff_srex_min_remaining_secs)
            && !elite_pool.is_empty()
        {
            last_srex_iter = iteration;
            let mut idx = rng.random_range(0..elite_pool.len());
            if elite_pool.len() > 1 && elite_pool[idx].dist == best_dist {
                idx = (idx + 1) % elite_pool.len();
            }
            let elite_parent = &elite_pool[idx].sol;
            let swap_routes = rng.random_range(hp.srex_swap_routes_min..=hp.srex_swap_routes_max);
            if hp.op_log {
                eprintln!(
                    "  [iter {}] TRY SREX-elite (dstag={}, swap_routes={}, elite_idx={}, pool={})",
                    iteration,
                    diversify_stag,
                    swap_routes,
                    idx,
                    elite_pool.len()
                );
            }

            if let Some(mut child) =
                crossover::srex_route_exchange(inst, &best, elite_parent, swap_routes, &mut rng)
            {
                let mut srex_infos: Vec<RouteInfo> = Vec::new();
                lns::regret_insertion(
                    inst,
                    &mut child,
                    4,
                    0.0,
                    true,
                    &mut rng,
                    &mut srex_infos,
                    0,
                    false,
                );
                local_search::local_search_full(
                    inst,
                    &mut child,
                    &disabled_penalty,
                    if hp.op_log { Some(&mut ls_stats) } else { None },
                );
                let child_dist = child.total_distance(inst);
                let child_is_better =
                    is_better_than_with_cached_best(&child, &best, best_dist, inst);

                // Always move to this new point to escape stagnation.
                sol = child.clone();
                try_insert_elite(&mut elite_pool, &child, child_dist, node_base, &hp);

                if child_is_better && child.is_feasible(inst) {
                    best = child;
                    best_dist = child_dist;
                    best_improved_this_iter = true;
                    main_improved_this_iter = true;
                    no_improve_count = 0;
                    eprintln!(
                        "  [iter {}] SREX-elite improved: {} vehicles, dist={:.2} (t={:.1}s, pool={})",
                        iteration,
                        best.num_vehicles(),
                        best_dist,
                        start.elapsed().as_secs_f64(),
                        elite_pool.len(),
                    );
                } else {
                    eprintln!(
                        "  [iter {}] SREX-elite perturb: veh={}, dist={:.2} (t={:.1}s, pool={})",
                        iteration,
                        sol.num_vehicles(),
                        child_dist,
                        start.elapsed().as_secs_f64(),
                        elite_pool.len(),
                    );
                }
            } else if hp.op_log {
                eprintln!("  [iter {iteration}] SREX-elite produced no feasible child");
            }
        }

        // Tabu kick: periodically retry while stagnation persists.
        if no_improve_count >= eff_tabu_stagnation_trigger
            && iteration.saturating_sub(last_tabu_iter) >= hp.tabu_retry_gap_iters
            && Instant::now() < deadline
        {
            last_tabu_iter = iteration;
            let tabu_start = Instant::now();
            let remaining = deadline.saturating_duration_since(tabu_start);
            let tabu_slice_ms = if no_improve_count >= 14 {
                hp.tabu_slice_ms_hi
            } else if no_improve_count >= 10 {
                hp.tabu_slice_ms_mid
            } else {
                hp.tabu_slice_ms_lo
            };
            let tabu_deadline = tabu_start + remaining.min(Duration::from_millis(tabu_slice_ms));
            if hp.op_log {
                eprintln!(
                    "  [iter {iteration}] TRY TABU (stag={no_improve_count}, slice={tabu_slice_ms}ms)"
                );
            }

            if lns::tabu_search_stagnation(inst, &mut sol, &mut best, &mut rng, tabu_deadline) {
                best_dist = best.total_distance(inst);
                best_improved_this_iter = true;
                main_improved_this_iter = true;
                no_improve_count = 0;
                eprintln!(
                    "  [iter {}] TABU improved: {} vehicles, dist={:.2} (t={:.1}s, tabu={}ms)",
                    iteration,
                    best.num_vehicles(),
                    best_dist,
                    start.elapsed().as_secs_f64(),
                    tabu_start.elapsed().as_millis(),
                );
            } else if hp.op_log {
                eprintln!(
                    "  [iter {}] TABU no best improvement (elapsed={}ms)",
                    iteration,
                    tabu_start.elapsed().as_millis()
                );
            }
        }

        // Infeasible-space tabu search (Goeke 2019): single-vertex operators
        // with self-adjusting penalty multipliers. Disabled by default — benchmarks
        // showed no improvement on 800-node PDPTW instances. Set INFEAS_TABU=1 to enable.
        if std::env::var("INFEAS_TABU").ok().is_some_and(|v| v == "1")
            && no_improve_count >= 12
            && (no_improve_count - 12).is_multiple_of(8)
            && Instant::now() < deadline
        {
            let itabu_start = Instant::now();
            let itabu_budget = if short_run { 500u64 } else { 1200u64 };
            if hp.op_log {
                eprintln!(
                    "  [iter {iteration}] TRY InfeasTabu (stag={no_improve_count}, budget={itabu_budget}ms)"
                );
            }
            let improved =
                infeasible_tabu::infeasible_tabu_search(&best, inst, &mut rng, itabu_budget);
            if improved.is_feasible(inst) && improved.is_better_than(&best, inst) {
                best = improved;
                best_dist = best.total_distance(inst);
                sol.clone_from(&best);
                accept_state.reset();
                alns_state.reset();
                best_improved_this_iter = true;
                main_improved_this_iter = true;
                no_improve_count = 0;
                try_insert_elite(&mut elite_pool, &best, best_dist, node_base, &hp);
                eprintln!(
                    "  [iter {}] InfeasTabu improved: {} vehicles, dist={:.2} (t={:.1}s, itabu={}ms)",
                    iteration,
                    best.num_vehicles(),
                    best_dist,
                    start.elapsed().as_secs_f64(),
                    itabu_start.elapsed().as_millis(),
                );
            } else if hp.op_log {
                eprintln!(
                    "  [iter {}] InfeasTabu no improvement (elapsed={}ms)",
                    iteration,
                    itabu_start.elapsed().as_millis(),
                );
            }
        }

        // VLSN: cyclic exchanges during stagnation (composite moves LS can't find).
        // Only useful with tight TW (few feasible edges → cheap graph).
        // Try once when best changes, and periodically during deep stagnation
        // since sol may have shifted via shaking/tabu giving a different graph.
        let vlsn_due = (best_dist - vlsn_tried_at_dist).abs() > 1e-6
            || (diversify_stag >= 8 && diversify_stag.is_multiple_of(4));
        if diversify_stag >= 4 && vlsn_due && Instant::now() < deadline {
            vlsn_tried_at_dist = best_dist;
            let vlsn_start = Instant::now();
            if hp.op_log {
                eprintln!(
                    "  [iter {iteration}] TRY VLSN (dstag={diversify_stag}, best_dist={best_dist:.2})"
                );
            }
            let vlsn_improved = vlsn::vlsn(inst, &mut sol);
            let vlsn_ms = vlsn_start.elapsed().as_millis();
            if vlsn_improved {
                local_search::local_search(inst, &mut sol, &disabled_penalty, None);
                if is_better_than_with_cached_best(&sol, &best, best_dist, inst) {
                    best.clone_from(&sol);
                    best_dist = best.total_distance(inst);
                    best_improved_this_iter = true;
                    main_improved_this_iter = true;
                    no_improve_count = 0;
                    vlsn_tried_at_dist = f64::MAX; // reset so we try again on new best
                    eprintln!(
                        "  [iter {}] VLSN improved: {} vehicles, dist={:.2} (t={:.1}s, vlsn={}ms)",
                        iteration,
                        best.num_vehicles(),
                        best_dist,
                        start.elapsed().as_secs_f64(),
                        vlsn_ms,
                    );
                } else if hp.op_log {
                    eprintln!(
                        "  [iter {}] VLSN changed solution but no best improvement (vlsn={}ms, veh={}, dist={:.2})",
                        iteration,
                        vlsn_ms,
                        sol.num_vehicles(),
                        sol.total_distance(inst),
                    );
                }
            } else if hp.op_log {
                eprintln!("  [iter {iteration}] VLSN no move (vlsn={vlsn_ms}ms)");
            }
        }

        // Collect routes for set covering pool (only from feasible solutions
        // to prevent infeasible routes from polluting the pool).
        if sol.is_feasible(inst) {
            route_pool.add_solution(inst, &sol);
        }
        // Also feed best solution routes after each improvement for better pool diversity.
        if best_improved_this_iter {
            route_pool.add_solution(inst, &best);
            arc_history.add_solution(&best);
            // Decay arc frequencies every 20 solutions to prevent old patterns dominating.
            if arc_history.num_solutions() >= 20 {
                arc_history.decay();
            }
        }

        // Solve set covering MIP when pool has new routes.
        // In distance mode, use relaxed thresholds: more frequent SC solves
        // exploit the richer pool from LNS candidate feeding.
        let eff_sc_min_new = if in_distance_mode {
            hp.sc_min_new_routes.min(16)
        } else {
            hp.sc_min_new_routes
        };
        let eff_sc_max_iters = if in_distance_mode {
            hp.sc_max_iters_between_solves.min(3)
        } else {
            hp.sc_max_iters_between_solves
        };
        let sc_pending = route_pool.new_routes_since_solve();
        let sc_due_to_batch = sc_pending >= eff_sc_min_new;
        let sc_due_to_staleness = iteration.saturating_sub(last_sc_iter) >= eff_sc_max_iters;
        let sc_due_to_stagnation = diversify_stag >= hp.sc_stagnation_trigger;
        if route_pool.has_new_routes()
            && (sc_due_to_batch || sc_due_to_staleness || sc_due_to_stagnation)
            && Instant::now() < deadline
        {
            last_sc_iter = iteration;
            if let Some(mut sc_sol) = route_pool.solve(inst, &best) {
                // Full LS on SC result: SC assembles pool routes that may not be
                // locally optimal together — full LS extracts maximum improvement.
                // SC fires rarely so the extra cost of full LS is amortized.
                local_search::local_search_full(inst, &mut sc_sol, &disabled_penalty, None);
                // Feed SC+LS routes back to pool: LS may have improved
                // individual routes, making them useful for future SC solves.
                if sc_sol.is_feasible(inst) {
                    route_pool.add_solution(inst, &sc_sol);
                }
                if !sc_sol.is_feasible(inst) {
                    eprintln!("  SC: infeasible solution discarded");
                } else if is_better_than_with_cached_best(&sc_sol, &best, best_dist, inst) {
                    best.clone_from(&sc_sol);
                    best_dist = best.total_distance(inst);
                    sol.clone_from(&best);
                    accept_state.reset();
                    alns_state.reset();
                    best_improved_this_iter = true;
                    main_improved_this_iter = true;
                    try_insert_elite(&mut elite_pool, &best, best_dist, node_base, &hp);
                    eprintln!(
                        "  [iter {}] SC+LS improved: {} vehicles, dist={:.2} (t={:.1}s, pool={})",
                        iteration,
                        best.num_vehicles(),
                        best_dist,
                        start.elapsed().as_secs_f64(),
                        route_pool.len(),
                    );
                }
            }
        }

        // Periodically apply full local search (including 2-opt, or-opt, 2k-opt)
        // to the best solution. For instances with long routes (loose TW), these
        // intra-route operators find distance improvements that the regular LS misses.
        // Only worthwhile when routes are long (few vehicles) — for many-vehicle
        // solutions the routes are short and regular LS already optimizes them well.
        if iteration.is_multiple_of(5) && Instant::now() < deadline {
            let full_ls_start = Instant::now();
            let mut full_ls_candidate = best.clone();
            local_search::local_search_full(inst, &mut full_ls_candidate, &disabled_penalty, None);
            if full_ls_candidate.is_feasible(inst)
                && is_better_than_with_cached_best(&full_ls_candidate, &best, best_dist, inst)
            {
                best.clone_from(&full_ls_candidate);
                best_dist = best.total_distance(inst);
                // Don't reset sol or LAHC: let LNS continue its current
                // exploration trajectory. Only update best for reference.
                best_improved_this_iter = true;
                main_improved_this_iter = true;
                if hp.op_log {
                    eprintln!(
                        "  [iter {}] Periodic full LS improved: dist={:.2} ({}ms)",
                        iteration,
                        best_dist,
                        full_ls_start.elapsed().as_millis(),
                    );
                }
            }
        }

        // Keep elite pool updated with latest best, including LS/LNS/VLSN improvements.
        try_insert_elite(&mut elite_pool, &best, best_dist, node_base, &hp);

        // Update acceptance criterion adaptive weights
        accept_state.record_outcome(best_improved_this_iter);

        // Update stagnation counter
        if best_improved_this_iter {
            no_improve_count = 0;
        } else {
            no_improve_count += 1;
        }
        // Diversification stagnation: only reset when main-thread operations
        // (LS, LNS, SC, Shake, CREX, SREX, Tabu) find improvements. Worker-thread
        // improvements don't reset this, ensuring diversification operators fire.
        if main_improved_this_iter {
            diversify_stag = 0;
        } else {
            diversify_stag += 1;
        }

        // Track stable vehicle count for distance-mode tuning
        let current_veh = best.num_vehicles();
        if current_veh == last_vehicle_count {
            stable_vehicle_iters += 1;
        } else {
            stable_vehicle_iters = 0;
            last_vehicle_count = current_veh;
        }

        // (Removed: partial ruin-recreate and diverse construction approaches
        // both waste compute without improving results. LNS with 50% destruction
        // in distance mode is more effective for structural reorganization.)

        if hp.op_log && ls_stats.total_moves() > 0 {
            eprintln!("  [iter {}] LS: {}", iteration, ls_stats.summary_line());
        }

        if iteration.is_multiple_of(10) {
            eprintln!(
                "  [iter {}] best: {} vehicles, dist={:.2} (t={:.1}s, stag={}, accept={})",
                iteration,
                best.num_vehicles(),
                best_dist,
                start.elapsed().as_secs_f64(),
                no_improve_count,
                accept_state.weights_str(),
            );
        }
    }

    // Final SP solve with extended time limit
    if let Some(mut sc_sol) = route_pool.solve_final(inst, &best) {
        local_search::local_search_full(inst, &mut sc_sol, &disabled_penalty, None);
        if !sc_sol.is_feasible(inst) {
            eprintln!("  Final SC: infeasible solution discarded");
        } else if is_better_than_with_cached_best(&sc_sol, &best, best_dist, inst) {
            best.clone_from(&sc_sol);
            best_dist = best.total_distance(inst);
            eprintln!(
                "  Final SC+LS improved: {} vehicles, dist={:.2}",
                best.num_vehicles(),
                best_dist,
            );
        }
    }

    eprintln!(
        "  Final: {} vehicles, dist={:.2}, {} iterations in {:.1}s",
        best.num_vehicles(),
        best_dist,
        iteration,
        start.elapsed().as_secs_f64()
    );
    eprintln!(
        "  GES_TOTAL: runs={}, inner_iters={}, route_attempts={}, direct_calls={}",
        ges_total_runs,
        ges_total_stats.inner_iterations,
        ges_total_stats.route_attempts,
        ges_total_stats.direct_insertion_calls
    );

    best
}

#[cfg(test)]
mod tests {
    use super::InitialSolutionMethod;

    #[test]
    fn initial_method_parse_aliases() {
        assert_eq!(
            InitialSolutionMethod::from_str("legacy"),
            Some(InitialSolutionMethod::Legacy)
        );
        assert_eq!(
            InitialSolutionMethod::from_str("old"),
            Some(InitialSolutionMethod::Legacy)
        );
        assert_eq!(
            InitialSolutionMethod::from_str("lu2006"),
            Some(InitialSolutionMethod::LuDessouky2006)
        );
        assert_eq!(
            InitialSolutionMethod::from_str("ropke2006"),
            Some(InitialSolutionMethod::RopkePisinger2006)
        );
        assert_eq!(
            InitialSolutionMethod::from_str("cluster"),
            Some(InitialSolutionMethod::Cluster2004)
        );
        assert_eq!(
            InitialSolutionMethod::from_str("hosny2012"),
            Some(InitialSolutionMethod::Hosny2012)
        );
        assert_eq!(InitialSolutionMethod::from_str("unknown"), None);
    }
}