pdp_lns 0.1.1

use crate::instance::{Instance, arc_bit};
use crate::solution::{
    RouteInfo, Solution, best_pair_insertion_penalized, best_pair_insertion_with_info,
    build_route_with_pair, build_route_with_pair_buf, compute_all_infos, compute_all_infos_light,
    compute_info_with_pair, compute_info_without_pair, first_feasible_insertion_with_info, fmax,
    fmin, is_route_feasible, is_single_pair_feasible, route_without_pair, route_without_pair_buf,
};
use rand::RngExt;
use rand::prelude::*;
use std::time::{Duration, Instant};

/// Fast range reduction using Lemire's nearly-divisionless method.
/// Returns a value in [0, range) with negligible bias (< 1/2^32).
/// ~4× faster than `rng.random_range()` which goes through non-inlined trait dispatch.
#[inline(always)]
fn fast_range(rng: &mut impl Rng, range: usize) -> usize {
    debug_assert!(range > 0);
    let x = rng.next_u32();
    (u64::from(x).wrapping_mul(range as u64) >> 32) as usize
}

/// Fast coin flip: single bit test instead of generic random_bool dispatch.
#[inline(always)]
fn fast_coin(rng: &mut impl Rng) -> bool {
    rng.next_u32() & 1 != 0
}

/// Fast modulo for values in [0, 2*n). Avoids expensive hardware division.
#[inline(always)]
fn fast_mod(val: usize, n: usize) -> usize {
    if val >= n { val - n } else { val }
}

/// GES ejection mode: "bf" = brute-force (git baseline), "dfs" = lexicographic DFS (optimized).
/// Controlled by GES_MODE environment variable, default = "dfs".
fn ges_use_dfs() -> bool {
    !matches!(std::env::var("GES_MODE"), Ok(v) if v.eq_ignore_ascii_case("bf"))
}

fn ges_profile_enabled() -> bool {
    matches!(std::env::var("GES_PROFILE"),
        Ok(v) if v == "1" || v.eq_ignore_ascii_case("true") || v.eq_ignore_ascii_case("yes"))
}

fn ges_use_easiest_pick() -> bool {
    !matches!(std::env::var("GES_PICK_MIX"),
        Ok(v) if v == "0" || v.eq_ignore_ascii_case("false") || v.eq_ignore_ascii_case("no"))
}

fn ges_use_conflict_masks() -> bool {
    !matches!(std::env::var("GES_CONFLICT_MASKS"),
        Ok(v) if v == "0" || v.eq_ignore_ascii_case("false") || v.eq_ignore_ascii_case("no"))
}

fn ges_use_penalty_decay() -> bool {
    matches!(std::env::var("GES_PENALTY_DECAY"),
        Ok(v) if v == "1" || v.eq_ignore_ascii_case("true") || v.eq_ignore_ascii_case("yes"))
}

fn ges_use_dynamic_k() -> bool {
    !matches!(std::env::var("GES_DYNAMIC_K"),
        Ok(v) if v == "0" || v.eq_ignore_ascii_case("false") || v.eq_ignore_ascii_case("no"))
}

fn ges_use_cross_route_chain() -> bool {
    !matches!(std::env::var("GES_CROSS_ROUTE_CHAIN"),
        Ok(v) if v == "0" || v.eq_ignore_ascii_case("false") || v.eq_ignore_ascii_case("no"))
}

fn ges_use_temporal_perturb() -> bool {
    !matches!(std::env::var("GES_TEMPORAL_PERTURB"),
        Ok(v) if v == "0" || v.eq_ignore_ascii_case("false") || v.eq_ignore_ascii_case("no"))
}

pub fn ges_use_starvation() -> bool {
    !matches!(std::env::var("GES_STARVATION"),
        Ok(v) if v == "0" || v.eq_ignore_ascii_case("false") || v.eq_ignore_ascii_case("no"))
}

pub fn ges_use_sc_duals() -> bool {
    !matches!(std::env::var("GES_SC_DUALS"),
        Ok(v) if v == "0" || v.eq_ignore_ascii_case("false") || v.eq_ignore_ascii_case("no"))
}

fn ges_starvation_target() -> usize {
    std::env::var("GES_STARVATION_TARGET")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(5)
}

fn ges_use_infeasible() -> bool {
    !matches!(std::env::var("GES_INFEASIBLE"),
        Ok(v) if v == "0" || v.eq_ignore_ascii_case("false") || v.eq_ignore_ascii_case("no"))
}
#[derive(Clone, Copy, PartialEq, Eq)]
enum RouteSelectionMode {
    ShortTop3,
    SpatialTemporalTop3,
}

fn ges_route_selection_mode() -> RouteSelectionMode {
    match std::env::var("GES_ROUTE_SELECT") {
        Ok(v) if v.eq_ignore_ascii_case("short") => RouteSelectionMode::ShortTop3,
        Ok(v)
            if v.eq_ignore_ascii_case("st")
                || v.eq_ignore_ascii_case("spatiotemporal")
                || v.eq_ignore_ascii_case("hybrid")
                || v.eq_ignore_ascii_case("overlap") =>
        {
            RouteSelectionMode::SpatialTemporalTop3
        }
        _ => RouteSelectionMode::SpatialTemporalTop3,
    }
}

#[inline]
fn route_idle_ratio(inst: &Instance, route: &[usize], info: &RouteInfo) -> f64 {
    if route.len() <= 2 || info.arrival.len() != route.len() {
        return 0.0;
    }
    let mut svc_sum = 0.0;
    for &node in &route[..route.len() - 1] {
        svc_sum += inst.svc(node);
    }
    let duration = (info.arrival[route.len() - 1] - info.arrival[0]).max(0.0);
    if duration <= 1e-9 {
        return 0.0;
    }
    let moving = info.distance + svc_sum;
    ((duration - moving).max(0.0) / duration).clamp(0.0, 1.0)
}

#[inline]
fn temporal_overlap_ratio(a0: f64, a1: f64, b0: f64, b1: f64) -> f64 {
    let lo = fmax(a0.min(a1), b0.min(b1));
    let hi = fmin(a0.max(a1), b0.max(b1));
    let inter = (hi - lo).max(0.0);
    let uni_lo = a0.min(a1).min(b0.min(b1));
    let uni_hi = a0.max(a1).max(b0.max(b1));
    let union = (uni_hi - uni_lo).max(1e-9);
    (inter / union).clamp(0.0, 1.0)
}

#[inline]
fn route_anchor_min_dist(inst: &Instance, ra: &[usize], rb: &[usize]) -> f64 {
    if ra.len() <= 2 || rb.len() <= 2 {
        return 0.0;
    }
    let a_first = 1;
    let a_last = ra.len() - 2;
    let a_mid = usize::midpoint(a_first, a_last);
    let b_first = 1;
    let b_last = rb.len() - 2;
    let b_mid = usize::midpoint(b_first, b_last);
    let a_nodes = [ra[a_first], ra[a_mid], ra[a_last]];
    let b_nodes = [rb[b_first], rb[b_mid], rb[b_last]];
    let mut best = f64::MAX;
    for &u in &a_nodes {
        for &v in &b_nodes {
            best = fmin(best, inst.dist(u, v));
        }
    }
    best
}

#[inline]
fn effective_k_limit(unassigned_len: usize, consec_failures: u64, dynamic_k: bool) -> usize {
    if !dynamic_k {
        return K_BASE;
    }
    if unassigned_len <= 1 {
        return K_HARD_MAX;
    }
    if unassigned_len <= 2 {
        return (K_BASE + 1).min(K_HARD_MAX);
    }
    if unassigned_len <= 4 && consec_failures >= 96 {
        return (K_BASE + 1).min(K_HARD_MAX);
    }
    K_BASE
}

#[derive(Default)]
struct GesProfile {
    attempts: usize,
    direct_success: u64,
    direct_fail: u64,
    ejection_success: u64,
    ejection_fail: u64,
    squeeze_attempts: u64,
    squeeze_success: u64,
    perturb_calls: u64,
    infos_recompute: u64,
    t_direct: Duration,
    t_ejection: Duration,
    t_squeeze: Duration,
    t_perturb: Duration,
    t_infos: Duration,
    repair_attempts: u64,
    repair_successes: u64,
    infeasible_entries: u64,
    infeasible_repair_ok: u64,
    t_repair: Duration,
}

#[derive(Default, Clone, Debug)]
pub struct GesRunStats {
    pub inner_iterations: u64,
    pub route_attempts: u64,
    pub direct_insertion_calls: u64,
    pub direct_successes: u64,
    pub direct_failures: u64,
    pub direct_route_checks: u64,
    pub direct_cache_hits: u64,
    pub direct_arc_tp_rejects: u64,
    pub direct_arc_d_rejects: u64,
    pub direct_arc_td_rejects: u64,
    pub direct_conflict_rejects: u64,
    pub ejection_calls: u64,
    pub ejection_successes: u64,
    pub squeeze_attempts: u64,
    pub squeeze_successes: u64,
    pub perturb_calls: u64,
    pub pickup_easiest_picks: u64,
    pub pickup_hardest_picks: u64,
    pub infeasible_entries: u64,
    pub infeasible_repair_successes: u64,
    pub infeasible_insertions: u64,
}

impl GesRunStats {
    pub fn merge_from(&mut self, other: &Self) {
        self.inner_iterations += other.inner_iterations;
        self.route_attempts += other.route_attempts;
        self.direct_insertion_calls += other.direct_insertion_calls;
        self.direct_successes += other.direct_successes;
        self.direct_failures += other.direct_failures;
        self.direct_route_checks += other.direct_route_checks;
        self.direct_cache_hits += other.direct_cache_hits;
        self.direct_arc_tp_rejects += other.direct_arc_tp_rejects;
        self.direct_arc_d_rejects += other.direct_arc_d_rejects;
        self.direct_arc_td_rejects += other.direct_arc_td_rejects;
        self.direct_conflict_rejects += other.direct_conflict_rejects;
        self.ejection_calls += other.ejection_calls;
        self.ejection_successes += other.ejection_successes;
        self.squeeze_attempts += other.squeeze_attempts;
        self.squeeze_successes += other.squeeze_successes;
        self.perturb_calls += other.perturb_calls;
        self.pickup_easiest_picks += other.pickup_easiest_picks;
        self.pickup_hardest_picks += other.pickup_hardest_picks;
        self.infeasible_entries += other.infeasible_entries;
        self.infeasible_repair_successes += other.infeasible_repair_successes;
        self.infeasible_insertions += other.infeasible_insertions;
    }
}

impl GesProfile {
    fn print(&self, elapsed: Duration, routes_before: usize, routes_after: usize, removed: bool) {
        eprintln!(
            "[GES_PROFILE] t={:.3}s routes:{}->{} removed={} attempts={} direct(ok/fail)={}/{} ejection(ok/fail)={}/{} squeeze(ok/try)={}/{} perturb={} infos_recompute={} repair(ok/try)={}/{} infeasible(entries/repair_ok)={}/{} time(direct/ejection/squeeze/perturb/infos/repair)={:.3}/{:.3}/{:.3}/{:.3}/{:.3}/{:.3}s",
            elapsed.as_secs_f64(),
            routes_before,
            routes_after,
            removed,
            self.attempts,
            self.direct_success,
            self.direct_fail,
            self.ejection_success,
            self.ejection_fail,
            self.squeeze_success,
            self.squeeze_attempts,
            self.perturb_calls,
            self.infos_recompute,
            self.repair_successes,
            self.repair_attempts,
            self.infeasible_entries,
            self.infeasible_repair_ok,
            self.t_direct.as_secs_f64(),
            self.t_ejection.as_secs_f64(),
            self.t_squeeze.as_secs_f64(),
            self.t_perturb.as_secs_f64(),
            self.t_infos.as_secs_f64(),
            self.t_repair.as_secs_f64(),
        );
    }
}

/// Maximum perturbation counter for the first GES attempt (Section 3.2, Table 1).
const MAX_PERTURBS: u64 = 1_250_000;
/// Reduced perturbation budget for subsequent attempts after a successful vehicle reduction.
const MAX_PERTURBS_SUBSEQUENT: u64 = 350_000;
/// Number of random moves per perturbation within GES (paper I_rand = 150).
const PERTURBS_WITHIN_GES: usize = 206;
/// Shuffle the unassigned stack after this many consecutive perturbations without progress.
const SHUFFLE_INTERVAL: u64 = 25;
/// Try a different squeeze target every this many consecutive failures.
const SQUEEZE_INTERVAL: u64 = 18;
/// Apply squeeze only near completion (many early squeezes are wasted).
const SQUEEZE_MAX_UNASSIGNED: usize = 6;
/// During long stalls, periodically switch from hardest-first to easiest-first
/// to quickly empty the stack and unlock a different feasibility basin.
const EASIEST_PICK_PERIOD: u64 = 28;
const EASIEST_PICK_STALL: u64 = 104;
/// Multiplicative penalty decay applied periodically during stagnation.
const PENALTY_DECAY_NUM: u64 = 15;
const PENALTY_DECAY_DEN: u64 = 22;
/// Infeasible space search: enter after this many consecutive failures.
const INFEASIBLE_THRESHOLD: u64 = 150;
/// Only enter infeasible mode when few pairs remain.
const INFEASIBLE_MAX_UNASSIGNED: usize = 5;
/// Initial penalty for TW and capacity in infeasible insertion mode.
const INFEASIBLE_PENALTY_INIT: f64 = 100.0;
/// Max iterations for squeeze targeted repair.
const REPAIR_MAX_ITERS: usize = 60;
/// Max iterations for post-infeasible-mode repair.
const REPAIR_MAX_ITERS_FULL: usize = 120;

/// Adaptive perturbation size inside one GES attempt.
/// Goal: keep strong diversification near success, but reduce wasted churn on long stalls.
fn adaptive_perturb_moves(unassigned_len: usize, consec_failures: u64) -> usize {
    if unassigned_len <= 2 {
        // Near completion, strong perturbation helps unlock the last pairs.
        return 3 * PERTURBS_WITHIN_GES;
    }
    if unassigned_len <= 4 {
        return PERTURBS_WITHIN_GES;
    }
    // During deep stalls, shrink move count to increase state transitions/sec.
    // Keep periodic full-size pulses for diversification.
    if consec_failures.is_multiple_of(64) {
        return PERTURBS_WITHIN_GES;
    }
    match consec_failures {
        0..=31 => PERTURBS_WITHIN_GES,
        32..=127 => 96,
        128..=511 => 64,
        512..=2047 => 40,
        _ => 24,
    }
}

/// Select next pickup from unassigned stack.
/// Default: hardest-first (max penalty). Under deep stagnation: periodic easiest-first.
fn select_unassigned_idx(
    unassigned_stack: &[usize],
    penalties: &[u64],
    consec_failures: u64,
    allow_easiest: bool,
) -> (usize, bool) {
    let want_easiest = allow_easiest
        && (unassigned_stack.len() <= 2
            || (consec_failures >= EASIEST_PICK_STALL
                && consec_failures.is_multiple_of(EASIEST_PICK_PERIOD)));
    let mut best_idx = unassigned_stack.len() - 1;
    let mut best_pen = penalties[unassigned_stack[best_idx]];
    for i in 0..unassigned_stack.len() - 1 {
        let p = penalties[unassigned_stack[i]];
        let better = if want_easiest {
            p < best_pen
        } else {
            p > best_pen
        };
        if better {
            best_pen = p;
            best_idx = i;
        }
    }
    (best_idx, want_easiest)
}

/// Pre-allocated buffers for ejection search and perturbation, reused across calls.
struct EjSearchBufs {
    rwp_buf: Vec<usize>,
    rwp_buf2: Vec<usize>,
    info_buf: RouteInfo,
    z_buf: Vec<f64>,
    smp_buf: Vec<u64>,
    rollback_routes: Vec<Vec<usize>>,
    squeeze_routes: Vec<Vec<usize>>,
    /// Reusable buffers for perturbation (avoids per-call allocation).
    perturb_buf1: Vec<usize>,
    perturb_buf2: Vec<usize>,
    perturb_route_buf: Vec<usize>,
    /// Reusable buffers for ejection search (avoids per-call allocation).
    route_conflicts_buf: Vec<u8>,
    phase1_cands_buf: Vec<(u64, usize, usize)>,
    ph2_routes_buf: Vec<(usize, u64, u64)>,
    /// Reusable buffer for route selection in GES outer loop.
    route_indices_buf: Vec<usize>,
    route_scores_buf: Vec<(usize, f64, f64, f64, f64)>, // (route_idx, overlap, idle, short, sc_raw)
    pair_scores_buf: Vec<f64>,
}

impl EjSearchBufs {
    fn new() -> Self {
        Self {
            rwp_buf: Vec::new(),
            rwp_buf2: Vec::new(),
            info_buf: RouteInfo::new(),
            z_buf: Vec::new(),
            smp_buf: Vec::new(),
            rollback_routes: Vec::new(),
            squeeze_routes: Vec::new(),
            perturb_buf1: Vec::new(),
            perturb_buf2: Vec::new(),
            perturb_route_buf: Vec::new(),
            route_conflicts_buf: Vec::new(),
            phase1_cands_buf: Vec::new(),
            ph2_routes_buf: Vec::new(),
            route_indices_buf: Vec::new(),
            route_scores_buf: Vec::new(),
            pair_scores_buf: Vec::new(),
        }
    }
}

/// Random-start circular scan for a pickup in a route. O(1) average since ~50%
/// of non-depot nodes are pickups. Returns 0 if no pickups found.
fn sample_pickup(inst: &Instance, route: &[usize], rng: &mut impl Rng) -> usize {
    let inner = route.len().saturating_sub(2); // exclude both depots
    if inner == 0 {
        return 0;
    }
    let start = fast_range(rng, inner);
    for i in 0..inner {
        let v = route[1 + fast_mod(start + i, inner)];
        if v != 0 && inst.is_pickup(v) {
            return v;
        }
    }
    0
}

/// Adaptive Guided Ejection Search.
/// Tries to reduce the number of vehicles in the solution.
/// Returns true if at least one route was successfully removed.
/// (Section 3.2 of the paper, Fig. 2)
///
/// After successfully removing a route, loops back to pick another route
/// and tries to remove it too (Fig. 2, line 10: "go to 1."). The perturbation
/// counter is shared across all attempts.
pub fn guided_ejection_search(
    inst: &Instance,
    sol: &mut Solution,
    rng: &mut impl Rng,
    deadline: Instant,
    stats_out: Option<&mut GesRunStats>,
    sc_duals: Option<&[f64]>,
) -> bool {
    let fn_start = Instant::now();
    let routes_before = sol.routes.len();
    let profile_enabled = ges_profile_enabled();
    let mut profile = GesProfile::default();
    let mut any_removed = false;
    let mut run_stats = GesRunStats::default();
    let use_pick_mix = ges_use_easiest_pick();
    let use_conflict_masks = ges_use_conflict_masks();
    let use_penalty_decay = ges_use_penalty_decay();
    let use_dynamic_k = ges_use_dynamic_k();
    let use_cross_route_chain = ges_use_cross_route_chain();
    let route_select_mode = ges_route_selection_mode();
    // 50/50 coin flip: use SC dual scoring or original heuristic
    let use_sc = sc_duals.is_some() && fast_coin(rng);

    // Pre-allocated buffers for ejection search — avoids per-call allocations
    let mut bufs = EjSearchBufs::new();

    // Pre-allocated buffer for direct insertion (Step 5)
    let mut route_buf: Vec<usize> = Vec::new();

    // Cached RouteInfo for all routes — avoids recomputing unchanged routes
    let mut infos: Vec<RouteInfo> = Vec::new();
    let mut infos_dirty: bool = true;
    // True when infos lack prefix_tw/suffix_tw (after compute_light in perturbation).
    // Must be resolved before squeeze_insertion which needs penalized insertion.
    let mut infos_tw_dirty: bool = false;
    // Cached bit-mask of requests present in each route for O(stride) conflict checks.
    let mut route_req_masks: Vec<Vec<u64>> = Vec::new();
    let mut route_masks_dirty: bool = true;

    // Penalty counters per PD pair (paper: p[h] := 1, incremented on failed direct insertion)
    let mut penalties: Vec<u64> = vec![1; inst.n + 1];

    // Track which routes have been tried (by original index) to avoid repeating
    let mut tried_routes: Vec<bool> = vec![false; sol.routes.len()];
    let mut attempts: usize = 0;

    // Versioned memoization: track per-route version and cache failed insertions
    let mut route_versions: Vec<u64> = Vec::new();
    let mut failed_inserts: Vec<Vec<u64>> = Vec::new();

    // Outer loop: after successfully removing a route, try another (Fig. 2, line 10)
    // On failure: try a different route instead of immediately giving up.
    loop {
        if sol.routes.len() <= 1 || Instant::now() >= deadline {
            break;
        }

        // Per-attempt perturbation limit: full budget for first attempt,
        // reduced budget after a successful reduction (second attempt usually fails).
        let max_perturbs = if any_removed {
            MAX_PERTURBS_SUBSEQUENT
        } else if attempts > 0 {
            MAX_PERTURBS / 3
        } else {
            MAX_PERTURBS
        };
        let mut perturb_counter: u64 = 0;

        // Snapshot before route removal — exact rollback on failure (Step 1)
        bufs.rollback_routes.clone_from(&sol.routes);
        let rollback_tried = tried_routes.clone();

        // Step 1: Select a route to remove and un-assign all PD-pairs from it.
        // Modes:
        // - short: original top-3 shortest routes
        // - st (default): top-3 by spatial/temporal overlap + idle + shortness
        let route_idx = {
            let nr = sol.routes.len();
            bufs.route_indices_buf.clear();
            for (i, &tried) in tried_routes.iter().enumerate().take(nr) {
                if !tried {
                    bufs.route_indices_buf.push(i);
                }
            }
            if bufs.route_indices_buf.is_empty() {
                break;
            }
            if route_select_mode == RouteSelectionMode::ShortTop3 {
                bufs.route_indices_buf
                    .sort_unstable_by_key(|&i| sol.routes[i].len());
                let top = bufs.route_indices_buf.len().min(3);
                bufs.route_indices_buf[fast_range(rng, top)]
            } else {
                if infos_dirty {
                    compute_all_infos_light(inst, &sol.routes, &mut infos);
                    infos_dirty = false;
                    infos_tw_dirty = true;
                }
                // Pass 1: compute raw overlap, idle, short, sc_raw for each candidate
                bufs.route_scores_buf.clear();
                for &i in &bufs.route_indices_buf {
                    let ri = &sol.routes[i];
                    let ii = &infos[i];
                    let idle = route_idle_ratio(inst, ri, ii);
                    let short = 1.0 / (ri.len().saturating_sub(2).max(1) as f64);
                    let a0 = if ri.len() > 2 { ii.arrival[1] } else { 0.0 };
                    let a1 = if ri.len() > 2 {
                        ii.arrival[ri.len() - 2]
                    } else {
                        0.0
                    };

                    bufs.pair_scores_buf.clear();
                    let scale_i = ii.distance / (ri.len().saturating_sub(1).max(1) as f64);
                    for (rj, ij) in infos.iter().enumerate() {
                        if rj == i {
                            continue;
                        }
                        let rj_route = &sol.routes[rj];
                        let b0 = if rj_route.len() > 2 {
                            ij.arrival[1]
                        } else {
                            0.0
                        };
                        let b1 = if rj_route.len() > 2 {
                            ij.arrival[rj_route.len() - 2]
                        } else {
                            0.0
                        };
                        let temporal = temporal_overlap_ratio(a0, a1, b0, b1);
                        let anchor_dist = route_anchor_min_dist(inst, ri, rj_route);
                        let scale_j =
                            ij.distance / (rj_route.len().saturating_sub(1).max(1) as f64);
                        let spatial = 1.0 / (1.0 + anchor_dist / (1.0 + 0.5 * (scale_i + scale_j)));
                        let st = 0.7 * spatial + 0.3 * temporal;
                        bufs.pair_scores_buf.push(st);
                    }
                    let overlap = if bufs.pair_scores_buf.is_empty() {
                        0.0
                    } else {
                        bufs.pair_scores_buf.sort_unstable_by(|a, b| b.total_cmp(a));
                        let k = bufs.pair_scores_buf.len().min(3);
                        let mut s = 0.0;
                        for &v in bufs.pair_scores_buf.iter().take(k) {
                            s += v;
                        }
                        s / k as f64
                    };
                    // Compute average dual value per pickup for route i
                    let sc_raw = if let Some(duals) = sc_duals {
                        let mut dual_sum = 0.0;
                        let mut np = 0usize;
                        for &v in ri {
                            if v != 0 && inst.is_pickup(v) {
                                dual_sum += duals[v];
                                np += 1;
                            }
                        }
                        if np > 0 { dual_sum / np as f64 } else { 0.0 }
                    } else {
                        0.0
                    };
                    bufs.route_scores_buf
                        .push((i, overlap, idle, short, sc_raw));
                }
                // Pass 2: normalize SC scores and compute blended score
                if use_sc {
                    let max_sc = bufs
                        .route_scores_buf
                        .iter()
                        .map(|e| e.4)
                        .fold(0.0_f64, f64::max);
                    // Sort by blended score with SC dual component
                    bufs.route_scores_buf.sort_unstable_by(|a, b| {
                        let sc_a = if max_sc > 1e-12 {
                            1.0 - (a.4 / max_sc).clamp(0.0, 1.0)
                        } else {
                            0.0
                        };
                        let sc_b = if max_sc > 1e-12 {
                            1.0 - (b.4 / max_sc).clamp(0.0, 1.0)
                        } else {
                            0.0
                        };
                        let score_a = 0.45 * a.1 + 0.2 * a.2 + 0.1 * a.3 + 0.25 * sc_a;
                        let score_b = 0.45 * b.1 + 0.2 * b.2 + 0.1 * b.3 + 0.25 * sc_b;
                        score_b.total_cmp(&score_a)
                    });
                } else {
                    // Original: 0.6 * overlap + 0.3 * idle + 0.1 * short
                    bufs.route_scores_buf.sort_unstable_by(|a, b| {
                        let score_a = 0.6 * a.1 + 0.3 * a.2 + 0.1 * a.3;
                        let score_b = 0.6 * b.1 + 0.3 * b.2 + 0.1 * b.3;
                        score_b.total_cmp(&score_a)
                    });
                }
                let top = bufs.route_scores_buf.len().min(3);
                bufs.route_scores_buf[fast_range(rng, top)].0
            }
        };
        tried_routes[route_idx] = true;
        attempts += 1;
        // swap_remove: O(1) instead of O(R) shift. Infos/masks at other indices
        // stay valid; the swapped-in element at route_idx keeps its correct info.
        let removed_route = sol.routes.swap_remove(route_idx);
        run_stats.route_attempts += 1;
        tried_routes.swap_remove(route_idx);
        if route_idx < infos.len() {
            infos.swap_remove(route_idx);
        }
        if use_conflict_masks && route_idx < route_req_masks.len() {
            route_req_masks.swap_remove(route_idx);
        }
        // No infos_dirty: swap_remove preserves index↔info correspondence.
        // Masks already correct after swap_remove.

        // Collect PD pairs from the removed route (as pickup IDs)
        let mut unassigned_stack: Vec<usize> = Vec::new();
        for &v in &removed_route {
            if v != 0 && inst.is_pickup(v) {
                unassigned_stack.push(v);
            }
        }
        unassigned_stack.shuffle(rng);

        let mut min_unassigned = unassigned_stack.len();
        let mut route_removed = false;
        let mut consec_failures: u64 = 0; // consecutive perturbations without progress
        // Initialize version cache for this attempt
        route_versions.clear();
        route_versions.resize(sol.routes.len(), 0);
        failed_inserts.clear();
        failed_inserts.resize_with(sol.routes.len(), || vec![u64::MAX; inst.n + 1]);

        // Revert-to-best: save solution + unassigned stack at each new minimum.
        // On long stalls, revert to the best-known state and try a different
        // perturbation trajectory instead of drifting in a poor region.
        let mut best_state_routes: Vec<Vec<usize>> = sol.routes.clone();
        let mut best_state_stack: Vec<usize> = unassigned_stack.clone();
        let mut reverts_remaining: u32 = 5; // limit reverts to avoid infinite loops
        const REVERT_INTERVAL: u64 = 320;

        let use_dfs = ges_use_dfs();

        // Infeasible space search state
        let use_infeasible = ges_use_infeasible();
        let mut infeasible_mode = false;
        let mut infeasible_insertions: u32 = 0;
        let mut infeasible_snapshot_routes: Vec<Vec<usize>> = Vec::new();
        let mut infeasible_snapshot_stack: Vec<usize> = Vec::new();
        let mut infeasible_pen: f64 = INFEASIBLE_PENALTY_INIT;

        // Batched deadline check: Instant::now() costs ~20ns (vDSO syscall).
        // At 57M iterations, checking every iteration wastes ~1.1s.
        // Check every 256 iterations instead (~0.005s granularity).
        let mut deadline_counter: u32 = 0;

        // Inner loop: try to reinsert all unassigned pairs (Algorithm 1)
        while !unassigned_stack.is_empty() && perturb_counter < max_perturbs {
            deadline_counter = deadline_counter.wrapping_add(1);
            if deadline_counter.trailing_zeros() >= 8 && Instant::now() >= deadline {
                break;
            }
            run_stats.inner_iterations += 1;
            // Ensure cached infos are up-to-date (light info suffices for
            // direct insertion + ejection; squeeze upgrades to full when needed).
            if infos_dirty {
                let t0 = if profile_enabled {
                    Some(Instant::now())
                } else {
                    None
                };
                compute_all_infos_light(inst, &sol.routes, &mut infos);
                if let Some(t0) = t0 {
                    profile.infos_recompute += 1;
                    profile.t_infos += t0.elapsed();
                }
                infos_dirty = false;
                infos_tw_dirty = true; // tw fields not computed
                route_masks_dirty = true;
            }
            if use_conflict_masks && route_masks_dirty {
                if route_req_masks.len() != sol.routes.len() {
                    route_req_masks.resize_with(sol.routes.len(), Vec::new);
                }
                for (ri, route) in sol.routes.iter().enumerate() {
                    inst.build_route_req_mask(route, &mut route_req_masks[ri]);
                }
                route_masks_dirty = false;
            }

            // Select hardest pair: pick the one with highest penalty (most failed insertions).
            // This prioritizes hard-to-insert pairs while routes still have capacity,
            // instead of LIFO which may waste good slots on easy pairs first.
            let (best_idx, easiest_pick) =
                select_unassigned_idx(&unassigned_stack, &penalties, consec_failures, use_pick_mix);
            if easiest_pick {
                run_stats.pickup_easiest_picks += 1;
            } else {
                run_stats.pickup_hardest_picks += 1;
            }
            let last = unassigned_stack.len() - 1;
            unassigned_stack.swap(best_idx, last);
            let pickup = *unassigned_stack.last().unwrap();
            let delivery = inst.delivery_of(pickup);

            // Try to insert the pair into any existing route (random order)
            // Partial Fisher-Yates: stop as soon as we find a feasible route.
            let t_direct0 = if profile_enabled {
                Some(Instant::now())
            } else {
                None
            };
            run_stats.direct_insertion_calls += 1;
            let mut insertion = None;
            let mut best_infeasible_cost = f64::MAX;
            let nr_order = sol.routes.len();
            let start = fast_range(rng, nr_order);
            for off in 0..nr_order {
                let ri = fast_mod(start + off, nr_order);
                run_stats.direct_route_checks += 1;
                // Version cache: skip route if insertion already failed for current version
                if ri < failed_inserts.len() && failed_inserts[ri][pickup] == route_versions[ri] {
                    run_stats.direct_cache_hits += 1;
                    continue;
                }
                // Conflict check: skip if inserting pair conflicts with any request in route
                let has_conflict = if use_conflict_masks {
                    run_stats.direct_cache_hits += 1;
                    inst.has_conflict_with_mask(&route_req_masks[ri], pickup)
                } else {
                    inst.has_route_conflict(&sol.routes[ri], pickup)
                };
                if has_conflict {
                    run_stats.direct_conflict_rejects += 1;
                    continue;
                }
                if !infeasible_mode {
                    if let Some((_, ep, ed)) = first_feasible_insertion_with_info(
                        inst,
                        &sol.routes[ri],
                        &infos[ri],
                        pickup,
                        delivery,
                    ) {
                        insertion = Some((ri, ep, ed));
                        break;
                    }
                    // Mark this insertion as failed for current route version
                    if ri < failed_inserts.len() {
                        failed_inserts[ri][pickup] = route_versions[ri];
                    }
                } else {
                    // Infeasible mode: penalized insertion, scan all routes for best
                    if let Some((cost, ep, ed)) = best_pair_insertion_penalized(
                        inst,
                        &sol.routes[ri],
                        &infos[ri],
                        pickup,
                        delivery,
                        infeasible_pen,
                        infeasible_pen,
                    ) && (insertion.is_none() || cost < best_infeasible_cost)
                    {
                        insertion = Some((ri, ep, ed));
                        best_infeasible_cost = cost;
                    }
                }
            }
            if let Some(t0) = t_direct0 {
                profile.t_direct += t0.elapsed();
            }

            if let Some((ri, ep, ed)) = insertion {
                // Direct insertion succeeded
                if infeasible_mode {
                    infeasible_insertions += 1;
                    run_stats.infeasible_insertions += 1;
                    // Full compute so infos stay complete for penalized insertion
                    build_route_with_pair_buf(
                        &sol.routes[ri],
                        pickup,
                        ep,
                        delivery,
                        ed,
                        &mut route_buf,
                    );
                    std::mem::swap(&mut sol.routes[ri], &mut route_buf);
                    infos[ri].compute(inst, &sol.routes[ri]);
                    infos_tw_dirty = false;
                    // Increase penalty to push future insertions toward less violation
                    infeasible_pen *= 1.3;
                } else {
                    // Fused: build route with pair + compute info in one pass
                    compute_info_with_pair(
                        inst,
                        &sol.routes[ri],
                        pickup,
                        ep,
                        delivery,
                        ed,
                        &mut route_buf,
                        &mut infos[ri],
                    );
                    std::mem::swap(&mut sol.routes[ri], &mut route_buf);
                    infos_tw_dirty = true;
                }
                if use_conflict_masks {
                    // O(1) bit set instead of O(n) mask rebuild: we only added one pair.
                    let req = inst.pickup_to_req[pickup];
                    if req != usize::MAX {
                        route_req_masks[ri][req >> 6] |= 1u64 << (req & 63);
                    }
                }
                if ri < route_versions.len() {
                    route_versions[ri] = route_versions[ri].wrapping_add(1);
                }
                unassigned_stack.pop();
                consec_failures = 0;

                // Check progress
                if unassigned_stack.len() < min_unassigned {
                    min_unassigned = unassigned_stack.len();
                    perturb_counter = 0;
                    // Save best state for potential revert
                    best_state_routes.clone_from(&sol.routes);
                    best_state_stack.clone_from(&unassigned_stack);
                }

                if unassigned_stack.is_empty() {
                    if infeasible_mode && infeasible_insertions > 0 {
                        // Some insertions were infeasible — try repair
                        let repaired = targeted_repair(
                            inst,
                            sol,
                            &mut infos,
                            &mut route_buf,
                            &mut bufs.info_buf,
                            rng,
                            REPAIR_MAX_ITERS_FULL,
                        );
                        if repaired {
                            route_removed = true;
                            run_stats.infeasible_repair_successes += 1;
                        } else {
                            // Repair failed — rollback to pre-infeasible state
                            sol.routes.clone_from(&infeasible_snapshot_routes);
                            unassigned_stack.clone_from(&infeasible_snapshot_stack);
                            unassigned_stack.shuffle(rng);
                            infeasible_mode = false;
                            infeasible_insertions = 0;
                            infos_dirty = true;
                            route_masks_dirty = true;
                            for v in &mut route_versions {
                                *v = v.wrapping_add(1);
                            }
                            consec_failures = 0;
                        }
                    } else {
                        route_removed = true;
                    }
                    if route_removed {
                        break;
                    }
                }
                if profile_enabled {
                    profile.direct_success += 1;
                }
                run_stats.direct_successes += 1;
                continue;
            }
            if profile_enabled {
                profile.direct_fail += 1;
            }
            run_stats.direct_failures += 1;

            // Direct insertion failed — increment penalty (paper Algorithm 1, line 10)
            penalties[pickup] += 1;

            // Try ejection search (paper Algorithm 1, line 11)
            let t_ej0 = if profile_enabled {
                Some(Instant::now())
            } else {
                None
            };
            run_stats.ejection_calls += 1;
            let k_limit = effective_k_limit(unassigned_stack.len(), consec_failures, use_dynamic_k);
            let ejection = if use_dfs {
                find_best_ejection_dfs(
                    inst,
                    sol,
                    pickup,
                    delivery,
                    &penalties,
                    &infos,
                    &mut bufs,
                    &route_req_masks,
                    k_limit,
                )
            } else {
                find_best_ejection_bf(inst, sol, pickup, delivery, &penalties, &mut bufs, k_limit)
            };
            if let Some(t0) = t_ej0 {
                profile.t_ejection += t0.elapsed();
            }

            if let Some((ri, new_route, ejected_pickups, num_ejected)) = ejection {
                sol.routes[ri] = new_route;
                infos[ri].compute_light(inst, &sol.routes[ri]);
                infos_tw_dirty = true;
                if ri < route_versions.len() {
                    route_versions[ri] = route_versions[ri].wrapping_add(1);
                }
                if use_conflict_masks {
                    // O(k) incremental update: add inserted pickup, remove ejected pickups.
                    let req = inst.pickup_to_req[pickup];
                    if req != usize::MAX {
                        route_req_masks[ri][req >> 6] |= 1u64 << (req & 63);
                    }
                    for &ep in &ejected_pickups[..num_ejected] {
                        let req = inst.pickup_to_req[ep];
                        if req != usize::MAX {
                            route_req_masks[ri][req >> 6] &= !(1u64 << (req & 63));
                        }
                    }
                }
                // Remove empty routes created by two-pair ejection from same route.
                // Use swap_remove to avoid O(R) shifts and keep infos/masks valid.
                let mut idx = 0usize;
                while idx < sol.routes.len() {
                    if sol.routes[idx].len() <= 2 {
                        sol.routes.swap_remove(idx);
                        if idx < infos.len() {
                            infos.swap_remove(idx);
                        }
                        if use_conflict_masks && idx < route_req_masks.len() {
                            route_req_masks.swap_remove(idx);
                        }
                        if idx < route_versions.len() {
                            route_versions.swap_remove(idx);
                        }
                        if idx < failed_inserts.len() {
                            failed_inserts.swap_remove(idx);
                        }
                        // Don't increment: check the swapped-in element
                    } else {
                        idx += 1;
                    }
                }
                unassigned_stack.pop(); // remove h_in
                for &ep in &ejected_pickups[..num_ejected] {
                    let inserted_elsewhere = use_cross_route_chain
                        && try_cross_route_insert(
                            inst,
                            sol,
                            ep,
                            &mut infos,
                            &mut route_buf,
                            &mut route_req_masks,
                            use_conflict_masks,
                            &mut route_versions,
                            rng,
                        );
                    if !inserted_elsewhere {
                        unassigned_stack.push(ep);
                    }
                }
                consec_failures = 0;
                if profile_enabled {
                    profile.ejection_success += 1;
                }
                run_stats.ejection_successes += 1;
            } else if consec_failures > 0
                && consec_failures.is_multiple_of(SQUEEZE_INTERVAL)
                && unassigned_stack.len() <= SQUEEZE_MAX_UNASSIGNED
            {
                if profile_enabled {
                    profile.ejection_fail += 1;
                }
                // Squeeze operation (cuopt-inspired): when ejection fails repeatedly,
                // force-insert the pair at the position minimizing total excess,
                // then perturb to try to repair feasibility.
                if infos_dirty || infos_tw_dirty {
                    let t0 = if profile_enabled {
                        Some(Instant::now())
                    } else {
                        None
                    };
                    compute_all_infos(inst, &sol.routes, &mut infos);
                    if let Some(t0) = t0 {
                        profile.infos_recompute += 1;
                        profile.t_infos += t0.elapsed();
                    }
                    infos_dirty = false;
                    infos_tw_dirty = false;
                }
                let t_sq0 = if profile_enabled {
                    Some(Instant::now())
                } else {
                    None
                };
                if profile_enabled {
                    profile.squeeze_attempts += 1;
                }
                run_stats.squeeze_attempts += 1;
                if let Some((ri, ep, ed)) =
                    squeeze_insertion(inst, &sol.routes, &infos, pickup, delivery)
                {
                    // Snapshot for rollback
                    bufs.squeeze_routes.clone_from(&sol.routes);
                    // Force-insert the pair
                    build_route_with_pair_buf(
                        &sol.routes[ri],
                        pickup,
                        ep,
                        delivery,
                        ed,
                        &mut route_buf,
                    );
                    std::mem::swap(&mut sol.routes[ri], &mut route_buf);
                    // Targeted repair: relocate pairs to reduce violations
                    let all_feasible = targeted_repair(
                        inst,
                        sol,
                        &mut infos,
                        &mut route_buf,
                        &mut bufs.info_buf,
                        rng,
                        REPAIR_MAX_ITERS,
                    );
                    if all_feasible {
                        // Squeeze + repair succeeded
                        unassigned_stack.pop();
                        consec_failures = 0;
                        infos_dirty = true;
                        for v in &mut route_versions {
                            *v = v.wrapping_add(1);
                        }
                        if profile_enabled {
                            profile.squeeze_success += 1;
                        }
                        run_stats.squeeze_successes += 1;
                    } else {
                        // Revert to pre-squeeze state
                        sol.routes.clone_from(&bufs.squeeze_routes);
                        infos_dirty = true;
                        for v in &mut route_versions {
                            *v = v.wrapping_add(1);
                        }
                    }
                }
                if let Some(t0) = t_sq0 {
                    profile.t_squeeze += t0.elapsed();
                }
            } else if profile_enabled {
                profile.ejection_fail += 1;
            }
            // If no ejection found: h_in stays on stack (equivalent to "eject h_in itself")

            // Check progress
            if unassigned_stack.len() < min_unassigned {
                min_unassigned = unassigned_stack.len();
                perturb_counter = 0;
                consec_failures = 0;
                // Save best state for potential revert
                best_state_routes.clone_from(&sol.routes);
                best_state_stack.clone_from(&unassigned_stack);
            }
            if unassigned_stack.is_empty() {
                if infeasible_mode && infeasible_insertions > 0 {
                    let repaired = targeted_repair(
                        inst,
                        sol,
                        &mut infos,
                        &mut route_buf,
                        &mut bufs.info_buf,
                        rng,
                        REPAIR_MAX_ITERS_FULL,
                    );
                    if repaired {
                        route_removed = true;
                        run_stats.infeasible_repair_successes += 1;
                    } else {
                        sol.routes.clone_from(&infeasible_snapshot_routes);
                        unassigned_stack.clone_from(&infeasible_snapshot_stack);
                        unassigned_stack.shuffle(rng);
                        infeasible_mode = false;
                        infeasible_insertions = 0;
                        infos_dirty = true;
                        for v in &mut route_versions {
                            *v = v.wrapping_add(1);
                        }
                        consec_failures = 0;
                    }
                } else {
                    route_removed = true;
                }
                if route_removed {
                    break;
                }
            }

            // Enter infeasible mode when deeply stuck
            if use_infeasible
                && !infeasible_mode
                && consec_failures >= INFEASIBLE_THRESHOLD
                && unassigned_stack.len() <= INFEASIBLE_MAX_UNASSIGNED
            {
                infeasible_mode = true;
                infeasible_insertions = 0;
                infeasible_pen = INFEASIBLE_PENALTY_INIT;
                infeasible_snapshot_routes.clone_from(&sol.routes);
                infeasible_snapshot_stack.clone_from(&unassigned_stack);
                run_stats.infeasible_entries += 1;
                // Upgrade to full infos for penalized insertion
                compute_all_infos(inst, &sol.routes, &mut infos);
                infos_dirty = false;
            }

            // Perturb (paper Algorithm 1, line 13: always after failed direct insertion)
            if infos_dirty {
                let t0 = if profile_enabled {
                    Some(Instant::now())
                } else {
                    None
                };
                if infeasible_mode {
                    compute_all_infos(inst, &sol.routes, &mut infos);
                } else {
                    compute_all_infos_light(inst, &sol.routes, &mut infos);
                }
                if let Some(t0) = t0 {
                    profile.infos_recompute += 1;
                    profile.t_infos += t0.elapsed();
                }
                infos_dirty = false;
            }
            // Adaptive perturbation: keep strong moves near completion, but
            // reduce work during long stalls to increase search throughput.
            let n_moves = adaptive_perturb_moves(unassigned_stack.len(), consec_failures);
            let t_pert0 = if profile_enabled {
                Some(Instant::now())
            } else {
                None
            };
            perturb_with_infos_bufs(
                inst,
                sol,
                n_moves,
                &mut infos,
                &mut bufs.info_buf,
                &mut bufs.perturb_buf1,
                &mut bufs.perturb_buf2,
                &mut bufs.perturb_route_buf,
                rng,
                pickup,
                &mut route_req_masks,
            );
            if infeasible_mode {
                // Infeasible mode needs full infos for penalized insertion
                compute_all_infos(inst, &sol.routes, &mut infos);
                infos_tw_dirty = false;
            } else {
                // Perturbation uses compute_light — tw/sector fields are stale.
                infos_tw_dirty = true;
            }
            // Masks maintained incrementally during perturbation if they were valid.
            // Only dirty if perturbation invalidated them (route count changed unexpectedly).
            route_masks_dirty = route_req_masks.len() != sol.routes.len();
            // Invalidate version cache: perturbation changed routes
            for v in &mut route_versions {
                *v = v.wrapping_add(1);
            }
            route_versions.resize(sol.routes.len(), 0);
            if let Some(t0) = t_pert0 {
                profile.perturb_calls += 1;
                profile.t_perturb += t0.elapsed();
            }
            run_stats.perturb_calls += 1;
            perturb_counter += 1;
            consec_failures += 1;

            // Shuffle-on-stagnation: shuffle the unassigned stack to try different
            // insertion order when stuck. Inspired by cuopt's shuffle interval.
            if consec_failures.is_multiple_of(SHUFFLE_INTERVAL) {
                unassigned_stack.shuffle(rng);
                if use_penalty_decay {
                    for p in &mut penalties {
                        if *p > 1 {
                            *p = ((*p * PENALTY_DECAY_NUM) / PENALTY_DECAY_DEN).max(1);
                        }
                    }
                }
            }

            // Revert-to-best: after long stall, revert to best-known state.
            // This prevents the search from drifting too far from productive regions.
            // The shuffled unassigned stack ensures a different insertion sequence.
            if reverts_remaining > 0
                && consec_failures > 0
                && consec_failures.is_multiple_of(REVERT_INTERVAL)
                && best_state_stack.len() < unassigned_stack.len()
            {
                sol.routes.clone_from(&best_state_routes);
                unassigned_stack.clone_from(&best_state_stack);
                unassigned_stack.shuffle(rng);
                infos_dirty = true;
                route_masks_dirty = true;
                // Invalidate version cache after revert
                for v in &mut route_versions {
                    *v = v.wrapping_add(1);
                }
                route_versions.resize(sol.routes.len(), 0);
                failed_inserts.resize_with(sol.routes.len(), || vec![u64::MAX; inst.n + 1]);
                reverts_remaining -= 1;
                if infeasible_mode {
                    infeasible_mode = false;
                    infeasible_insertions = 0;
                }
            }
            // Penalties keep long-term insertion difficulty memory, but optional
            // periodic decay prevents early failures from dominating too strongly.
        }

        if route_removed {
            any_removed = true;
            // Return control to main loop after each successful removal.
            // LS between removals restructures the solution, potentially enabling
            // the next removal. Avoids wasting time on a failed follow-up attempt.
            break;
        }
        // Failed to remove the route — exact rollback via snapshot (Step 1)
        sol.routes.clone_from(&bufs.rollback_routes);
        tried_routes.clone_from(&rollback_tried);
        tried_routes[route_idx] = true; // Mark this route as tried
        infos_dirty = true; // Routes restored — infos stale
        route_masks_dirty = true;
        // Reset penalties for next route attempt (penalties are route-specific)
        penalties.fill(1);
        // Continue to try a different route instead of giving up
    }

    if profile_enabled {
        profile.attempts = attempts;
        profile.print(
            fn_start.elapsed(),
            routes_before,
            sol.routes.len(),
            any_removed,
        );
    }
    if let Some(out) = stats_out {
        out.merge_from(&run_stats);
    }

    any_removed
}

/// Baseline number of ejected requests (old fixed behavior).
const K_BASE: usize = 3;
/// Absolute upper bound for dynamic ejection depth.
const K_HARD_MAX: usize = 6;

/// Best ejection candidate recorded during the lexicographic search.
struct EjBest {
    ri: usize,
    ep: usize,
    d_orig_pos: usize,
    ej_pickups: [usize; K_HARD_MAX],
    num_ej_pickups: usize,
    ej_nodes: [usize; 2 * K_HARD_MAX],
    num_ej_nodes: usize,
}

/// Constant context shared across all DFS calls for one ejection search.
/// Passed by reference to avoid copying 15+ parameters per recursive call.
struct EjCtx<'a> {
    inst: &'a Instance,
    route: &'a [usize],
    info: &'a RouteInfo,
    z: &'a [f64],
    smp: &'a [u64],
    penalties: &'a [u64],
    delivery: usize,
    demand_p: i32,
    demand_d: i32,
    cap: i32,
    ri: usize,
    ep: usize,
    n: usize,
    k_limit: usize,
    depth_limit: usize,
    best_p_sum: u64,
    best: Option<EjBest>,
    /// Precomputed arc_feas_t row for delivery (tests: is prev_nd→delivery feasible?)
    arc_feas_t_d: *const u64,
    /// Precomputed arc_feas row for delivery (tests: is delivery→node feasible?)
    arc_feas_d: *const u64,
    /// Total penalty of mandatory-eject pickups (conflicting with new pair).
    /// Any valid ejection set must include these. Used for lookahead pruning.
    mandatory_cost: u64,
    /// Pre-cached conflict row for pickup (avoids repeated pickup_to_req lookup in DFS).
    pickup_conflict_row: *const u64,
}

impl EjCtx<'_> {
    /// Fast conflict check using pre-cached conflict row for the target pickup.
    #[inline(always)]
    fn node_conflicts_with_pickup(&self, node: usize) -> bool {
        if self.pickup_conflict_row.is_null() {
            return false;
        }
        let rj = self.inst.pickup_to_req[node];
        if rj == usize::MAX {
            return false;
        }
        let word = unsafe { *self.pickup_conflict_row.add(rj >> 6) };
        (word >> (rj & 63)) & 1 != 0
    }
}

/// Count how many pickups in the route conflict with the given pickup.
#[inline]
fn count_route_conflicts(inst: &Instance, route: &[usize], pickup: usize) -> usize {
    let ri = inst.pickup_to_req[pickup];
    if ri == usize::MAX {
        return 0;
    }
    let conflict_row = &inst.conflict[ri * inst.conflict_stride..];
    let mut count = 0usize;
    for &v in route {
        if v == 0 {
            continue;
        }
        if inst.demand[v] <= 0 {
            continue;
        }
        let rj = unsafe { *inst.pickup_to_req.get_unchecked(v) };
        if rj == usize::MAX {
            continue;
        }
        let word = unsafe { *conflict_row.get_unchecked(rj >> 6) };
        if (word >> (rj & 63)) & 1 != 0 {
            count += 1;
            if count >= 3 {
                return count;
            }
        }
    }
    count
}

/// Lexicographic insertion-ejection search (Nagata & Kobayashi 2010, Section 3.3).
#[allow(clippy::needless_range_loop)]
fn find_best_ejection_dfs(
    inst: &Instance,
    sol: &Solution,
    pickup: usize,
    delivery: usize,
    penalties: &[u64],
    infos: &[RouteInfo],
    bufs: &mut EjSearchBufs,
    route_masks: &[Vec<u64>],
    k_limit: usize,
) -> Option<(usize, Vec<usize>, [usize; K_HARD_MAX], usize)> {
    let demand_p = inst.demand[pickup];
    let demand_d = inst.demand[delivery];
    let cap = inst.capacity;
    let mut best_p_sum = penalties[pickup];
    let mut best_result: Option<(usize, Vec<usize>, [usize; K_HARD_MAX], usize)> = None;

    // Precompute conflict counts per route (reused across Phase 1 and Phase 2).
    let use_masks = route_masks.len() == sol.routes.len();
    bufs.route_conflicts_buf.clear();
    if use_masks {
        bufs.route_conflicts_buf.extend(
            route_masks
                .iter()
                .map(|mask| inst.count_conflicts_with_mask(mask, pickup, 3) as u8),
        );
    } else {
        bufs.route_conflicts_buf.extend(
            sol.routes
                .iter()
                .map(|route| count_route_conflicts(inst, route, pickup) as u8),
        );
    }

    // Phase 1: Global penalty-ordered k=1 search with safe route conflict pre-filters.
    // We evaluate all eligible pickup candidates by increasing penalty and stop
    // once penalty exceeds current best_p_sum.
    struct Phase1Meta {
        ri: usize,
        node: usize,
        ep: usize,
        ed: usize,
        ins_cost: f64,
    }
    let mut phase1_best: Option<Phase1Meta> = None;
    bufs.phase1_cands_buf.clear();
    let arc_feas_t_p = inst.arc_feas_t_row(pickup);
    let arc_feas_t_d = inst.arc_feas_t_row(delivery);
    let arc_feas_d = inst.arc_feas_row(delivery);
    for ri in 0..sol.routes.len() {
        let route = &sol.routes[ri];
        let n = route.len();
        if n <= 2 {
            continue;
        }
        // Route-level conflict pre-filter: if 2+ conflicts, no k=1 ejection resolves all.
        let num_conflicts = bufs.route_conflicts_buf[ri] as usize;
        if num_conflicts >= 2 {
            continue;
        }
        // Route-level arc pre-filter: single pass for both arc_feas_t_p and arc_feas_d
        let mut has_t_p = false;
        let mut has_d = false;
        for &v in route {
            if !has_t_p && arc_bit(arc_feas_t_p, v) {
                has_t_p = true;
            }
            if !has_d && arc_bit(arc_feas_d, v) {
                has_d = true;
            }
            if has_t_p && has_d {
                break;
            }
        }
        if !has_t_p || !has_d {
            continue;
        }
        for &node in &route[1..n - 1] {
            if !inst.is_pickup(node) {
                continue;
            }
            // Pre-filter: candidates with penalty > best_p_sum will be skipped
            // by the sorted iteration's break condition anyway.
            if penalties[node] > best_p_sum {
                continue;
            }
            if num_conflicts == 1 && !inst.requests_conflict(pickup, node) {
                continue;
            }
            bufs.phase1_cands_buf.push((penalties[node], ri, node));
        }
    }
    bufs.phase1_cands_buf
        .sort_unstable_by_key(|&(pen, _, _)| pen);
    for &(pen, ri, node) in &bufs.phase1_cands_buf {
        if pen > best_p_sum {
            break;
        }
        let route = &sol.routes[ri];
        let d_node = inst.delivery_of(node);
        compute_info_without_pair(
            inst,
            route,
            node,
            d_node,
            &mut bufs.rwp_buf,
            &mut bufs.info_buf,
        );
        // Note: no conflict check needed here — Phase 1 pre-filtering guarantees
        // that ejecting `node` resolves all conflicts (0-conflict routes have none,
        // 1-conflict routes only consider the conflicting node as candidate).
        if let Some((ins_cost, ep, ed)) = first_feasible_insertion_with_info(
            inst,
            &bufs.rwp_buf,
            &bufs.info_buf,
            pickup,
            delivery,
        ) {
            if pen < best_p_sum {
                best_p_sum = pen;
                phase1_best = Some(Phase1Meta {
                    ri,
                    node,
                    ep,
                    ed,
                    ins_cost,
                });
                // Penalty=1 is minimum; k=2 always costs >=2, so return immediately.
                if best_p_sum <= 1 {
                    let new_route = build_route_with_pair(&bufs.rwp_buf, pickup, ep, delivery, ed);
                    let mut ejected = [0usize; K_HARD_MAX];
                    ejected[0] = node;
                    return Some((ri, new_route, ejected, 1));
                }
            } else if pen == best_p_sum {
                let better = match &phase1_best {
                    None => true,
                    Some(b) => ins_cost < b.ins_cost,
                };
                if better {
                    phase1_best = Some(Phase1Meta {
                        ri,
                        node,
                        ep,
                        ed,
                        ins_cost,
                    });
                }
            }
        }
    }

    // Build route for final Phase 1 best only (deferred construction)
    if let Some(b) = phase1_best {
        let d_node = inst.delivery_of(b.node);
        route_without_pair_buf(&sol.routes[b.ri], b.node, d_node, &mut bufs.rwp_buf);
        let new_route = build_route_with_pair(&bufs.rwp_buf, pickup, b.ep, delivery, b.ed);
        let mut ejected = [0usize; K_HARD_MAX];
        ejected[0] = b.node;
        best_result = Some((b.ri, new_route, ejected, 1));
    }

    // Phase 2: DFS for k>=2 ejections (Phase 1 exhaustively covers k=1).
    if best_p_sum <= 2 {
        return best_result;
    }
    let mut best_dfs: Option<EjBest> = None;

    // Pre-filter and sort routes by effective minimum cost for k>=2.
    bufs.ph2_routes_buf.clear();
    for ri in 0..sol.routes.len() {
        let route = &sol.routes[ri];
        let n = route.len();
        if n <= 2 {
            continue;
        }
        let num_conflicts = bufs.route_conflicts_buf[ri] as usize;
        if num_conflicts >= 3 {
            continue;
        }
        // Single pass: arc constraints + two-cheapest pickups + conflict cost.
        let mut has_t_p = false;
        let mut has_d = false;
        let mut p1 = u64::MAX;
        let mut p2 = u64::MAX;
        let mut conflict_cost = 0u64;
        let mut min_nonconflict = u64::MAX;
        for &v in route {
            if !has_t_p && arc_bit(arc_feas_t_p, v) {
                has_t_p = true;
            }
            if !has_d && arc_bit(arc_feas_d, v) {
                has_d = true;
            }
            if v != 0 && inst.is_pickup(v) {
                let p = penalties[v];
                if p < p1 {
                    p2 = p1;
                    p1 = p;
                } else if p < p2 {
                    p2 = p;
                }
                if num_conflicts > 0 && inst.requests_conflict(pickup, v) {
                    conflict_cost = conflict_cost.saturating_add(p);
                } else {
                    min_nonconflict = min_nonconflict.min(p);
                }
            }
        }
        if !has_t_p || !has_d {
            continue;
        }
        // k>=2 lower bounds: c=0 → p1+p2, c=1 → conflict+min_nonconflict, c=2 → conflict_cost
        let route_min_cost = match num_conflicts {
            0 => p1.saturating_add(p2),
            1 => conflict_cost.saturating_add(min_nonconflict),
            _ => conflict_cost,
        };
        if route_min_cost >= best_p_sum {
            continue;
        }
        bufs.ph2_routes_buf
            .push((ri, route_min_cost, conflict_cost));
    }
    bufs.ph2_routes_buf.sort_unstable_by_key(|&(_, mp, _)| mp);

    for &(ri, min_pen, mandatory_cost) in &bufs.ph2_routes_buf {
        if min_pen >= best_p_sum {
            break;
        }

        let route = &sol.routes[ri];
        let n = route.len();
        let info = &infos[ri];

        bufs.z_buf.resize(n, 0.0);
        bufs.z_buf[n - 1] = inst.late(0);
        for i in (0..n - 1).rev() {
            bufs.z_buf[i] = fmin(
                inst.late(route[i]),
                bufs.z_buf[i + 1] - inst.svc(route[i]) - info.edge_dists[i],
            );
        }

        bufs.smp_buf.resize(n, u64::MAX);
        for i in (0..n - 1).rev() {
            let node = route[i];
            bufs.smp_buf[i] = if node != 0 && inst.is_pickup(node) {
                std::cmp::min(penalties[node], bufs.smp_buf[i + 1])
            } else {
                bufs.smp_buf[i + 1]
            };
        }

        // Route-level smp pre-check
        if bufs.smp_buf[0] >= best_p_sum {
            continue;
        }

        // Find the earliest mandatory (conflicting) pickup in the route.
        // The DFS can only eject nodes at positions > ep, so ep must be
        // before all mandatory pickups. Otherwise they'd stay in the prefix,
        // creating a route with conflicting pairs.
        let pickup_req = inst.pickup_to_req[pickup];
        let pickup_conflict_row = if pickup_req != usize::MAX {
            unsafe {
                inst.conflict
                    .as_ptr()
                    .add(pickup_req * inst.conflict_stride)
            }
        } else {
            std::ptr::null()
        };

        // Find ep_limit: first mandatory (conflicting) pickup position.
        let mut ep_limit = n - 1;
        if mandatory_cost > 0 && !pickup_conflict_row.is_null() {
            for pos in 1..n - 1 {
                let v = route[pos];
                if inst.is_pickup(v) {
                    let rj = inst.pickup_to_req[v];
                    if rj != usize::MAX {
                        let word = unsafe { *pickup_conflict_row.add(rj >> 6) };
                        if (word >> (rj & 63)) & 1 != 0 {
                            ep_limit = pos;
                            break;
                        }
                    }
                }
            }
        }

        let mut ctx = EjCtx {
            inst,
            route,
            info,
            z: &bufs.z_buf,
            smp: &bufs.smp_buf,
            penalties,
            delivery,
            demand_p,
            demand_d,
            cap,
            ri,
            ep: 0,
            n,
            k_limit,
            depth_limit: 2 * k_limit + 1,
            best_p_sum,
            best: None,
            arc_feas_t_d,
            arc_feas_d,
            mandatory_cost,
            pickup_conflict_row,
        };

        for ep in 0..ep_limit {
            if info.load[ep] + demand_p > cap {
                continue;
            }
            if !arc_bit(arc_feas_t_p, route[ep]) {
                continue;
            }
            // Skip: no pickup after ep has penalty < best_p_sum
            if bufs.smp_buf[ep + 1] >= ctx.best_p_sum {
                continue;
            }
            let a_p = fmax(
                info.arrival[ep] + inst.svc(route[ep]) + inst.dist(route[ep], pickup),
                inst.early(pickup),
            );
            if a_p > inst.late(pickup) {
                continue;
            }

            ctx.ep = ep;
            ej_scan(
                &mut ctx,
                ep + 1,
                a_p,
                pickup,
                0,
                0,
                false,
                0,
                0,
                [0; K_HARD_MAX],
                0,
                [0; 2 * K_HARD_MAX],
                0,
                0,
                0, // paid_mandatory starts at 0
            );
        }

        if ctx.best_p_sum < best_p_sum {
            best_p_sum = ctx.best_p_sum;
            best_dfs = ctx.best;
        }

        if best_p_sum <= 2 {
            break;
        }
    }

    // Reconstruct from DFS result
    if let Some(b) = best_dfs {
        let route = &sol.routes[b.ri];
        let n = route.len();
        let mut new_route = Vec::with_capacity(n + 2);
        new_route.extend_from_slice(&route[..=b.ep]);
        new_route.push(pickup);
        for orig_pos in b.ep + 1..n {
            if orig_pos == b.d_orig_pos {
                new_route.push(delivery);
            }
            let node = route[orig_pos];
            let mut ejected = false;
            for k in 0..b.num_ej_nodes {
                if b.ej_nodes[k] == node {
                    ejected = true;
                    break;
                }
            }
            if !ejected {
                new_route.push(node);
            }
        }
        best_result = Some((b.ri, new_route, b.ej_pickups, b.num_ej_pickups));
    }

    best_result
}

/// Recursive lexicographic DFS for the ejection search.
/// Uses EjCtx for constant context to minimize stack frame size.
#[allow(clippy::too_many_arguments, clippy::needless_range_loop)]
fn ej_scan(
    ctx: &mut EjCtx,
    orig_start: usize,
    a_pred: f64,
    pred_node: usize,
    q_adj: i32,
    p_sum: u64,
    d_inserted: bool,
    depth: usize,
    num_ej_pickups: usize,
    ej_pickups: [usize; K_HARD_MAX],
    num_ej_nodes: usize,
    ej_nodes: [usize; 2 * K_HARD_MAX],
    d_orig_pos: usize,
    pending: usize,
    paid_mandatory: u64,
) {
    let inst = ctx.inst;
    let route = ctx.route;
    let n = ctx.n;
    let delivery = ctx.delivery;
    let demand_p = ctx.demand_p;
    let cap = ctx.cap;

    let mut a_prev = a_pred;
    let mut prev_nd = pred_node;

    // Early prune: penalty so far + remaining mandatory ejections >= best.
    let remaining_mandatory = ctx.mandatory_cost - paid_mandatory;
    if p_sum + remaining_mandatory >= ctx.best_p_sum {
        return;
    }

    for orig_pos in orig_start..n - 1 {
        let node = route[orig_pos];

        // Prune: if delivery not yet inserted and arrival already exceeds
        // delivery's latest time, no future position can accommodate it.
        if !d_inserted && a_prev > inst.late(delivery) {
            return;
        }

        // --- Option 1: Insert delivery before this position ---
        if !d_inserted
            && depth < ctx.depth_limit
            && arc_bit(ctx.arc_feas_t_d, prev_nd)
            && arc_bit(ctx.arc_feas_d, node)
        {
            let a_d = fmax(
                a_prev + inst.svc(prev_nd) + inst.dist(prev_nd, delivery),
                inst.early(delivery),
            );
            if a_d <= inst.late(delivery) {
                let q_adj_d = q_adj - ctx.demand_d;
                ej_scan(
                    ctx,
                    orig_pos,
                    a_d,
                    delivery,
                    q_adj_d,
                    p_sum,
                    true,
                    depth + 1,
                    num_ej_pickups,
                    ej_pickups,
                    num_ej_nodes,
                    ej_nodes,
                    orig_pos,
                    pending,
                    paid_mandatory,
                );
            }
        }

        // --- Option 2: Eject this node ---
        if depth < ctx.depth_limit && node != 0 {
            if inst.is_pickup(node) {
                if num_ej_pickups < ctx.k_limit {
                    let new_psum = p_sum + ctx.penalties[node];
                    // Mandatory-cost lookahead: account for remaining mandatory
                    // ejections when checking if this branch can improve.
                    let is_mandatory = ctx.node_conflicts_with_pickup(node);
                    let new_paid = if is_mandatory {
                        paid_mandatory + ctx.penalties[node]
                    } else {
                        paid_mandatory
                    };
                    let new_remaining = ctx.mandatory_cost - new_paid;
                    if new_psum + new_remaining < ctx.best_p_sum {
                        let mut ep2 = ej_pickups;
                        ep2[num_ej_pickups] = node;
                        let mut en2 = ej_nodes;
                        en2[num_ej_nodes] = node;
                        ej_scan(
                            ctx,
                            orig_pos + 1,
                            a_prev,
                            prev_nd,
                            q_adj + inst.demand[node],
                            new_psum,
                            d_inserted,
                            depth + 1,
                            num_ej_pickups + 1,
                            ep2,
                            num_ej_nodes + 1,
                            en2,
                            d_orig_pos,
                            pending + 1,
                            new_paid,
                        );
                    }
                }
            } else if inst.pair_pickup[node] != 0 {
                let its_p = inst.pickup_of(node);
                if (0..num_ej_pickups).any(|k| ej_pickups[k] == its_p) {
                    let mut en2 = ej_nodes;
                    en2[num_ej_nodes] = node;
                    ej_scan(
                        ctx,
                        orig_pos + 1,
                        a_prev,
                        prev_nd,
                        q_adj + inst.demand[node],
                        p_sum,
                        d_inserted,
                        depth + 1,
                        num_ej_pickups,
                        ej_pickups,
                        num_ej_nodes + 1,
                        en2,
                        d_orig_pos,
                        pending - 1,
                        paid_mandatory,
                    );
                }
            }
        }

        // --- Option 3: Keep this node ---
        // Condition (v): delivery whose pickup was ejected cannot be kept
        if node != 0 && inst.pair_pickup[node] != 0 {
            let its_p = inst.pickup_of(node);
            if (0..num_ej_pickups).any(|k| ej_pickups[k] == its_p) {
                return;
            }
        }

        // Conflict prune: a pickup that conflicts with the new pair can never be
        // kept — the two requests are mutually exclusive in any route.
        if node != 0 && inst.is_pickup(node) && ctx.node_conflicts_with_pickup(node) {
            return;
        }

        // Arc feasibility: prev_nd → node
        if !inst.is_arc_feasible(prev_nd, node) {
            return;
        }

        // Compute arrival at this position
        let a_pos = fmax(
            a_prev + inst.svc(prev_nd) + inst.dist(prev_nd, node),
            inst.early(node),
        );

        if a_pos > inst.late(node) {
            return;
        }

        // Capacity check at this kept position
        let load_here = ctx.info.load[orig_pos] + demand_p - q_adj;
        if orig_pos < n - 1 && load_here > cap {
            return;
        }
        // Strong suffix capacity lower bound: once delivery is inserted and pickup-ejection
        // budget is exhausted, remaining load shift is fixed; if suffix max load still exceeds
        // capacity, this branch can never become feasible.
        if d_inserted
            && num_ej_pickups >= ctx.k_limit
            && ctx.info.suffix_max_load[orig_pos] + demand_p - q_adj > cap
        {
            return;
        }

        // Terminal check: delivery inserted, all ejected pairs complete
        if d_inserted && pending == 0 {
            if a_pos <= ctx.z[orig_pos] {
                let cap_ok = orig_pos >= n - 1
                    || ctx.info.suffix_max_load[orig_pos] + demand_p - q_adj <= cap;
                if cap_ok {
                    if p_sum < ctx.best_p_sum {
                        ctx.best_p_sum = p_sum;
                        ctx.best = Some(EjBest {
                            ri: ctx.ri,
                            ep: ctx.ep,
                            d_orig_pos,
                            ej_pickups,
                            num_ej_pickups,
                            ej_nodes,
                            num_ej_nodes,
                        });
                    }
                    return;
                }
            }
            if num_ej_pickups >= ctx.k_limit {
                return;
            }
            // Suffix min penalty pruning with mandatory cost
            if orig_pos + 1 < n {
                if p_sum
                    .saturating_add(remaining_mandatory)
                    .saturating_add(ctx.smp[orig_pos + 1])
                    >= ctx.best_p_sum
                {
                    return;
                }
            } else {
                return;
            }
        }

        a_prev = a_pos;
        prev_nd = node;
    }

    // --- Post-loop: insert delivery before final depot ---
    if !d_inserted
        && depth < ctx.depth_limit
        && inst.is_arc_feasible(prev_nd, delivery)
        && inst.is_arc_feasible(delivery, 0)
    {
        let a_d = fmax(
            a_prev + inst.svc(prev_nd) + inst.dist(prev_nd, delivery),
            inst.early(delivery),
        );
        if a_d <= inst.late(delivery) {
            let a_depot = fmax(
                a_d + inst.svc(delivery) + inst.dist(delivery, 0),
                inst.early(0),
            );
            if a_depot <= inst.late(0) && pending == 0 && p_sum < ctx.best_p_sum {
                ctx.best_p_sum = p_sum;
                ctx.best = Some(EjBest {
                    ri: ctx.ri,
                    ep: ctx.ep,
                    d_orig_pos: n - 1,
                    ej_pickups,
                    num_ej_pickups,
                    ej_nodes,
                    num_ej_nodes,
                });
            }
        }
    }

    // --- Post-loop: delivery already inserted, check depot feasibility ---
    if d_inserted && pending == 0 && inst.is_arc_feasible(prev_nd, 0) {
        let a_depot = fmax(
            a_prev + inst.svc(prev_nd) + inst.dist(prev_nd, 0),
            inst.early(0),
        );
        if a_depot <= inst.late(0) && p_sum < ctx.best_p_sum {
            ctx.best_p_sum = p_sum;
            ctx.best = Some(EjBest {
                ri: ctx.ri,
                ep: ctx.ep,
                d_orig_pos,
                ej_pickups,
                num_ej_pickups,
                ej_nodes,
                num_ej_nodes,
            });
        }
    }
}

/// Brute-force ejection search (original algorithm).
/// Searches routes × requests for k=1..2 ejections (k_limit is accepted for API parity).
#[allow(clippy::too_many_arguments)]
fn find_best_ejection_bf(
    inst: &Instance,
    sol: &Solution,
    pickup: usize,
    delivery: usize,
    penalties: &[u64],
    bufs: &mut EjSearchBufs,
    _k_limit: usize,
) -> Option<(usize, Vec<usize>, [usize; K_HARD_MAX], usize)> {
    let mut best_p_sum = penalties[pickup];

    // Deferred construction: store indices, build route only for final best.
    let mut best_ri: usize = 0;
    let mut best_ej1: usize = 0;
    let mut best_ej2: usize = 0; // 0 = k=1 ejection (pickup IDs are always >= 1)
    let mut best_ep: usize = 0;
    let mut best_ed: usize = 0;
    let mut best_empty = false; // true if ejection empties the route
    let mut found = false;

    let nr = sol.routes.len();

    // Pre-compute arc feasibility rows outside the per-route loop
    let arc_feas_t_p = inst.arc_feas_t_row(pickup);
    let arc_feas_d = inst.arc_feas_row(delivery);

    for ri in 0..nr {
        let route = &sol.routes[ri];

        // Route-level arc feasibility pre-filter
        if !route.iter().any(|&v| arc_bit(arc_feas_t_p, v)) {
            continue;
        }
        if !route.iter().any(|&v| arc_bit(arc_feas_d, v)) {
            continue;
        }

        let mut pickups_in_route: Vec<usize> = route
            .iter()
            .copied()
            .filter(|&v| v != 0 && inst.is_pickup(v))
            .collect();
        pickups_in_route.sort_unstable_by_key(|&p| penalties[p]);

        // k=1: try ejecting each single request
        for &req_p in &pickups_in_route {
            let p_sum = penalties[req_p];
            if p_sum >= best_p_sum {
                break; // sorted ascending — all remaining have >= penalty
            }
            let req_d = inst.delivery_of(req_p);
            compute_info_without_pair(
                inst,
                route,
                req_p,
                req_d,
                &mut bufs.rwp_buf,
                &mut bufs.info_buf,
            );

            // Conflict check: skip if inserting pair conflicts with remaining route
            if bufs.rwp_buf.len() >= 2 && inst.has_route_conflict(&bufs.rwp_buf, pickup) {
                continue;
            }

            if bufs.rwp_buf.len() < 2 {
                if is_single_pair_feasible(inst, pickup, delivery) {
                    best_p_sum = p_sum;
                    best_ri = ri;
                    best_ej1 = req_p;
                    best_ej2 = 0;
                    best_empty = true;
                    found = true;
                }
                continue;
            }

            if let Some((_, ep, ed)) = first_feasible_insertion_with_info(
                inst,
                &bufs.rwp_buf,
                &bufs.info_buf,
                pickup,
                delivery,
            ) {
                best_p_sum = p_sum;
                best_ri = ri;
                best_ej1 = req_p;
                best_ej2 = 0;
                best_ep = ep;
                best_ed = ed;
                best_empty = false;
                found = true;
            }
        }

        // k=2: try ejecting each pair of requests from this route
        for i in 0..pickups_in_route.len() {
            let req_p1 = pickups_in_route[i];
            // Early termination: if penalties[req_p1] alone >= best_p_sum,
            // then any pair including req_p1 (or later) exceeds threshold.
            if penalties[req_p1] >= best_p_sum {
                break;
            }
            for &req_p2 in &pickups_in_route[i + 1..] {
                let p_sum = penalties[req_p1] + penalties[req_p2];
                if p_sum >= best_p_sum {
                    break; // sorted — remaining j have higher penalty
                }
                let req_d1 = inst.delivery_of(req_p1);
                let req_d2 = inst.delivery_of(req_p2);
                route_without_pair_buf(route, req_p1, req_d1, &mut bufs.rwp_buf);
                route_without_pair_buf(&bufs.rwp_buf, req_p2, req_d2, &mut bufs.rwp_buf2);

                // Conflict check: skip if inserting pair conflicts with remaining route
                if bufs.rwp_buf2.len() >= 2 && inst.has_route_conflict(&bufs.rwp_buf2, pickup) {
                    continue;
                }

                if bufs.rwp_buf2.len() < 2 {
                    if is_single_pair_feasible(inst, pickup, delivery) {
                        best_p_sum = p_sum;
                        best_ri = ri;
                        best_ej1 = req_p1;
                        best_ej2 = req_p2;
                        best_empty = true;
                        found = true;
                    }
                    continue;
                }

                bufs.info_buf.compute(inst, &bufs.rwp_buf2);
                if let Some((_, ep, ed)) = first_feasible_insertion_with_info(
                    inst,
                    &bufs.rwp_buf2,
                    &bufs.info_buf,
                    pickup,
                    delivery,
                ) {
                    best_p_sum = p_sum;
                    best_ri = ri;
                    best_ej1 = req_p1;
                    best_ej2 = req_p2;
                    best_ep = ep;
                    best_ed = ed;
                    best_empty = false;
                    found = true;
                }
            }
        }
    }

    if !found {
        return None;
    }

    // Build the result route only for the final best ejection.
    let (ejected, num_ejected) = if best_ej2 == 0 {
        let mut ep = [0usize; K_HARD_MAX];
        ep[0] = best_ej1;
        (ep, 1)
    } else {
        let mut ep = [0usize; K_HARD_MAX];
        ep[0] = best_ej1;
        ep[1] = best_ej2;
        (ep, 2)
    };

    if best_empty {
        return Some((best_ri, vec![0, pickup, delivery, 0], ejected, num_ejected));
    }

    // Rebuild route-without-pair(s), then insert pickup/delivery at stored positions.
    let route = &sol.routes[best_ri];
    let ej1_d = inst.delivery_of(best_ej1);
    if best_ej2 == 0 {
        route_without_pair_buf(route, best_ej1, ej1_d, &mut bufs.rwp_buf);
        let new_route = build_route_with_pair(&bufs.rwp_buf, pickup, best_ep, delivery, best_ed);
        Some((best_ri, new_route, ejected, num_ejected))
    } else {
        let ej2_d = inst.delivery_of(best_ej2);
        route_without_pair_buf(route, best_ej1, ej1_d, &mut bufs.rwp_buf);
        route_without_pair_buf(&bufs.rwp_buf, best_ej2, ej2_d, &mut bufs.rwp_buf2);
        let new_route = build_route_with_pair(&bufs.rwp_buf2, pickup, best_ep, delivery, best_ed);
        Some((best_ri, new_route, ejected, num_ejected))
    }
}

/// Limited cross-route ejection chaining: immediately try to place an ejected pair
/// into another route before sending it back to unassigned stack.
fn try_cross_route_insert(
    inst: &Instance,
    sol: &mut Solution,
    pickup: usize,
    infos: &mut [RouteInfo],
    route_buf: &mut Vec<usize>,
    route_req_masks: &mut [Vec<u64>],
    use_conflict_masks: bool,
    route_versions: &mut [u64],
    rng: &mut impl Rng,
) -> bool {
    if sol.routes.is_empty() {
        return false;
    }
    let delivery = inst.delivery_of(pickup);
    let nr = sol.routes.len();
    let start = fast_range(rng, nr);
    for off in 0..nr {
        let ri = fast_mod(start + off, nr);
        let has_conflict = if use_conflict_masks && route_req_masks.len() == nr {
            inst.has_conflict_with_mask(&route_req_masks[ri], pickup)
        } else {
            inst.has_route_conflict(&sol.routes[ri], pickup)
        };
        if has_conflict {
            continue;
        }
        if let Some((_, ep, ed)) =
            first_feasible_insertion_with_info(inst, &sol.routes[ri], &infos[ri], pickup, delivery)
        {
            compute_info_with_pair(
                inst,
                &sol.routes[ri],
                pickup,
                ep,
                delivery,
                ed,
                route_buf,
                &mut infos[ri],
            );
            std::mem::swap(&mut sol.routes[ri], route_buf);
            if use_conflict_masks && route_req_masks.len() == nr {
                let req = inst.pickup_to_req[pickup];
                if req != usize::MAX {
                    route_req_masks[ri][req >> 6] |= 1u64 << (req & 63);
                }
            }
            if ri < route_versions.len() {
                route_versions[ri] = route_versions[ri].wrapping_add(1);
            }
            return true;
        }
    }
    false
}

/// Squeeze insertion: find the route + position that minimizes total excess
/// (TW violation + capacity violation), ignoring feasibility. Uses very high
/// penalty coefficients so the search prioritizes low-violation positions.
/// Returns (route_index, ep, ed) or None if no route can host the pair.
fn squeeze_insertion(
    inst: &Instance,
    routes: &[Vec<usize>],
    infos: &[RouteInfo],
    pickup: usize,
    delivery: usize,
) -> Option<(usize, usize, usize)> {
    // High penalty: 1000× distance scale to heavily penalize violations.
    // This makes the search prefer positions with near-zero excess.
    let squeeze_pen = 10_000.0;
    let mut best_cost = f64::MAX;
    let mut best_ri = 0;
    let mut best_ep = 0;
    let mut best_ed = 0;
    let mut found = false;

    for ri in 0..routes.len() {
        if inst.has_route_conflict(&routes[ri], pickup) {
            continue;
        }
        if let Some((cost, ep, ed)) = best_pair_insertion_penalized(
            inst,
            &routes[ri],
            &infos[ri],
            pickup,
            delivery,
            squeeze_pen,
            squeeze_pen,
        ) && cost < best_cost
        {
            best_cost = cost;
            best_ri = ri;
            best_ep = ep;
            best_ed = ed;
            found = true;
        }
    }

    if found {
        Some((best_ri, best_ep, best_ed))
    } else {
        None
    }
}

/// Pre-GES evacuation: move easy pairs out of the victim route so GES only
/// handles the hard-to-place residual pairs.  Returns the number of pairs evacuated.
pub fn starve_route(
    inst: &Instance,
    sol: &mut Solution,
    rng: &mut impl Rng,
    deadline: Instant,
    sc_duals: Option<&[f64]>,
) -> usize {
    let use_sc = sc_duals.is_some() && fast_coin(rng);
    let target = ges_starvation_target();
    let nr = sol.routes.len();
    if nr < 2 {
        return 0;
    }

    // --- 1. Compute infos & select victim via SpatialTemporalTop3 scoring (deterministic) ---
    let mut infos: Vec<RouteInfo> = Vec::with_capacity(nr);
    for _ in 0..nr {
        infos.push(RouteInfo::new());
    }
    compute_all_infos_light(inst, &sol.routes, &mut infos);

    // Pass 1: compute raw scores for each route
    struct RouteScore {
        idx: usize,
        overlap: f64,
        idle: f64,
        short: f64,
        sc_raw: f64,
    }
    let mut route_scores: Vec<RouteScore> = Vec::with_capacity(nr);

    for i in 0..nr {
        let ri = &sol.routes[i];
        let ii = &infos[i];
        let idle = route_idle_ratio(inst, ri, ii);
        let short = 1.0 / (ri.len().saturating_sub(2).max(1) as f64);
        let a0 = if ri.len() > 2 { ii.arrival[1] } else { 0.0 };
        let a1 = if ri.len() > 2 {
            ii.arrival[ri.len() - 2]
        } else {
            0.0
        };

        let scale_i = ii.distance / (ri.len().saturating_sub(1).max(1) as f64);
        let mut pair_scores: Vec<f64> = Vec::new();
        for (rj, ij) in infos.iter().enumerate() {
            if rj == i {
                continue;
            }
            let rj_route = &sol.routes[rj];
            let b0 = if rj_route.len() > 2 {
                ij.arrival[1]
            } else {
                0.0
            };
            let b1 = if rj_route.len() > 2 {
                ij.arrival[rj_route.len() - 2]
            } else {
                0.0
            };
            let temporal = temporal_overlap_ratio(a0, a1, b0, b1);
            let anchor_dist = route_anchor_min_dist(inst, ri, rj_route);
            let scale_j = ij.distance / (rj_route.len().saturating_sub(1).max(1) as f64);
            let spatial = 1.0 / (1.0 + anchor_dist / (1.0 + 0.5 * (scale_i + scale_j)));
            pair_scores.push(0.7 * spatial + 0.3 * temporal);
        }
        let overlap = if pair_scores.is_empty() {
            0.0
        } else {
            pair_scores.sort_unstable_by(|a, b| b.total_cmp(a));
            let k = pair_scores.len().min(3);
            let mut s = 0.0;
            for &v in pair_scores.iter().take(k) {
                s += v;
            }
            s / k as f64
        };
        let sc_raw = if let Some(duals) = sc_duals {
            let mut dual_sum = 0.0;
            let mut np = 0usize;
            for &v in ri {
                if v != 0 && inst.is_pickup(v) {
                    dual_sum += duals[v];
                    np += 1;
                }
            }
            if np > 0 { dual_sum / np as f64 } else { 0.0 }
        } else {
            0.0
        };
        route_scores.push(RouteScore {
            idx: i,
            overlap,
            idle,
            short,
            sc_raw,
        });
    }

    // Pass 2: normalize SC and pick best victim
    let mut best_victim = 0usize;
    let mut best_score = f64::NEG_INFINITY;

    if use_sc {
        let max_sc = route_scores
            .iter()
            .map(|e| e.sc_raw)
            .fold(0.0_f64, f64::max);
        for rs in &route_scores {
            let sc_replaceable = if max_sc > 1e-12 {
                1.0 - (rs.sc_raw / max_sc).clamp(0.0, 1.0)
            } else {
                0.0
            };
            let score = 0.45 * rs.overlap + 0.2 * rs.idle + 0.1 * rs.short + 0.25 * sc_replaceable;
            if score > best_score {
                best_score = score;
                best_victim = rs.idx;
            }
        }
    } else {
        for rs in &route_scores {
            let score = 0.6 * rs.overlap + 0.3 * rs.idle + 0.1 * rs.short;
            if score > best_score {
                best_score = score;
                best_victim = rs.idx;
            }
        }
    }

    // Count PD pairs in victim route
    let num_pairs = sol.routes[best_victim]
        .iter()
        .filter(|&&v| v != 0 && inst.is_pickup(v))
        .count();
    if num_pairs <= target {
        return 0;
    }

    // --- 2. Compute evacuability for each pair in victim ---
    struct PairEvac {
        pickup: usize,
        delivery: usize,
        accepting_routes: usize,
        best_cost: f64,
    }

    let mut pair_evacs: Vec<PairEvac> = Vec::new();
    let victim_route = &sol.routes[best_victim];
    for &v in victim_route {
        if v == 0 || !inst.is_pickup(v) {
            continue;
        }
        let p = v;
        let d = inst.delivery_of(p);
        let mut accepting = 0usize;
        let mut best_cost = f64::MAX;
        for (ri, info_ri) in infos.iter().enumerate() {
            if ri == best_victim {
                continue;
            }
            if let Some((cost, _, _)) =
                best_pair_insertion_with_info::<false>(inst, &sol.routes[ri], info_ri, p, d)
            {
                accepting += 1;
                if cost < best_cost {
                    best_cost = cost;
                }
            }
        }
        pair_evacs.push(PairEvac {
            pickup: p,
            delivery: d,
            accepting_routes: accepting,
            best_cost,
        });
    }

    // Sort by evacuability descending: most accepting routes first, then lowest cost
    pair_evacs.sort_unstable_by(|a, b| {
        b.accepting_routes
            .cmp(&a.accepting_routes)
            .then_with(|| a.best_cost.total_cmp(&b.best_cost))
    });

    // --- 3. Greedy evacuation loop ---
    let mut evacuated = 0usize;
    let mut virt_route_buf: Vec<usize> = Vec::new();
    let mut dest_route_buf: Vec<usize> = Vec::new();
    let mut info_buf = RouteInfo::new();

    for pe in &pair_evacs {
        if Instant::now() >= deadline {
            break;
        }

        // Count current PD pairs in victim
        let current_pairs = sol.routes[best_victim]
            .iter()
            .filter(|&&v| v != 0 && inst.is_pickup(v))
            .count();
        if current_pairs <= target {
            break;
        }

        if pe.accepting_routes == 0 {
            continue; // leave for GES
        }

        let p = pe.pickup;
        let d = pe.delivery;

        // Verify pair is still in victim (defensive)
        if !sol.routes[best_victim].contains(&p) {
            continue;
        }

        // Find best insertion across all other routes (recompute with current state)
        let mut best_ri = usize::MAX;
        let mut best_ep = 0usize;
        let mut best_ed = 0usize;
        let mut best_insert_cost = f64::MAX;

        for (ri, info_ri) in infos.iter_mut().enumerate() {
            if ri == best_victim {
                continue;
            }
            // Recompute info for target route (may have changed from earlier insertions)
            info_ri.compute_light(inst, &sol.routes[ri]);
            if let Some((cost, ep, ed)) =
                best_pair_insertion_with_info::<false>(inst, &sol.routes[ri], info_ri, p, d)
                && cost < best_insert_cost
            {
                best_insert_cost = cost;
                best_ri = ri;
                best_ep = ep;
                best_ed = ed;
            }
        }

        if best_ri == usize::MAX {
            continue; // no feasible insertion found
        }

        // Remove pair from victim: compute_info_without_pair
        compute_info_without_pair(
            inst,
            &sol.routes[best_victim],
            p,
            d,
            &mut virt_route_buf,
            &mut info_buf,
        );
        std::mem::swap(&mut sol.routes[best_victim], &mut virt_route_buf);
        std::mem::swap(&mut infos[best_victim], &mut info_buf);

        // Insert pair into target: compute_info_with_pair
        compute_info_with_pair(
            inst,
            &sol.routes[best_ri],
            p,
            best_ep,
            d,
            best_ed,
            &mut dest_route_buf,
            &mut info_buf,
        );
        std::mem::swap(&mut sol.routes[best_ri], &mut dest_route_buf);
        std::mem::swap(&mut infos[best_ri], &mut info_buf);

        evacuated += 1;
    }

    // --- 4. Cleanup: remove empty victim route if fully evacuated ---
    let victim_pairs = sol.routes[best_victim]
        .iter()
        .filter(|&&v| v != 0 && inst.is_pickup(v))
        .count();
    if victim_pairs == 0 && sol.routes[best_victim].len() <= 2 {
        sol.routes.swap_remove(best_victim);
        // infos vector doesn't need cleanup — it's local
    }

    evacuated
}

/// Targeted repair: relocate pairs from infeasible routes to feasible positions,
/// reducing total violation at each step. Returns true if all routes become feasible.
fn targeted_repair(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    route_buf: &mut Vec<usize>,
    info_buf: &mut RouteInfo,
    rng: &mut impl Rng,
    max_iters: usize,
) -> bool {
    compute_all_infos(inst, &sol.routes, infos);

    let mut order_buf: Vec<usize> = Vec::new();
    for _iter in 0..max_iters {
        // Collect infeasible route indices
        let mut infeasible: Vec<usize> = Vec::new();
        for (ri, info_ri) in infos.iter().enumerate() {
            if info_ri.tw_excess > 1e-10 || info_ri.cap_excess > 1e-10 {
                infeasible.push(ri);
            }
        }
        if infeasible.is_empty() {
            return true;
        }

        // Pick random infeasible route
        let src_ri = infeasible[fast_range(rng, infeasible.len())];
        let pickup = sample_pickup(inst, &sol.routes[src_ri], rng);
        if pickup == 0 {
            continue;
        }
        let delivery = inst.delivery_of(pickup);

        let src_violation_before = infos[src_ri].tw_excess + infos[src_ri].cap_excess;

        // Build random route order for destination search
        order_buf.clear();
        for ri in 0..sol.routes.len() {
            if ri != src_ri {
                order_buf.push(ri);
            }
        }
        order_buf.shuffle(rng);

        let mut moved = false;

        // Try feasible relocation first
        for &dst_ri in &order_buf {
            if inst.has_route_conflict(&sol.routes[dst_ri], pickup) {
                continue;
            }
            if let Some((_, ep, ed)) = first_feasible_insertion_with_info(
                inst,
                &sol.routes[dst_ri],
                &infos[dst_ri],
                pickup,
                delivery,
            ) {
                // Compute src violation after removal
                compute_info_without_pair(
                    inst,
                    &sol.routes[src_ri],
                    pickup,
                    delivery,
                    route_buf,
                    info_buf,
                );
                let src_violation_after = info_buf.tw_excess + info_buf.cap_excess;
                if src_violation_after < src_violation_before - 1e-10 {
                    // Apply: insert into dst
                    build_route_with_pair_buf(
                        &sol.routes[dst_ri],
                        pickup,
                        ep,
                        delivery,
                        ed,
                        route_buf,
                    );
                    std::mem::swap(&mut sol.routes[dst_ri], route_buf);
                    infos[dst_ri].compute(inst, &sol.routes[dst_ri]);
                    // Remove from src
                    route_without_pair_buf(&sol.routes[src_ri], pickup, delivery, route_buf);
                    std::mem::swap(&mut sol.routes[src_ri], route_buf);
                    if sol.routes[src_ri].len() <= 2 {
                        sol.routes.swap_remove(src_ri);
                        infos.swap_remove(src_ri);
                    } else {
                        infos[src_ri].compute(inst, &sol.routes[src_ri]);
                    }
                    moved = true;
                    break;
                }
            }
        }
        if moved {
            continue;
        }

        // Fallback: penalized relocation (accept if total violation decreases)
        let mut best_pen_cost = f64::MAX;
        let mut best_dst = 0usize;
        let mut best_ep = 0usize;
        let mut best_ed = 0usize;
        let mut found_pen = false;
        for &dst_ri in &order_buf {
            if inst.has_route_conflict(&sol.routes[dst_ri], pickup) {
                continue;
            }
            if let Some((cost, ep, ed)) = best_pair_insertion_penalized(
                inst,
                &sol.routes[dst_ri],
                &infos[dst_ri],
                pickup,
                delivery,
                10_000.0,
                10_000.0,
            ) && cost < best_pen_cost
            {
                best_pen_cost = cost;
                best_dst = dst_ri;
                best_ep = ep;
                best_ed = ed;
                found_pen = true;
            }
        }
        if found_pen {
            // Compute dst violation after insertion
            compute_info_with_pair(
                inst,
                &sol.routes[best_dst],
                pickup,
                best_ep,
                delivery,
                best_ed,
                route_buf,
                info_buf,
            );
            let dst_violation_after = info_buf.tw_excess + info_buf.cap_excess;
            // Compute src violation after removal
            compute_info_without_pair(
                inst,
                &sol.routes[src_ri],
                pickup,
                delivery,
                &mut order_buf,
                info_buf,
            );
            let src_violation_after = info_buf.tw_excess + info_buf.cap_excess;
            let dst_violation_before = infos[best_dst].tw_excess + infos[best_dst].cap_excess;
            let total_before = src_violation_before + dst_violation_before;
            let total_after = src_violation_after + dst_violation_after;
            if total_after < total_before - 1e-10 {
                // Apply move: insert into dst
                build_route_with_pair_buf(
                    &sol.routes[best_dst],
                    pickup,
                    best_ep,
                    delivery,
                    best_ed,
                    route_buf,
                );
                std::mem::swap(&mut sol.routes[best_dst], route_buf);
                infos[best_dst].compute(inst, &sol.routes[best_dst]);
                // Remove from src
                route_without_pair_buf(&sol.routes[src_ri], pickup, delivery, route_buf);
                std::mem::swap(&mut sol.routes[src_ri], route_buf);
                if sol.routes[src_ri].len() <= 2 {
                    sol.routes.swap_remove(src_ri);
                    infos.swap_remove(src_ri);
                } else {
                    infos[src_ri].compute(inst, &sol.routes[src_ri]);
                }
            }
        }
    }

    // Final check
    sol.routes.iter().all(|r| is_route_feasible(inst, r))
}

/// Perturbation: randomly apply `PairMove` or `SwapMove` operations.
pub fn perturb_solution(
    inst: &Instance,
    sol: &mut Solution,
    num_moves: usize,
    info_buf: &mut RouteInfo,
    rng: &mut impl Rng,
) {
    for _ in 0..num_moves {
        if sol.routes.len() < 2 {
            break;
        }
        if fast_coin(rng) {
            pair_move(inst, sol, info_buf, rng);
        } else {
            swap_move(inst, sol, info_buf, rng);
        }
    }
}

/// Perturbation with cached infos: maintains infos incrementally instead of
/// invalidating the entire cache. Saves ~50 RouteInfo::compute calls per cycle.
#[allow(dead_code)]
pub fn perturb_with_infos(
    inst: &Instance,
    sol: &mut Solution,
    num_moves: usize,
    infos: &mut Vec<RouteInfo>,
    info_buf: &mut RouteInfo,
    rng: &mut impl Rng,
) {
    let mut buf1: Vec<usize> = Vec::new();
    let mut buf2: Vec<usize> = Vec::new();
    let mut route_buf: Vec<usize> = Vec::new();
    let mut masks: Vec<Vec<u64>> = Vec::new();
    perturb_with_infos_bufs(
        inst,
        sol,
        num_moves,
        infos,
        info_buf,
        &mut buf1,
        &mut buf2,
        &mut route_buf,
        rng,
        0,
        &mut masks,
    );
}

/// Inner perturbation with externally-provided buffers (avoids per-call allocation).
/// `target_pickup`: if nonzero, occasionally do conflict-resolution moves for this pair.
/// `masks`: if non-empty, use bitmask conflict checks and maintain incrementally.
fn perturb_with_infos_bufs(
    inst: &Instance,
    sol: &mut Solution,
    num_moves: usize,
    infos: &mut Vec<RouteInfo>,
    info_buf: &mut RouteInfo,
    buf1: &mut Vec<usize>,
    buf2: &mut Vec<usize>,
    route_buf: &mut Vec<usize>,
    rng: &mut impl Rng,
    target_pickup: usize,
    masks: &mut Vec<Vec<u64>>,
) {
    let use_masks = masks.len() == sol.routes.len();
    let use_temporal = ges_use_temporal_perturb();
    for i in 0..num_moves {
        if sol.routes.len() < 2 {
            break;
        }
        // Every 8th move: try to resolve a conflict for the target pair.
        if target_pickup != 0 && i.trailing_zeros() >= 3 {
            conflict_move_cached(inst, sol, infos, masks, route_buf, rng, target_pickup);
            continue;
        }
        // Every 8th move (offset 3, avoiding conflict at offset 0): intra-route
        // temporal perturbation. Relocates a PD pair within the same route to
        // change arrival times, creating slack for TW-constrained insertions.
        if use_temporal && i & 7 == 3 {
            intra_route_reorder(inst, sol, infos, buf1, route_buf, rng);
            continue;
        }
        if fast_coin(rng) {
            pair_move_cached(inst, sol, infos, masks, route_buf, rng);
        } else {
            swap_move_cached(
                inst, sol, infos, masks, info_buf, buf1, buf2, route_buf, rng,
            );
        }
    }
    // If masks were not initially available, don't try to maintain them
    if use_masks && masks.len() != sol.routes.len() {
        // Route count changed during perturbation; caller should rebuild
        masks.clear();
    }
}

/// `PairMove`: randomly move a PD pair from one route to another.
fn pair_move(inst: &Instance, sol: &mut Solution, info_buf: &mut RouteInfo, rng: &mut impl Rng) {
    let nr = sol.routes.len();
    if nr < 2 {
        return;
    }

    let r1_idx = rng.random_range(0..nr);
    let pickup = sample_pickup(inst, &sol.routes[r1_idx], rng);
    if pickup == 0 {
        return;
    }
    let delivery = inst.delivery_of(pickup);

    let r2_idx = loop {
        let r = rng.random_range(0..nr);
        if r != r1_idx {
            break r;
        }
    };

    info_buf.compute(inst, &sol.routes[r2_idx]);
    // Biased Perturbation (Sartori & Buriol 2020, Section 4.5.2):
    // Pick best feasible insertion instead of first, biasing toward better positions.
    if let Some((_, ep, ed)) = best_pair_insertion_with_info::<false>(
        inst,
        &sol.routes[r2_idx],
        info_buf,
        pickup,
        delivery,
    ) {
        let new_r1 = route_without_pair(&sol.routes[r1_idx], pickup, delivery);
        let new_r2 = build_route_with_pair(&sol.routes[r2_idx], pickup, ep, delivery, ed);

        sol.routes[r2_idx] = new_r2;
        if new_r1.len() <= 2 {
            sol.routes.swap_remove(r1_idx);
        } else {
            sol.routes[r1_idx] = new_r1;
        }
    }
}

/// Conflict-resolution move: find a pair conflicting with `target_pickup` in a random
/// route and try to relocate it to another route. This clears conflicts to enable
/// future direct insertion of the target pair.
fn conflict_move_cached(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    masks: &mut Vec<Vec<u64>>,
    route_buf: &mut Vec<usize>,
    rng: &mut impl Rng,
    target_pickup: usize,
) {
    let nr = sol.routes.len();
    if nr < 2 {
        return;
    }
    let use_masks = masks.len() == nr;
    // Pick a random starting route, scan for a conflicting pair.
    // Use mask pre-filter when available: O(stride) per route vs O(n) node scan.
    let start_ri = fast_range(rng, nr);
    let mut conflict_p = 0usize;
    let mut source_ri = 0usize;
    for off in 0..nr {
        let ri = fast_mod(start_ri + off, nr);
        // Fast mask pre-filter: skip routes with no conflict at all
        if use_masks && !inst.has_conflict_with_mask(&masks[ri], target_pickup) {
            continue;
        }
        for &v in &sol.routes[ri] {
            if v != 0 && inst.is_pickup(v) && inst.requests_conflict(target_pickup, v) {
                conflict_p = v;
                break;
            }
        }
        if conflict_p != 0 {
            source_ri = ri;
            break;
        }
    }
    if conflict_p == 0 {
        return;
    }
    let conflict_d = inst.delivery_of(conflict_p);
    // Try to move the conflicting pair to another route (random order).
    let start_r2 = fast_range(rng, nr);
    for off in 0..nr {
        let r2 = fast_mod(start_r2 + off, nr);
        if r2 == source_ri {
            continue;
        }
        let has_conflict = if use_masks {
            inst.has_conflict_with_mask(&masks[r2], conflict_p)
        } else {
            inst.has_route_conflict(&sol.routes[r2], conflict_p)
        };
        if has_conflict {
            continue;
        }
        if let Some((_, ep, ed)) = first_feasible_insertion_with_info(
            inst,
            &sol.routes[r2],
            &infos[r2],
            conflict_p,
            conflict_d,
        ) {
            // Fused: build route with pair + compute info in one pass
            compute_info_with_pair(
                inst,
                &sol.routes[r2],
                conflict_p,
                ep,
                conflict_d,
                ed,
                route_buf,
                &mut infos[r2],
            );
            std::mem::swap(&mut sol.routes[r2], route_buf);
            if use_masks {
                let req = inst.pickup_to_req[conflict_p];
                if req != usize::MAX {
                    masks[r2][req >> 6] |= 1u64 << (req & 63);
                }
            }
            // Fused: remove pair from source route + recompute info in one pass
            compute_info_without_pair(
                inst,
                &sol.routes[source_ri],
                conflict_p,
                conflict_d,
                route_buf,
                &mut infos[source_ri],
            );
            if route_buf.len() <= 2 {
                sol.routes.swap_remove(source_ri);
                infos.swap_remove(source_ri);
                if use_masks && source_ri < masks.len() {
                    masks.swap_remove(source_ri);
                }
            } else {
                std::mem::swap(&mut sol.routes[source_ri], route_buf);
                if use_masks {
                    let req = inst.pickup_to_req[conflict_p];
                    if req != usize::MAX {
                        masks[source_ri][req >> 6] &= !(1u64 << (req & 63));
                    }
                }
            }
            return;
        }
    }
}

/// `PairMove` with cached infos: uses infos for insertion check, updates after move.
/// When `masks` is non-empty, uses O(stride) bitmask conflict check and maintains masks.
fn pair_move_cached(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    masks: &mut Vec<Vec<u64>>,
    route_buf: &mut Vec<usize>,
    rng: &mut impl Rng,
) {
    let nr = sol.routes.len();
    if nr < 2 {
        return;
    }
    let use_masks = masks.len() == nr;

    let r1_idx = fast_range(rng, nr);
    let pickup = sample_pickup(inst, &sol.routes[r1_idx], rng);
    if pickup == 0 {
        return;
    }
    let delivery = inst.delivery_of(pickup);

    let r2_idx = fast_mod(r1_idx + 1 + fast_range(rng, nr - 1), nr);

    // Conflict pre-filter: skip if pair conflicts with any pair in target route.
    let has_conflict = if use_masks {
        inst.has_conflict_with_mask(&masks[r2_idx], pickup)
    } else {
        inst.has_route_conflict(&sol.routes[r2_idx], pickup)
    };
    if has_conflict {
        return;
    }
    // Arc pre-filters: skip if no node in target route can precede pickup
    // or if delivery cannot be followed by any node in target route.
    {
        let r2 = &sol.routes[r2_idx];
        let arc_t_p = inst.arc_feas_t_row(pickup);
        let arc_d = inst.arc_feas_row(delivery);
        let mut has_t_p = false;
        let mut has_d = false;
        for &v in r2 {
            if !has_t_p && arc_bit(arc_t_p, v) {
                has_t_p = true;
            }
            if !has_d && arc_bit(arc_d, v) {
                has_d = true;
            }
            if has_t_p && has_d {
                break;
            }
        }
        if !has_t_p || !has_d {
            return;
        }
    }

    // Biased Perturbation (Sartori & Buriol 2020, Section 4.5.2):
    // Pick best feasible insertion among sparse-arc candidates. SPARSE=true reduces
    // O(n²) insertion cost by focusing on positions with good arc structure, while
    // changing the perturbation trajectory to explore different solution regions.
    if let Some((_, ep, ed)) = best_pair_insertion_with_info::<true>(
        inst,
        &sol.routes[r2_idx],
        &infos[r2_idx],
        pickup,
        delivery,
    ) {
        // Fused: build route with pair + compute info in one pass
        // (saves one O(n) pass vs. separate build + compute_light)
        compute_info_with_pair(
            inst,
            &sol.routes[r2_idx],
            pickup,
            ep,
            delivery,
            ed,
            route_buf,
            &mut infos[r2_idx],
        );
        std::mem::swap(&mut sol.routes[r2_idx], route_buf);
        if use_masks {
            // Incremental: add pickup's request bit to dest mask
            let req = inst.pickup_to_req[pickup];
            if req != usize::MAX {
                masks[r2_idx][req >> 6] |= 1u64 << (req & 63);
            }
        }

        // Fused: remove pair from source route + recompute info in one pass
        compute_info_without_pair(
            inst,
            &sol.routes[r1_idx],
            pickup,
            delivery,
            route_buf,
            &mut infos[r1_idx],
        );
        if route_buf.len() <= 2 {
            sol.routes.swap_remove(r1_idx);
            infos.swap_remove(r1_idx);
            if use_masks && r1_idx < masks.len() {
                masks.swap_remove(r1_idx);
            }
        } else {
            std::mem::swap(&mut sol.routes[r1_idx], route_buf);
            if use_masks {
                // Incremental: clear pickup's request bit from source mask
                let req = inst.pickup_to_req[pickup];
                if req != usize::MAX {
                    masks[r1_idx][req >> 6] &= !(1u64 << (req & 63));
                }
            }
        }
    }
}

/// `SwapMove`: swap PD pairs between two routes at best positions.
#[allow(clippy::similar_names)]
fn swap_move(inst: &Instance, sol: &mut Solution, info_buf: &mut RouteInfo, rng: &mut impl Rng) {
    let nr = sol.routes.len();
    if nr < 2 {
        return;
    }

    let r1_idx = rng.random_range(0..nr);
    let r2_idx = loop {
        let r = rng.random_range(0..nr);
        if r != r1_idx {
            break r;
        }
    };

    let p1 = sample_pickup(inst, &sol.routes[r1_idx], rng);
    let p2 = sample_pickup(inst, &sol.routes[r2_idx], rng);
    if p1 == 0 || p2 == 0 {
        return;
    }

    let d1 = inst.delivery_of(p1);
    let d2 = inst.delivery_of(p2);

    let r1_without = route_without_pair(&sol.routes[r1_idx], p1, d1);
    let r2_without = route_without_pair(&sol.routes[r2_idx], p2, d2);

    if r1_without.len() < 2 || r2_without.len() < 2 {
        return;
    }

    info_buf.compute(inst, &r2_without);
    let swap_ins1 = first_feasible_insertion_with_info(inst, &r2_without, info_buf, p1, d1);
    let Some((_, ep1, ed1)) = swap_ins1 else {
        return;
    };
    info_buf.compute(inst, &r1_without);
    if let Some((_, ep2, ed2)) =
        first_feasible_insertion_with_info(inst, &r1_without, info_buf, p2, d2)
    {
        sol.routes[r1_idx] = build_route_with_pair(&r1_without, p2, ep2, d2, ed2);
        sol.routes[r2_idx] = build_route_with_pair(&r2_without, p1, ep1, d1, ed1);
    }
}

/// `SwapMove` with cached infos: uses pre-allocated buffers and first_feasible for speed.
/// When `masks` is non-empty, uses O(stride) bitmask conflict check and maintains masks.
#[allow(clippy::similar_names)]
fn swap_move_cached(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut [RouteInfo],
    masks: &mut [Vec<u64>],
    info_buf: &mut RouteInfo,
    buf1: &mut Vec<usize>,
    buf2: &mut Vec<usize>,
    route_buf: &mut Vec<usize>,
    rng: &mut impl Rng,
) {
    let nr = sol.routes.len();
    if nr < 2 {
        return;
    }
    let use_masks = masks.len() == nr;

    let r1_idx = fast_range(rng, nr);
    let r2_idx = fast_mod(r1_idx + 1 + fast_range(rng, nr - 1), nr);

    let p1 = sample_pickup(inst, &sol.routes[r1_idx], rng);
    let p2 = sample_pickup(inst, &sol.routes[r2_idx], rng);
    if p1 == 0 || p2 == 0 {
        return;
    }

    let d1 = inst.delivery_of(p1);
    let d2 = inst.delivery_of(p2);

    // Check BOTH directions' conflict pre-filters before expensive work.
    // For mask path: compute both checks without short-circuit, branch once.
    // Reduces branch misprediction (single 50/50 branch → one ~75% branch).
    if use_masks {
        let c1 = inst.has_conflict_with_mask(&masks[r2_idx], p1);
        let c2 = inst.has_conflict_with_mask(&masks[r1_idx], p2);
        if c1 | c2 {
            return;
        }
    } else {
        if inst.has_route_conflict(&sol.routes[r2_idx], p1) {
            return;
        }
        if inst.has_route_conflict(&sol.routes[r1_idx], p2) {
            return;
        }
    }
    {
        let arc_t_p1 = inst.arc_feas_t_row(p1);
        let arc_d1 = inst.arc_feas_row(d1);
        let arc_t_p2 = inst.arc_feas_t_row(p2);
        let arc_d2 = inst.arc_feas_row(d2);
        // p1→r2: need predecessor for p1 and successor for d1 in r2
        let (mut ok_t_p1, mut ok_d1) = (false, false);
        for &v in &sol.routes[r2_idx] {
            if !ok_t_p1 && arc_bit(arc_t_p1, v) {
                ok_t_p1 = true;
            }
            if !ok_d1 && arc_bit(arc_d1, v) {
                ok_d1 = true;
            }
            if ok_t_p1 && ok_d1 {
                break;
            }
        }
        if !ok_t_p1 || !ok_d1 {
            return;
        }
        // p2→r1: need predecessor for p2 and successor for d2 in r1
        let (mut ok_t_p2, mut ok_d2) = (false, false);
        for &v in &sol.routes[r1_idx] {
            if !ok_t_p2 && arc_bit(arc_t_p2, v) {
                ok_t_p2 = true;
            }
            if !ok_d2 && arc_bit(arc_d2, v) {
                ok_d2 = true;
            }
            if ok_t_p2 && ok_d2 {
                break;
            }
        }
        if !ok_t_p2 || !ok_d2 {
            return;
        }
    }

    // Fused: build buf2 (r2 without p2/d2) + compute its RouteInfo in one pass
    compute_info_without_pair(inst, &sol.routes[r2_idx], p2, d2, buf2, info_buf);
    if buf2.len() < 2 {
        return;
    }
    let swap_ins1 = first_feasible_insertion_with_info(inst, buf2, info_buf, p1, d1);
    let Some((_, ep1, ed1)) = swap_ins1 else {
        return;
    };

    // Fused: build buf1 (r1 without p1/d1) + compute its RouteInfo in one pass
    compute_info_without_pair(inst, &sol.routes[r1_idx], p1, d1, buf1, info_buf);
    if buf1.len() < 2 {
        return;
    }
    if let Some((_, ep2, ed2)) = first_feasible_insertion_with_info(inst, buf1, info_buf, p2, d2) {
        // Fused: build route with pair + compute info in one pass (2× saves)
        compute_info_with_pair(inst, buf1, p2, ep2, d2, ed2, route_buf, &mut infos[r1_idx]);
        std::mem::swap(&mut sol.routes[r1_idx], route_buf);
        compute_info_with_pair(inst, buf2, p1, ep1, d1, ed1, route_buf, &mut infos[r2_idx]);
        std::mem::swap(&mut sol.routes[r2_idx], route_buf);
        if use_masks {
            // Incremental: p1 moved r1→r2, p2 moved r2→r1
            let req1 = inst.pickup_to_req[p1];
            let req2 = inst.pickup_to_req[p2];
            if req1 != usize::MAX {
                masks[r1_idx][req1 >> 6] &= !(1u64 << (req1 & 63));
                masks[r2_idx][req1 >> 6] |= 1u64 << (req1 & 63);
            }
            if req2 != usize::MAX {
                masks[r2_idx][req2 >> 6] &= !(1u64 << (req2 & 63));
                masks[r1_idx][req2 >> 6] |= 1u64 << (req2 & 63);
            }
        }
    }
}

/// Intra-route temporal perturbation: relocate a PD pair to a random feasible
/// position within the SAME route, changing arrival times to create "slack".
/// Does not change route membership (masks unchanged), only node ordering.
#[inline(never)]
fn intra_route_reorder(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut [RouteInfo],
    reduced_buf: &mut Vec<usize>,
    new_route_buf: &mut Vec<usize>,
    rng: &mut impl Rng,
) {
    let nr = sol.routes.len();
    if nr == 0 {
        return;
    }

    let ri = fast_range(rng, nr);
    let route = &sol.routes[ri];
    // Need at least 2 PD pairs (route = [0, p1, d1, p2, d2, 0] = len 6)
    if route.len() < 6 {
        return;
    }

    let pickup = sample_pickup(inst, route, rng);
    if pickup == 0 {
        return;
    }
    let delivery = inst.delivery_of(pickup);

    // Remove pair from route into reduced_buf
    route_without_pair_buf(route, pickup, delivery, reduced_buf);
    let reduced_len = reduced_buf.len();
    if reduced_len < 4 {
        return;
    }

    let num_edges = reduced_len - 1;

    // Try random positions (up to 8 attempts)
    for _ in 0..8 {
        let ep = fast_range(rng, num_edges);
        let remaining = num_edges - ep;
        let ed = ep + fast_range(rng, remaining);

        // Build candidate route
        build_route_with_pair_buf(reduced_buf, pickup, ep, delivery, ed, new_route_buf);

        // Check feasibility: TW + capacity (precedence guaranteed by ep <= ed)
        if is_tw_cap_feasible(inst, new_route_buf) {
            // Accept: recompute info and commit
            infos[ri].compute(inst, new_route_buf);
            std::mem::swap(&mut sol.routes[ri], new_route_buf);
            return;
        }
    }
}

/// Fast TW + capacity feasibility check for a complete route.
/// Precedence is NOT checked — use only when precedence is guaranteed by construction
/// (e.g., build_route_with_pair where pickup edge <= delivery edge).
#[inline(never)]
fn is_tw_cap_feasible(inst: &Instance, route: &[usize]) -> bool {
    let n = route.len();
    let mut load: i32 = 0;
    let mut time = inst.early(0);
    for i in 1..n {
        let prev = route[i - 1];
        let curr = route[i];
        time = fmax(
            time + inst.svc(prev) + inst.dist(prev, curr),
            inst.early(curr),
        );
        if time > inst.late(curr) {
            return false;
        }
        if i < n - 1 {
            load += inst.demand[curr];
            if load > inst.capacity || load < 0 {
                return false;
            }
        }
    }
    true
}