pdp_lns 0.1.1

use crate::instance::Instance;
use crate::local_search;
use crate::set_covering::RoutePool;
use crate::solution::{
    PenaltyState, RouteInfo, Solution, best_pair_insertion_with_info, build_route_with_pair,
    compute_all_infos, compute_info_with_pair, single_pair_feasible_distance,
};
use rand::prelude::*;
use rand::rngs::SmallRng;
use std::sync::OnceLock;
use std::time::{Duration, Instant};

/// LAHC array length — shorter list for faster adaptation to improving solutions.
const LHC_LEN: usize = 500;

/// Persistent LAHC state that survives across `lns()` calls.
/// When kept alive between iterations, the costs array accumulates
/// history, enabling genuine acceptance of worse solutions.
pub struct LahcState {
    pub costs: Vec<f64>,
    pub iter: usize,
    #[allow(dead_code)]
    pub current_cost: f64,
}

impl LahcState {
    pub fn new(initial_cost: f64) -> Self {
        Self {
            costs: vec![initial_cost; LHC_LEN],
            iter: 0,
            current_cost: initial_cost,
        }
    }

    /// Reset the LAHC state (called when sol is reset to best).
    pub fn reset(&mut self) {
        self.costs.clear();
    }
}

#[derive(Debug, Clone, Copy)]
#[repr(u8)]
pub enum DestroyOp {
    RouteElimination = 0,
    Shaw = 1,
    WorstDistance = 2,
    Cluster = 3,
    StringSisrs = 4,
    RcRemoval = 5,
    Random = 6,
    HistoryRemoval = 7,
    LateArrival = 8,
    DetourMisplacement = 9,
    ChainedRuin = 10,
}

impl DestroyOp {
    const COUNT: usize = 11;
    const ALL: [Self; 11] = [
        Self::RouteElimination,
        Self::Shaw,
        Self::WorstDistance,
        Self::Cluster,
        Self::StringSisrs,
        Self::RcRemoval,
        Self::Random,
        Self::HistoryRemoval,
        Self::LateArrival,
        Self::DetourMisplacement,
        Self::ChainedRuin,
    ];
}

#[derive(Debug, Clone, Copy)]
#[repr(u8)]
pub enum RepairOp {
    Regret2 = 0,
    Regret3 = 1,
    Regret45 = 2,
    RegretMax = 3,
    InsertionDifficulty = 4,
}

impl RepairOp {
    const COUNT: usize = 5;
    const ALL: [Self; 5] = [
        Self::Regret2,
        Self::Regret3,
        Self::Regret45,
        Self::RegretMax,
        Self::InsertionDifficulty,
    ];
    fn to_k(self, num_routes: usize, rng: &mut impl Rng) -> usize {
        match self {
            Self::Regret2 => 2,
            Self::Regret3 => 3,
            Self::Regret45 => {
                if rng.random_bool(0.5) {
                    4
                } else {
                    5
                }
            }
            Self::RegretMax => num_routes.max(2),
            // ID only needs best and 2nd-best, but k=2 ensures top-2 tracking
            Self::InsertionDifficulty => 2,
        }
    }
    fn use_insertion_difficulty(self) -> bool {
        matches!(self, Self::InsertionDifficulty)
    }
}

/// ALNS operator weights and scoring state (Ropke & Pisinger 2006).
pub struct AlnsState {
    destroy_weights: [f64; DestroyOp::COUNT],
    repair_weights: [f64; RepairOp::COUNT],
    destroy_scores: [f64; DestroyOp::COUNT],
    destroy_counts: [u32; DestroyOp::COUNT],
    repair_scores: [f64; RepairOp::COUNT],
    repair_counts: [u32; RepairOp::COUNT],
    segment_iter: usize,
}

/// Reaction factor: how much new segment performance affects weights.
const ALNS_REACTION_FACTOR: f64 = 0.1;
/// Number of iterations per ALNS segment before weight update.
const ALNS_SEGMENT_LENGTH: usize = 100;
/// Minimum operator weight (prevents full disabling).
const ALNS_MIN_WEIGHT: f64 = 0.01;

impl Default for AlnsState {
    fn default() -> Self {
        Self::new()
    }
}

impl AlnsState {
    pub fn new() -> Self {
        // Initial weights:
        // RouteElim=14%, Shaw=25%, Worst=12%, Cluster=7%, String=5%,
        // RC=7%, Random=9%, HistoryRemoval=2%, LateArrival=4%, DetourMisplacement=12%, ChainedRuin=5% (rounding: -1 Shaw, -1 Random, -1 Detour → +3 ChainedRuin ≈ 5%)
        let dw = [
            0.14, // RouteElimination
            0.25, // Shaw
            0.12, // WorstDistance
            0.07, // Cluster
            0.05, // StringSisrs
            0.07, // RcRemoval
            0.09, // Random
            0.02, // HistoryRemoval
            0.04, // LateArrival
            0.12, // DetourMisplacement
            0.05, // ChainedRuin
        ];
        Self {
            destroy_weights: dw,
            repair_weights: [1.0, 1.0, 1.0, 1.0, 1.0],
            destroy_scores: [0.0; DestroyOp::COUNT],
            destroy_counts: [0; DestroyOp::COUNT],
            repair_scores: [0.0; RepairOp::COUNT],
            repair_counts: [0; RepairOp::COUNT],
            segment_iter: 0,
        }
    }

    pub fn reset(&mut self) {
        *self = Self::new();
    }

    pub fn select_destroy(&self, rng: &mut impl Rng) -> DestroyOp {
        let total: f64 = self.destroy_weights.iter().sum();
        let mut r = rng.random::<f64>() * total;
        for (i, &w) in self.destroy_weights.iter().enumerate() {
            r -= w;
            if r <= 0.0 {
                return DestroyOp::ALL[i];
            }
        }
        DestroyOp::ALL[DestroyOp::COUNT - 1]
    }

    pub fn select_repair(&self, rng: &mut impl Rng) -> RepairOp {
        let total: f64 = self.repair_weights.iter().sum();
        let mut r = rng.random::<f64>() * total;
        for (i, &w) in self.repair_weights.iter().enumerate() {
            r -= w;
            if r <= 0.0 {
                return RepairOp::ALL[i];
            }
        }
        RepairOp::ALL[RepairOp::COUNT - 1]
    }

    fn record_outcome(&mut self, destroy: DestroyOp, repair: RepairOp, score: f64) {
        let di = destroy as usize;
        let ri = repair as usize;
        self.destroy_scores[di] += score;
        self.destroy_counts[di] += 1;
        self.repair_scores[ri] += score;
        self.repair_counts[ri] += 1;
        self.segment_iter += 1;
        if self.segment_iter >= ALNS_SEGMENT_LENGTH {
            self.update_weights();
        }
    }

    fn update_weights(&mut self) {
        let r = ALNS_REACTION_FACTOR;
        for i in 0..DestroyOp::COUNT {
            if self.destroy_counts[i] > 0 {
                let avg = self.destroy_scores[i] / f64::from(self.destroy_counts[i]);
                self.destroy_weights[i] = (1.0 - r) * self.destroy_weights[i] + r * avg;
            }
            self.destroy_weights[i] = self.destroy_weights[i].max(ALNS_MIN_WEIGHT);
        }
        for i in 0..RepairOp::COUNT {
            if self.repair_counts[i] > 0 {
                let avg = self.repair_scores[i] / f64::from(self.repair_counts[i]);
                self.repair_weights[i] = (1.0 - r) * self.repair_weights[i] + r * avg;
            }
            self.repair_weights[i] = self.repair_weights[i].max(ALNS_MIN_WEIGHT);
        }
        self.destroy_scores = [0.0; DestroyOp::COUNT];
        self.destroy_counts = [0; DestroyOp::COUNT];
        self.repair_scores = [0.0; RepairOp::COUNT];
        self.repair_counts = [0; RepairOp::COUNT];
        self.segment_iter = 0;
    }
}

/// Record-to-Record Travel acceptance state (AHGSLNS paper).
pub struct RrtState {
    pub record: f64,
    pub deviation: f64,
    d0: f64,
}

impl RrtState {
    pub fn new(initial_cost: f64) -> Self {
        let d0 = 0.02 * initial_cost;
        Self {
            record: initial_cost,
            deviation: d0,
            d0,
        }
    }

    /// Decay deviation and replenish if too small.
    pub fn decay(&mut self) {
        self.deviation *= 0.9985;
        if self.deviation < 0.02 * self.d0 {
            self.deviation = 0.25 * self.d0;
        }
    }

    /// Reset to sentinel (reinit on next lns() call).
    pub fn reset(&mut self) {
        self.deviation = 0.0;
        self.d0 = 0.0;
        self.record = f64::MAX;
    }

    pub fn init(&mut self, cost: f64) {
        self.d0 = 0.02 * cost;
        self.deviation = self.d0;
        self.record = cost;
    }
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[repr(u8)]
pub enum AcceptanceCriterion {
    Lahc = 0,
    Rrt = 1,
}

impl AcceptanceCriterion {
    const COUNT: usize = 2;
    const ALL: [Self; 2] = [Self::Lahc, Self::Rrt];
}

/// Adaptive acceptance criterion weights.
#[allow(dead_code)]
const ACCEPT_SEGMENT_LENGTH: usize = 6;
#[allow(dead_code)]
const ACCEPT_REACTION_FACTOR: f64 = 0.15;
#[allow(dead_code)]
const ACCEPT_MIN_WEIGHT: f64 = 0.05;
#[allow(dead_code)]
const ACCEPT_SCORE_BEST: f64 = 33.0;
/// Base score when criterion is selected but no best improvement.
/// Set to 1.0 (not 5.0) to avoid inflating weights of criteria that
/// accept worse solutions (RRT) without finding improvements.
#[allow(dead_code)]
const ACCEPT_SCORE_BASE: f64 = 1.0;

/// Wraps LAHC and RRT with adaptive roulette-wheel selection.
pub struct AcceptanceState {
    pub lahc: LahcState,
    pub rrt: RrtState,
    pub active: AcceptanceCriterion,
    pub current_cost: f64,
    weights: [f64; AcceptanceCriterion::COUNT],
    #[allow(dead_code)]
    scores: [f64; AcceptanceCriterion::COUNT],
    #[allow(dead_code)]
    counts: [u32; AcceptanceCriterion::COUNT],
    #[allow(dead_code)]
    segment_iter: usize,
}

impl AcceptanceState {
    pub fn new(initial_cost: f64) -> Self {
        Self {
            lahc: LahcState::new(initial_cost),
            rrt: RrtState::new(initial_cost),
            active: AcceptanceCriterion::Lahc,
            current_cost: initial_cost,
            // Fixed weights: LAHC=0.80, RRT=0.20.
            // LAHC is the most reliable for distance optimization (compares
            // absolute values). RRT inflates its weight via the adaptive
            // selector because LNS-thread improvements get misattributed.
            weights: [0.80, 0.20],
            scores: [0.0; AcceptanceCriterion::COUNT],
            counts: [0; AcceptanceCriterion::COUNT],
            segment_iter: 0,
        }
    }

    /// Roulette-wheel selection of acceptance criterion.
    pub fn select_criterion(&mut self, rng: &mut impl Rng) {
        let total: f64 = self.weights.iter().sum();
        let mut r = rng.random::<f64>() * total;
        for (i, &w) in self.weights.iter().enumerate() {
            r -= w;
            if r <= 0.0 {
                self.active = AcceptanceCriterion::ALL[i];
                return;
            }
        }
        self.active = AcceptanceCriterion::ALL[AcceptanceCriterion::COUNT - 1];
    }

    /// Record outcome. Weights are fixed (not adaptive) to prevent SA/RRT
    /// from inflating via misattributed LNS-thread improvements.
    #[allow(clippy::unused_self)]
    pub fn record_outcome(&mut self, _best_improved: bool) {
        // No-op: weights stay fixed at [0.80, 0.20].
    }

    #[allow(dead_code)]
    fn update_weights(&mut self) {
        let r = ACCEPT_REACTION_FACTOR;
        for i in 0..AcceptanceCriterion::COUNT {
            if self.counts[i] > 0 {
                let avg = self.scores[i] / f64::from(self.counts[i]);
                self.weights[i] = (1.0 - r) * self.weights[i] + r * avg;
            }
            self.weights[i] = self.weights[i].max(ACCEPT_MIN_WEIGHT);
        }
        self.scores = [0.0; AcceptanceCriterion::COUNT];
        self.counts = [0; AcceptanceCriterion::COUNT];
        self.segment_iter = 0;
    }

    /// Reset criteria states (called when sol is reset to best).
    /// Preserves adaptive weights across resets.
    pub fn reset(&mut self) {
        self.lahc.reset();
        self.rrt.reset();
    }

    /// Initialize or reinitialize all criteria for a given cost.
    pub fn init(&mut self, cost: f64) {
        self.lahc = LahcState::new(cost);
        self.rrt.init(cost);
        self.current_cost = cost;
    }

    pub fn weights_str(&self) -> String {
        format!("LAHC={:.2} RRT={:.2}", self.weights[0], self.weights[1])
    }
}

/// Arc frequency tracker for history removal (Masson et al. 2012, §3.2).
/// Counts how often each directed arc j→j' appears in best solutions.
pub struct ArcHistory {
    freq: Vec<u32>,
    stride: usize,
    num_solutions: u32,
}

impl ArcHistory {
    pub fn new(n: usize) -> Self {
        let stride = n + 1;
        Self {
            freq: vec![0u32; stride * stride],
            stride,
            num_solutions: 0,
        }
    }

    /// Add all arcs from a solution to the frequency counts.
    pub fn add_solution(&mut self, sol: &Solution) {
        for route in &sol.routes {
            for w in route.windows(2) {
                self.freq[w[0] * self.stride + w[1]] += 1;
            }
        }
        self.num_solutions += 1;
    }

    /// Halve all frequencies to prevent old solutions from dominating.
    pub fn decay(&mut self) {
        for f in &mut self.freq {
            *f >>= 1;
        }
        self.num_solutions = self.num_solutions.div_ceil(2);
    }

    #[inline]
    pub fn get(&self, from: usize, to: usize) -> u32 {
        self.freq[from * self.stride + to]
    }

    pub fn is_empty(&self) -> bool {
        self.num_solutions == 0
    }

    pub fn num_solutions(&self) -> u32 {
        self.num_solutions
    }
}

/// Min PD pairs removed by removal heuristic (Table 1).
const MIN_REMOVE: usize = 9;
/// Max PD pairs removed by removal heuristic (Table 1).
const MAX_REMOVE: usize = 92;
/// Probability of selecting route elimination removal.
const ROUTE_ELIM_PROB: f64 = 0.2;
/// Probability of selecting Shaw removal.
const SHAW_PROB: f64 = 0.40;
/// Probability of selecting worst-distance removal.
const WORST_PROB: f64 = 0.16;
/// Probability of selecting cluster removal (geographic proximity of PD pairs).
const CLUSTER_PROB: f64 = 0.0;
/// Probability of selecting SISRS string removal (Christiaens & Vanden Berghe 2020).
const STRING_PROB: f64 = 0.0;
/// Probability of selecting LP reduced-cost guided removal (used by tabu fixed-prob selection).
#[allow(dead_code)]
const RC_REMOVAL_PROB: f64 = 0.08;
/// Min consecutive iterations without improvement before stopping.
const MIN_NO_IMPROVE: usize = 266;
/// Alt stopping: iterations without improvement.
const ALT_NO_IMPROVE: usize = 402;
/// Maximum string cardinality for SISRS removal.
const SISRS_L_MAX: usize = 10;
/// Split rate for SISRS split string procedure.
const SISRS_ALPHA: f64 = 0.5;
/// Split depth for SISRS split string procedure.
const SISRS_BETA: f64 = 0.01;
/// Tabu search: max iterations in one stagnation kick.
const TABU_MAX_ITERS: usize = 160;
/// Random neighbors sampled per iteration.
const TABU_RANDOM_NEIGHBORS: usize = 8;
/// Previously promising moves to re-evaluate each iteration.
const TABU_CACHE_REEVAL: usize = 3;
/// Max number of cached removal moves.
const TABU_CACHE_SIZE: usize = 24;
/// Tabu tenure in iterations.
const TABU_TENURE: usize = 12;
/// Min tabu tenure for non-improving moves.
const TABU_TENURE_MIN: usize = 4;
/// Objective big-M used in LNS candidate comparison.
const LNS_OBJ_BIG_M: f64 = 1e9;
/// Probability of trying branch-and-bound repair instead of regret insertion.
const BNB_REPAIR_PROB: f64 = 0.065_344_595_364_038_75;
/// Apply BnB repair only when the number of removed pairs is small.
const BNB_REPAIR_MAX_PAIRS: usize = 8;
/// Maximum insertion options maintained per pair in BnB (paper-style top-y).
const BNB_REPAIR_TOP_Y: usize = 15;
/// Limited-discrepancy budget in BnB repair.
const BNB_REPAIR_MAX_DISCREPANCIES: usize = 3;
/// Per-repair wall-clock budget to keep LNS throughput stable.
const BNB_REPAIR_BUDGET_MS: u64 = 15;
/// Distance-mode stagnation: mid threshold.
#[allow(dead_code)]
const DIST_STAG_MID: usize = 183;
/// Distance-mode stagnation: high threshold.
#[allow(dead_code)]
const DIST_STAG_HIGH: usize = 314;
/// Focus stagnation trigger (iterations).
#[allow(dead_code)]
const FOCUS_STAG_TRIGGER: usize = 28;
/// Trust stagnation trigger (iterations).
#[allow(dead_code)]
const TRUST_STAG_TRIGGER: usize = 65;
/// Pool feeding period in distance mode (every N LNS iterations).
const POOL_ADD_PERIOD_DIST: usize = 26;
/// Pool feeding period in default (non-distance) mode.
const POOL_ADD_PERIOD_DEFAULT: usize = 48;

#[inline]
fn bnb_repair_enabled() -> bool {
    static ENABLED: OnceLock<bool> = OnceLock::new();
    *ENABLED.get_or_init(|| {
        std::env::var("BNB_REPAIR_ENABLE").ok().is_some_and(|s| {
            let v = s.trim().to_ascii_lowercase();
            !(v == "0" || v == "false" || v == "no" || v == "off")
        })
    })
}

#[derive(Clone, Copy)]
struct TabuEntry {
    signature: u64,
    expire_iter: usize,
}

#[derive(Clone)]
struct CachedRemoval {
    signature: u64,
    pickups: Vec<usize>,
    score: f64,
}

#[derive(Clone, Copy)]
enum BnbRepairMove {
    Existing {
        route_idx: usize,
        ep: usize,
        ed: usize,
        delta: f64,
    },
    NewRoute {
        delta: f64,
    },
}

impl BnbRepairMove {
    #[inline]
    fn delta(self) -> f64 {
        match self {
            Self::Existing { delta, .. } | Self::NewRoute { delta } => delta,
        }
    }
}

/// Minimal short-term memory for tabu search.
struct TabuList {
    entries: Vec<TabuEntry>,
}

impl TabuList {
    fn new(_tenure: usize) -> Self {
        Self {
            entries: Vec::new(),
        }
    }

    fn prune(&mut self, iter: usize) {
        self.entries.retain(|e| e.expire_iter > iter);
    }

    fn contains(&self, signature: u64, iter: usize) -> bool {
        self.entries
            .iter()
            .any(|e| e.signature == signature && e.expire_iter > iter)
    }

    fn add_with_tenure(&mut self, signature: u64, iter: usize, tenure: usize) {
        self.prune(iter);
        self.entries.retain(|e| e.signature != signature);
        self.entries.push(TabuEntry {
            signature,
            expire_iter: iter + tenure,
        });
    }
}

#[inline]
fn mix64(mut x: u64) -> u64 {
    x ^= x >> 30;
    x = x.wrapping_mul(0xbf58_476d_1ce4_e5b9);
    x ^= x >> 27;
    x = x.wrapping_mul(0x94d0_49bb_1331_11eb);
    x ^ (x >> 31)
}

/// Order-independent signature of a removal move.
#[inline]
fn tabu_move_signature(removed_pickups: &[usize]) -> u64 {
    let mut x = 0u64;
    let mut s = 0u64;
    for &p in removed_pickups {
        let m = mix64(p as u64);
        x ^= m;
        s = s.wrapping_add(m);
    }
    mix64(x.rotate_left(17) ^ s ^ (removed_pickups.len() as u64))
}

#[inline]
fn is_tabu_admissible(
    tabu: &TabuList,
    signature: u64,
    iter: usize,
    candidate_cost: f64,
    global_best_cost: f64,
) -> bool {
    !tabu.contains(signature, iter) || candidate_cost < global_best_cost - 1e-10
}

#[allow(clippy::too_many_arguments, clippy::cast_precision_loss)]
fn evaluate_tabu_candidate(
    inst: &Instance,
    base: &Solution,
    removed_template: &[usize],
    assigned_pickup: &[bool],
    removed_node_mask: &mut [bool],
    pickup_seen: &mut [bool],
    filtered_removed: &mut Vec<usize>,
    candidate: &mut Solution,
    rng: &mut impl Rng,
    regret_infos: &mut Vec<RouteInfo>,
    iter: usize,
) -> Option<(u64, f64)> {
    candidate.clone_from(base);
    filtered_removed.clear();

    for &pickup in removed_template {
        if pickup == 0
            || pickup > inst.n
            || !inst.is_pickup(pickup)
            || !assigned_pickup[pickup]
            || pickup_seen[pickup]
        {
            continue;
        }
        pickup_seen[pickup] = true;
        filtered_removed.push(pickup);
    }
    for &pickup in &*filtered_removed {
        pickup_seen[pickup] = false;
    }
    if filtered_removed.is_empty() {
        return None;
    }

    let signature = tabu_move_signature(filtered_removed);

    for &pickup in &*filtered_removed {
        removed_node_mask[pickup] = true;
        removed_node_mask[inst.delivery_of(pickup)] = true;
    }
    for route in &mut candidate.routes {
        route.retain(|&v| !removed_node_mask[v]);
    }
    for &pickup in &*filtered_removed {
        removed_node_mask[pickup] = false;
        removed_node_mask[inst.delivery_of(pickup)] = false;
    }

    candidate.routes.retain(|r| r.len() > 2);
    candidate
        .unassigned
        .extend(filtered_removed.iter().copied());

    let num_routes = candidate.routes.len();
    let k_options = [2, 3, 4, 5, num_routes.max(2)];
    let k = k_options[rng.random_range(0..k_options.len())];
    let noise_level = if rng.random_bool(0.60) { 0.025 } else { 0.0 };
    // Mostly sparse for speed, sometimes dense to escape A'^- confinement.
    let use_sparse = !iter.is_multiple_of(3);
    let total_dist = regret_insertion(
        inst,
        candidate,
        k,
        noise_level,
        use_sparse,
        rng,
        regret_infos,
        0,
        false,
    );
    let candidate_cost = candidate.routes.len() as f64 * 1e9 + total_dist;
    Some((signature, candidate_cost))
}

#[inline]
fn cache_upsert(cache: &mut Vec<CachedRemoval>, signature: u64, pickups: &[usize], score: f64) {
    if let Some(entry) = cache.iter_mut().find(|e| e.signature == signature) {
        // Exponential smoothing stabilizes noisy neighborhood estimates.
        entry.score = 0.65 * entry.score + 0.35 * score;
        entry.pickups.clear();
        entry.pickups.extend_from_slice(pickups);
        return;
    }

    cache.push(CachedRemoval {
        signature,
        pickups: pickups.to_vec(),
        score,
    });

    if cache.len() > TABU_CACHE_SIZE
        && let Some((worst_idx, _)) = cache
            .iter()
            .enumerate()
            .max_by(|a, b| a.1.score.partial_cmp(&b.1.score).unwrap())
    {
        cache.swap_remove(worst_idx);
    }
}

#[inline]
fn lns_obj_from_solution(inst: &Instance, sol: &Solution) -> f64 {
    #[allow(clippy::cast_precision_loss)]
    {
        sol.routes.len() as f64 * LNS_OBJ_BIG_M + sol.total_distance(inst)
    }
}

fn bnb_repair_options_for_pair(
    inst: &Instance,
    sol: &Solution,
    infos: &[RouteInfo],
    pickup: usize,
    top_y: usize,
) -> Vec<BnbRepairMove> {
    let delivery = inst.delivery_of(pickup);
    let mut options: Vec<BnbRepairMove> = Vec::new();

    for (ri, (route, info)) in sol.routes.iter().zip(infos.iter()).enumerate() {
        if inst.has_route_conflict(route, pickup) {
            continue;
        }
        if let Some((new_dist, ep, ed)) =
            best_pair_insertion_with_info::<false>(inst, route, info, pickup, delivery)
        {
            options.push(BnbRepairMove::Existing {
                route_idx: ri,
                ep,
                ed,
                delta: new_dist - info.distance,
            });
        }
    }

    if let Some(new_route_dist) = single_pair_feasible_distance(inst, pickup, delivery) {
        options.push(BnbRepairMove::NewRoute {
            delta: LNS_OBJ_BIG_M + new_route_dist,
        });
    }

    options.sort_by(|a, b| a.delta().partial_cmp(&b.delta()).unwrap());
    if options.len() > top_y {
        options.truncate(top_y);
    }
    options
}

fn bnb_select_pair_and_lb(
    inst: &Instance,
    sol: &Solution,
    remaining: &[usize],
    top_y: usize,
    infos: &mut Vec<RouteInfo>,
) -> Option<(usize, Vec<BnbRepairMove>, f64)> {
    compute_all_infos(inst, &sol.routes, infos);

    let mut chosen_idx: Option<usize> = None;
    let mut chosen_options: Vec<BnbRepairMove> = Vec::new();
    let mut worst_best_delta = f64::NEG_INFINITY;
    let mut lb_sum = 0.0;

    for (idx, &pickup) in remaining.iter().enumerate() {
        let options = bnb_repair_options_for_pair(inst, sol, infos, pickup, top_y);
        if options.is_empty() {
            return None;
        }

        let best_delta = options[0].delta();
        lb_sum += best_delta;

        if best_delta > worst_best_delta {
            worst_best_delta = best_delta;
            chosen_idx = Some(idx);
            chosen_options = options;
        }
    }

    chosen_idx.map(|idx| (idx, chosen_options, lb_sum))
}

#[allow(clippy::too_many_arguments)]
fn bnb_repair_dfs(
    inst: &Instance,
    sol: &Solution,
    remaining: &[usize],
    discrepancies: usize,
    max_discrepancies: usize,
    top_y: usize,
    deadline: Instant,
    infos: &mut Vec<RouteInfo>,
    best_cost: &mut f64,
    best_sol: &mut Option<Solution>,
) {
    if Instant::now() >= deadline {
        return;
    }

    if remaining.is_empty() {
        let cost = lns_obj_from_solution(inst, sol);
        if cost < *best_cost - 1e-10 {
            *best_cost = cost;
            *best_sol = Some(sol.clone());
        }
        return;
    }

    let current_cost = lns_obj_from_solution(inst, sol);
    if current_cost >= *best_cost - 1e-10 {
        return;
    }

    let Some((pair_idx, options, lb_delta)) =
        bnb_select_pair_and_lb(inst, sol, remaining, top_y, infos)
    else {
        return;
    };

    if current_cost + lb_delta >= *best_cost - 1e-10 {
        return;
    }

    let pickup = remaining[pair_idx];
    let delivery = inst.delivery_of(pickup);
    let best_delta = options[0].delta();

    for (rank, mv) in options.iter().enumerate() {
        if Instant::now() >= deadline {
            break;
        }

        let next_discrepancies = discrepancies + usize::from(rank > 0);
        if next_discrepancies > max_discrepancies {
            continue;
        }

        let optimistic = current_cost + mv.delta() + (lb_delta - best_delta);
        if optimistic >= *best_cost - 1e-10 {
            continue;
        }

        let mut child = sol.clone();
        match *mv {
            BnbRepairMove::Existing {
                route_idx, ep, ed, ..
            } => {
                child.routes[route_idx] =
                    build_route_with_pair(&child.routes[route_idx], pickup, ep, delivery, ed);
            }
            BnbRepairMove::NewRoute { .. } => {
                child.routes.push(vec![0, pickup, delivery, 0]);
            }
        }

        let mut child_remaining = remaining.to_vec();
        child_remaining.swap_remove(pair_idx);

        bnb_repair_dfs(
            inst,
            &child,
            &child_remaining,
            next_discrepancies,
            max_discrepancies,
            top_y,
            deadline,
            infos,
            best_cost,
            best_sol,
        );
    }
}

fn branch_and_bound_repair(inst: &Instance, sol: &mut Solution, deadline: Instant) -> bool {
    if sol.unassigned.is_empty() {
        return true;
    }

    let remaining = sol.unassigned.clone();
    let mut best_cost = f64::INFINITY;
    let mut best_sol: Option<Solution> = None;

    // Seed upper bound with deterministic regret repair for stronger pruning.
    let mut ub_sol = sol.clone();
    let mut ub_infos: Vec<RouteInfo> = Vec::new();
    let mut ub_rng = SmallRng::seed_from_u64(0xBADC_0FFE);
    regret_insertion(
        inst,
        &mut ub_sol,
        3,
        0.0,
        true,
        &mut ub_rng,
        &mut ub_infos,
        0,
        false,
    );
    if ub_sol.unassigned.is_empty() {
        best_cost = lns_obj_from_solution(inst, &ub_sol);
        best_sol = Some(ub_sol);
    }

    let mut root = sol.clone();
    root.unassigned.clear();
    let mut infos: Vec<RouteInfo> = Vec::new();
    bnb_repair_dfs(
        inst,
        &root,
        &remaining,
        0,
        BNB_REPAIR_MAX_DISCREPANCIES,
        BNB_REPAIR_TOP_Y,
        deadline,
        &mut infos,
        &mut best_cost,
        &mut best_sol,
    );

    if let Some(mut repaired) = best_sol {
        repaired.routes.retain(|r| r.len() > 2);
        repaired.unassigned.clear();
        *sol = repaired;
        true
    } else {
        false
    }
}

/// Large Neighbourhood Search with adaptive acceptance (LAHC/SA/RRT).
/// Section 3.3 of the paper, extended with AHGSLNS adaptive criterion selection.
/// Updates `best` if a better solution is found during LNS.
///
/// `accept_state`: persistent acceptance state across calls. Pass `None` for a
/// one-shot call (e.g. parent_b generation); pass `Some(&mut state)` in
/// the main loop so that acceptance history accumulates across iterations.
pub fn lns(
    inst: &Instance,
    sol: &mut Solution,
    best: &mut Solution,
    ges_duration: Duration,
    rng: &mut impl Rng,
    deadline: Instant,
    accept_state: Option<&mut AcceptanceState>,
    alns_state: Option<&mut AlnsState>,
    arc_history: Option<&ArcHistory>,
    mut penalty: Option<&mut PenaltyState>,
    mut route_pool: Option<&mut RoutePool>,
) -> usize {
    let start = Instant::now();
    let mut lns_iterations: usize = 0;
    let min_time = if inst.m >= 200 {
        (ges_duration * 5 / 4).max(Duration::from_secs(3))
    } else {
        ges_duration * 2
    };

    // Initialize or restore acceptance state (LAHC/SA/RRT).
    // SA/LAHC/RRT always operate on the distance scale. Vehicle count and
    // unassigned changes are handled by a hierarchical gate in the inner loop.
    let distance_mode = sol.unassigned.is_empty()
        && best.unassigned.is_empty()
        && sol.num_vehicles() == best.num_vehicles();
    let initial_cost = sol.total_distance(inst);
    let mut local_accept = AcceptanceState::new(initial_cost);
    let accept = match accept_state {
        Some(s) => {
            if s.lahc.costs.is_empty() {
                // First call or after reset: initialize all criteria
                s.init(initial_cost);
            } else {
                // Continuing: sync current_cost with the actual sol cost.
                s.current_cost = initial_cost;
                // Reheat RRT deviation if too small.
                let target_d0 = 0.02 * initial_cost;
                if s.rrt.deviation < target_d0 * 0.05 {
                    s.rrt.deviation = target_d0 * 0.5;
                    s.rrt.d0 = target_d0;
                }
                // Update RRT record to current best distance.
                if initial_cost < s.rrt.record {
                    s.rrt.record = initial_cost;
                }
            }
            s
        }
        None => {
            // One-shot call: use local state (defaults to LAHC)
            &mut local_accept
        }
    };
    // Initialize or restore ALNS state
    let mut local_alns = AlnsState::new();
    let alns = match alns_state {
        Some(s) => s,
        None => &mut local_alns,
    };
    let mut no_improve_count: usize = 0;
    let mut best_cost = { best.total_distance(inst) };

    // Pre-allocate candidate solution (reused across iterations via clone_from)
    let mut candidate = sol.clone();
    // Reusable bitset for pair removal
    let mut removed_set = vec![false; inst.n + 1];
    // Pre-allocate RouteInfo vec for regret insertion (reused across iterations)
    let mut regret_infos: Vec<RouteInfo> = Vec::new();

    // Adaptive min removal: for long routes (loose TW), small removals barely
    // change structure. Scale min_remove with avg route length.
    let total_customers: usize = sol.routes.iter().map(|r| r.len().saturating_sub(2)).sum();
    let avg_route_len = total_customers / sol.routes.len().max(1);
    let adaptive_min_remove = (avg_route_len / 8).max(MIN_REMOVE);

    loop {
        if Instant::now() >= deadline {
            break;
        }
        // Stopping conditions (Section 3.3):
        if no_improve_count >= MIN_NO_IMPROVE {
            break;
        }
        if no_improve_count >= ALT_NO_IMPROVE && start.elapsed() >= min_time {
            break;
        }

        lns_iterations += 1;

        // Reuse candidate allocation instead of sol.clone() each iteration
        candidate.clone_from(sol);

        // O(1) count of assigned pairs
        let num_assigned = inst.m - candidate.unassigned.len();
        // Scale max removal with instance size. In distance mode (vehicle count
        // stable), allow up to 25% pair removal for structural changes.
        // Otherwise cap at 33% to preserve vehicle count.
        let effective_max_remove = if distance_mode {
            MAX_REMOVE.max(inst.m / 4)
        } else {
            MAX_REMOVE.max(inst.m / 3)
        };
        let num_remove = rng.random_range(
            adaptive_min_remove..=effective_max_remove.min(num_assigned.max(adaptive_min_remove)),
        );

        let destroy_op = alns.select_destroy(rng);
        let removed = match destroy_op {
            DestroyOp::RouteElimination => {
                route_elimination_removal(inst, &candidate, num_remove, rng)
            }
            DestroyOp::Shaw => shaw_removal(inst, &candidate, num_remove, rng),
            DestroyOp::WorstDistance => worst_removal(inst, &candidate, num_remove, rng),
            DestroyOp::Cluster => cluster_removal(inst, &candidate, num_remove, rng),
            DestroyOp::StringSisrs => string_removal(inst, &candidate, num_remove, rng),
            DestroyOp::RcRemoval => rc_removal(inst, &candidate, num_remove, rng),
            DestroyOp::Random => random_removal(inst, &candidate, num_remove, rng),
            DestroyOp::HistoryRemoval => match arc_history {
                Some(h) if !h.is_empty() => history_removal(inst, &candidate, num_remove, rng, h),
                _ => shaw_removal(inst, &candidate, num_remove, rng),
            },
            DestroyOp::LateArrival => late_arrival_removal(inst, &candidate, num_remove, rng),
            DestroyOp::DetourMisplacement => {
                detour_misplacement_removal(inst, &candidate, num_remove, rng)
            }
            DestroyOp::ChainedRuin => chained_ruin_removal(inst, &candidate, num_remove, rng),
        };

        // Bitset-based removal: O(n) instead of O(removed * routes * route_len)
        for &pickup in &removed {
            removed_set[pickup] = true;
            removed_set[inst.delivery_of(pickup)] = true;
        }
        for route in &mut candidate.routes {
            route.retain(|&v| !removed_set[v]);
        }
        // Clear bitset for reuse (only clear entries we set)
        for &pickup in &removed {
            removed_set[pickup] = false;
            removed_set[inst.delivery_of(pickup)] = false;
        }
        // Clean up empty routes
        candidate.routes.retain(|r| r.len() > 2);
        candidate.unassigned.extend(removed);

        // Re-insert removed pairs using either regret heuristic or a bounded
        // branch-and-bound repair (paper-style top-y + limited discrepancy).
        let repair_op = alns.select_repair(rng);
        let repair_k = repair_op.to_k(candidate.routes.len().max(1), rng);
        let use_id = repair_op.use_insertion_difficulty();

        let use_bnb_repair = bnb_repair_enabled()
            && candidate.unassigned.len() >= 2
            && candidate.unassigned.len() <= BNB_REPAIR_MAX_PAIRS
            && rng.random_bool(BNB_REPAIR_PROB);

        let total_dist = if use_bnb_repair {
            let now = Instant::now();
            let bnb_budget = Duration::from_millis(BNB_REPAIR_BUDGET_MS)
                .min(deadline.saturating_duration_since(now));
            let bnb_deadline = now + bnb_budget;
            if bnb_budget > Duration::ZERO
                && branch_and_bound_repair(inst, &mut candidate, bnb_deadline)
            {
                candidate.total_distance(inst)
            } else {
                let noise_level = if rng.random_bool(0.5) { 0.10 } else { 0.0 };
                regret_insertion(
                    inst,
                    &mut candidate,
                    repair_k,
                    noise_level,
                    true,
                    rng,
                    &mut regret_infos,
                    0,
                    use_id,
                )
            }
        } else {
            let noise_level = if rng.random_bool(0.5) { 0.10 } else { 0.0 };
            regret_insertion(
                inst,
                &mut candidate,
                repair_k,
                noise_level,
                true,
                rng,
                &mut regret_infos,
                0,
                use_id,
            )
        };

        // Candidate cost: always on distance scale so SA/LAHC/RRT thresholds
        // remain meaningful. Vehicle/unassigned changes are handled by a
        // hierarchical gate before the acceptance criterion.
        // Capacity and TW penalties are added when the penalty mechanism is active.
        let candidate_cost = {
            let base = total_dist;
            if let Some(ref pen) = penalty {
                if pen.is_active() {
                    let cap_excess: f64 = regret_infos
                        .iter()
                        .take(candidate.routes.len())
                        .map(|info| info.cap_excess)
                        .sum();
                    let tw_excess: f64 = regret_infos
                        .iter()
                        .take(candidate.routes.len())
                        .map(|info| info.tw_excess)
                        .sum();
                    base + pen.penalty_cap * cap_excess + pen.penalty_tw * tw_excess
                } else {
                    base
                }
            } else {
                base
            }
        };

        // Save costs before acceptance for ALNS scoring
        let prev_current_cost = accept.current_cost;
        let prev_best_cost = best_cost;

        // Hierarchical gate: vehicle count and unassigned changes are decided
        // deterministically, without going through SA/LAHC/RRT (whose thresholds
        // are calibrated on the distance scale and would be nonsensical for the
        // composite vehicles*1e9 scale).
        let cand_veh = candidate.routes.len();
        let cand_unasgn = candidate.unassigned.len();
        let sol_veh = sol.routes.len();
        let sol_unasgn = sol.unassigned.len();

        // More unassigned → always reject.
        // Same unassigned but more vehicles → reject.
        let dominated =
            cand_unasgn > sol_unasgn || (cand_unasgn == sol_unasgn && cand_veh > sol_veh);

        let accepted = if dominated {
            false
        } else if cand_unasgn < sol_unasgn || cand_veh < sol_veh {
            // Strictly better on hierarchy → always accept
            true
        } else {
            // Same vehicles and unassigned: accept by LAHC/RRT on distance scale
            match accept.active {
                AcceptanceCriterion::Lahc => {
                    let x = accept.lahc.iter % LHC_LEN;
                    candidate_cost <= accept.lahc.costs[x] || candidate_cost <= accept.current_cost
                }
                AcceptanceCriterion::Rrt => {
                    candidate_cost <= accept.rrt.record + accept.rrt.deviation
                }
            }
        };
        if accepted {
            std::mem::swap(sol, &mut candidate);
            accept.current_cost = candidate_cost;
            // In distance mode, run fast intra-route LS (2-Opt + Or-Opt)
            // on every accepted solution. This fixes sub-optimal route
            // orderings left by regret insertion, keeping the working
            // solution polished for future destroy-repair cycles.
            if distance_mode && sol.unassigned.is_empty() {
                local_search::local_search_intra(inst, sol);
                accept.current_cost = sol.total_distance(inst);
            }
        }
        // Always update LAHC history (prevents stale history on criterion switch-back)
        let x = accept.lahc.iter % LHC_LEN;
        accept.lahc.costs[x] = accept.current_cost;
        accept.lahc.iter += 1;
        // RRT decay (applied every inner iteration)
        accept.rrt.decay();
        // Update RRT record on improvement
        if accept.current_cost < accept.rrt.record {
            accept.rrt.record = accept.current_cost;
        }

        // Update best solution.
        // In distance mode, apply basic LS when a new best is found to polish
        // the solution before recording it. This is rare (~1-5 times per round)
        // so the ~10ms LS cost is amortized.
        let hierarchy_ok = sol.unassigned.len() <= best.unassigned.len()
            && (sol.unassigned.len() < best.unassigned.len()
                || sol.num_vehicles() <= best.num_vehicles());
        if hierarchy_ok && accept.current_cost < best_cost && sol.is_feasible(inst) {
            if distance_mode {
                let disabled = PenaltyState::disabled();
                local_search::local_search(inst, sol, &disabled, None);
                let new_cost = sol.total_distance(inst);
                accept.current_cost = new_cost;
            }
            best.clone_from(sol);
            best_cost = accept.current_cost;
            no_improve_count = 0;
        } else {
            no_improve_count += 1;
        }

        // Adaptive penalty: record feasibility using cached regret_infos
        if let Some(ref mut pen) = penalty
            && pen.is_active()
        {
            let nr = sol.routes.len();
            let cap_feas = regret_infos
                .iter()
                .take(nr)
                .all(|info| info.cap_excess < 1e-10);
            let tw_feas = regret_infos
                .iter()
                .take(nr)
                .all(|info| info.tw_excess < 1e-10);
            pen.record_cap_feasibility(cap_feas);
            pen.record_tw_feasibility(tw_feas);
            pen.maybe_adjust_cap(100);
            pen.maybe_adjust_tw(100);
        }

        // ALNS scoring (Ropke & Pisinger 2006, sigma1/sigma2/sigma3)
        let score = if best_cost < prev_best_cost - 1e-10 {
            33.0 // new global best
        } else if accepted && candidate_cost < prev_current_cost - 1e-10 {
            9.0 // improved current solution
        } else if accepted {
            13.0 // accepted worse/equal
        } else {
            0.0 // rejected
        };
        alns.record_outcome(destroy_op, repair_op, score);

        // Periodically feed accepted candidate routes into SC pool.
        // Every 50 iterations, if the current solution is feasible, add its
        // routes to the pool. This gives SC access to thousands of diverse
        // routes that LNS produces but would otherwise discard.
        let pool_period = if distance_mode {
            POOL_ADD_PERIOD_DIST
        } else {
            POOL_ADD_PERIOD_DEFAULT
        };
        if accepted
            && lns_iterations.is_multiple_of(pool_period)
            && sol.is_feasible(inst)
            && let Some(ref mut pool) = route_pool
        {
            pool.add_solution(inst, sol);
        }
    }
    lns_iterations
}

/// Short tabu-search kick used only during stagnation.
/// Uses LNS-style neighborhood generation with tabu memory + aspiration.
pub fn tabu_search_stagnation(
    inst: &Instance,
    sol: &mut Solution,
    best: &mut Solution,
    rng: &mut impl Rng,
    deadline: Instant,
) -> bool {
    let mut improved = false;
    let mut tabu = TabuList::new(TABU_TENURE);
    let mut global_best_cost = best.cost(inst);
    let mut move_cache: Vec<CachedRemoval> = Vec::new();

    let mut candidate = sol.clone();
    let mut best_neighbor = sol.clone();
    let mut best_signature = 0u64;
    let mut removed_set = vec![false; inst.n + 1];
    let mut pickup_seen = vec![false; inst.n + 1];
    let mut assigned_pickup = vec![false; inst.n + 1];
    let mut filtered_removed: Vec<usize> = Vec::new();
    let mut regret_infos: Vec<RouteInfo> = Vec::new();

    for iter in 0..TABU_MAX_ITERS {
        if Instant::now() >= deadline {
            break;
        }

        assigned_pickup.fill(false);
        for route in &sol.routes {
            for &v in route {
                if v != 0 && inst.is_pickup(v) {
                    assigned_pickup[v] = true;
                }
            }
        }

        let current_cost = sol.cost(inst);
        tabu.prune(iter);
        let mut found_neighbor = false;
        let mut best_neighbor_cost = f64::MAX;

        move_cache.sort_by(|a, b| a.score.partial_cmp(&b.score).unwrap());
        let reevaluate = move_cache.len().min(TABU_CACHE_REEVAL);
        for entry in move_cache.iter_mut().take(reevaluate) {
            if Instant::now() >= deadline {
                break;
            }

            let cached_pickups = entry.pickups.clone();
            let Some((signature, candidate_cost)) = evaluate_tabu_candidate(
                inst,
                sol,
                &cached_pickups,
                &assigned_pickup,
                &mut removed_set,
                &mut pickup_seen,
                &mut filtered_removed,
                &mut candidate,
                rng,
                &mut regret_infos,
                iter,
            ) else {
                entry.score = f64::INFINITY;
                continue;
            };

            entry.signature = signature;
            entry.score = 0.65 * entry.score + 0.35 * candidate_cost;
            entry.pickups.clear();
            entry.pickups.extend_from_slice(&filtered_removed);

            if !is_tabu_admissible(&tabu, signature, iter, candidate_cost, global_best_cost) {
                continue;
            }
            if candidate_cost < best_neighbor_cost {
                best_neighbor.clone_from(&candidate);
                best_neighbor_cost = candidate_cost;
                best_signature = signature;
                found_neighbor = true;
            }
        }

        let num_assigned = inst.m - sol.unassigned.len();
        for _ in 0..TABU_RANDOM_NEIGHBORS {
            if Instant::now() >= deadline {
                break;
            }
            if num_assigned == 0 {
                continue;
            }

            let min_remove = MIN_REMOVE.max((inst.m / 40).max(6));
            let max_remove = MAX_REMOVE.max(min_remove);
            let upper = max_remove.min(num_assigned.max(min_remove));
            let num_remove = rng.random_range(min_remove..=upper);

            let r: f64 = rng.random();
            let removed = if r < ROUTE_ELIM_PROB {
                route_elimination_removal(inst, sol, num_remove, rng)
            } else if r < ROUTE_ELIM_PROB + SHAW_PROB {
                shaw_removal(inst, sol, num_remove, rng)
            } else if r < ROUTE_ELIM_PROB + SHAW_PROB + WORST_PROB {
                worst_removal(inst, sol, num_remove, rng)
            } else if r < ROUTE_ELIM_PROB + SHAW_PROB + WORST_PROB + CLUSTER_PROB {
                cluster_removal(inst, sol, num_remove, rng)
            } else if r < ROUTE_ELIM_PROB + SHAW_PROB + WORST_PROB + CLUSTER_PROB + STRING_PROB {
                string_removal(inst, sol, num_remove, rng)
            } else {
                random_removal(inst, sol, num_remove, rng)
            };
            let Some((signature, candidate_cost)) = evaluate_tabu_candidate(
                inst,
                sol,
                &removed,
                &assigned_pickup,
                &mut removed_set,
                &mut pickup_seen,
                &mut filtered_removed,
                &mut candidate,
                rng,
                &mut regret_infos,
                iter,
            ) else {
                continue;
            };

            cache_upsert(
                &mut move_cache,
                signature,
                &filtered_removed,
                candidate_cost,
            );

            if !is_tabu_admissible(&tabu, signature, iter, candidate_cost, global_best_cost) {
                continue;
            }
            if candidate_cost < best_neighbor_cost {
                best_neighbor.clone_from(&candidate);
                best_neighbor_cost = candidate_cost;
                best_signature = signature;
                found_neighbor = true;
            }
        }

        move_cache.retain(|e| e.score.is_finite() && !e.pickups.is_empty());

        if !found_neighbor {
            break;
        }

        let move_improves_current = best_neighbor_cost < current_cost - 1e-10;
        sol.clone_from(&best_neighbor);
        let tenure = if move_improves_current {
            rng.random_range(TABU_TENURE..=(TABU_TENURE * 2))
        } else {
            rng.random_range(TABU_TENURE_MIN..=TABU_TENURE)
        };
        tabu.add_with_tenure(best_signature, iter, tenure);

        if best_neighbor_cost < global_best_cost && sol.is_feasible(inst) {
            best.clone_from(sol);
            global_best_cost = best_neighbor_cost;
            improved = true;
        }
    }

    improved
}

/// Shaw removal: remove `count` similar PD pairs.
/// Similarity = phi1*dist + phi2*|`arrival_diff`| + phi3*|`demand_diff`|
/// (Ropke and Pisinger 2006, as cited in paper Section 3.3)
#[allow(clippy::similar_names)]
fn shaw_removal(inst: &Instance, sol: &Solution, count: usize, rng: &mut impl Rng) -> Vec<usize> {
    let phi1: f64 = 9.0;
    let phi2: f64 = 3.0;
    let phi3: f64 = 2.0;
    let phi4: f64 = 4.0; // LP reduced cost similarity weight

    // Collect all assigned pickups with their route arrival times and RC scores
    let mut assigned: Vec<usize> = Vec::new();
    let mut arrival_time: Vec<f64> = vec![0.0; inst.n + 1];
    let mut rc_score: Vec<f64> = vec![0.0; inst.n + 1]; // per-pickup RC score

    let mut pos = vec![0usize; inst.n + 1];
    for route in &sol.routes {
        for (i, &v) in route.iter().enumerate() {
            pos[v] = i;
        }
    }

    for route in &sol.routes {
        let mut time = inst.early(0);
        for i in 1..route.len() {
            let prev = route[i - 1];
            let curr = route[i];
            time = crate::solution::fmax(
                time + inst.svc(prev) + inst.dist(prev, curr),
                inst.early(curr),
            );
            arrival_time[curr] = time;
            if curr != 0 && inst.is_pickup(curr) {
                assigned.push(curr);
                // RC score: sum of reduced costs on arcs incident to pickup and delivery
                let d = inst.delivery_of(curr);
                let prev_p = route[pos[curr].saturating_sub(1)];
                let next_p = route[(pos[curr] + 1).min(route.len() - 1)];
                let prev_d = route[pos[d].saturating_sub(1)];
                let next_d = route[(pos[d] + 1).min(route.len() - 1)];
                let rc = inst.rc(prev_p, curr).min(1e12)
                    + inst.rc(curr, next_p).min(1e12)
                    + inst.rc(prev_d, d).min(1e12)
                    + inst.rc(d, next_d).min(1e12);
                rc_score[curr] = rc;
            }
        }
    }

    if assigned.is_empty() {
        return Vec::new();
    }

    // Normalize RC scores for comparability with other similarity terms
    let max_rc = rc_score.iter().copied().fold(0.0f64, f64::max);
    let rc_norm = if max_rc > 1e-10 { max_rc } else { 1.0 };

    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];

    // Randomly select first pair
    let first = assigned[rng.random_range(0..assigned.len())];
    removed.push(first);
    in_removed[first] = true;

    // Pre-allocate similarity buffer outside loop to avoid per-iteration allocation
    let mut similarities: Vec<(f64, usize)> = Vec::with_capacity(assigned.len());

    while removed.len() < count && removed.len() < assigned.len() {
        // Select a random pair from removed set as reference
        let ref_pickup = removed[rng.random_range(0..removed.len())];

        // Compute similarity to all non-removed assigned pairs
        similarities.clear();
        for &p in &assigned {
            if in_removed[p] {
                continue;
            }
            let dist_sim = inst.dist(ref_pickup, p);
            let time_sim = (arrival_time[ref_pickup] - arrival_time[p]).abs();
            let demand_sim = f64::from((inst.demand[ref_pickup] - inst.demand[p]).abs());
            // RC similarity: pairs on arcs with similar reduced costs should be removed together
            let rc_sim = (rc_score[ref_pickup] - rc_score[p]).abs() / rc_norm;
            let sim = phi1 * dist_sim + phi2 * time_sim + phi3 * demand_sim + phi4 * rc_sim;
            similarities.push((sim, p));
        }

        if similarities.is_empty() {
            break;
        }

        // Select with randomization: pick index = floor(random^p * len)
        // where p is the randomization parameter (typically 6 from Ropke & Pisinger)
        // Use select_nth_unstable_by (O(n)) instead of full sort (O(n log n)).
        let p_rand: f64 = 6.0;
        #[allow(
            clippy::cast_possible_truncation,
            clippy::cast_sign_loss,
            clippy::cast_precision_loss
        )]
        let idx = (rng.random::<f64>().powf(p_rand) * similarities.len() as f64) as usize;
        let idx = idx.min(similarities.len() - 1);
        similarities.select_nth_unstable_by(idx, |a, b| a.0.partial_cmp(&b.0).unwrap());
        let selected = similarities[idx].1;

        removed.push(selected);
        in_removed[selected] = true;
    }

    removed
}

/// Chained Ruin: multiple small focused ruins in different areas of the solution,
/// each expanded with Shaw-similar pairs. Creates spatially diverse but locally
/// coherent destruction patterns.
///
/// Algorithm:
/// 1. Pre-compute assigned pickups, positions, node_route, arrival_time arrays
/// 2. Pick 2-4 chain seeds from different routes (spatial diversity)
/// 3. Per chain: core ruin (60% budget) via random sub-method, then Shaw expansion (40% budget)
/// 4. Return Vec<usize> of pickup IDs
#[allow(
    clippy::cast_possible_truncation,
    clippy::cast_sign_loss,
    clippy::cast_precision_loss
)]
fn chained_ruin_removal(
    inst: &Instance,
    sol: &Solution,
    count: usize,
    rng: &mut impl Rng,
) -> Vec<usize> {
    use crate::solution::fmax;

    // Fallback for tiny instances
    if count < 4 {
        return shaw_removal(inst, sol, count, rng);
    }

    // --- Pre-compute arrays ---
    let mut assigned: Vec<usize> = Vec::new();
    let mut pos = vec![0usize; inst.n + 1];
    let mut node_route: Vec<u16> = vec![0; inst.n + 1];
    let mut arrival_time = vec![0.0f64; inst.n + 1];

    for (ri, route) in sol.routes.iter().enumerate() {
        let mut time = inst.early(0);
        for (i, &v) in route.iter().enumerate() {
            pos[v] = i;
            node_route[v] = ri as u16;
            if i > 0 {
                let prev = route[i - 1];
                time = fmax(time + inst.svc(prev) + inst.dist(prev, v), inst.early(v));
            }
            arrival_time[v] = time;
            if v != 0 && inst.is_pickup(v) {
                assigned.push(v);
            }
        }
    }

    if assigned.is_empty() {
        return Vec::new();
    }

    // Build per-route pickup lists for seed selection from different routes
    let num_routes = sol.routes.len();
    let mut route_pickups: Vec<Vec<usize>> = vec![Vec::new(); num_routes];
    for &p in &assigned {
        route_pickups[node_route[p] as usize].push(p);
    }
    // Non-empty route indices, shuffled for random seed selection
    let mut nonempty_routes: Vec<usize> = (0..num_routes)
        .filter(|&ri| !route_pickups[ri].is_empty())
        .collect();
    nonempty_routes.shuffle(rng);

    // Pick 2-4 chains from different routes
    let num_chains = rng.random_range(2..=4).min(nonempty_routes.len());
    if num_chains == 0 {
        return shaw_removal(inst, sol, count, rng);
    }

    let budget_per_chain = count / num_chains;
    let mut extra = count - budget_per_chain * num_chains; // distribute remainder

    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];

    // Pre-compute detour costs for worst-detour sub-method
    let mut detour_cost = vec![0.0f64; inst.n + 1];
    for route in &sol.routes {
        for &v in &route[1..route.len() - 1] {
            if v == 0 || !inst.is_pickup(v) {
                continue;
            }
            let d = inst.delivery_of(v);
            let ri = node_route[v] as usize;
            let prev_p = sol.routes[ri][pos[v] - 1];
            let next_p = sol.routes[ri][pos[v] + 1];
            let cost_p = inst.dist(prev_p, v) + inst.dist(v, next_p) - inst.dist(prev_p, next_p);
            let prev_d = sol.routes[ri][pos[d] - 1];
            let next_d = sol.routes[ri][pos[d] + 1];
            let cost_d = inst.dist(prev_d, d) + inst.dist(d, next_d) - inst.dist(prev_d, next_d);
            detour_cost[v] = cost_p + cost_d;
        }
    }

    // Pre-compute late-arrival scores
    let mut late_score = vec![0.0f64; inst.n + 1];
    for &p in &assigned {
        let d = inst.delivery_of(p);
        late_score[p] = (arrival_time[p] - inst.early(p)) + (arrival_time[d] - inst.early(d));
    }

    for &route_idx in nonempty_routes.iter().take(num_chains) {
        let chain_budget = budget_per_chain
            + if extra > 0 {
                extra -= 1;
                1
            } else {
                0
            };
        if chain_budget == 0 {
            continue;
        }

        // Pick seed from this chain's route
        let candidates: Vec<usize> = route_pickups[route_idx]
            .iter()
            .copied()
            .filter(|&p| !in_removed[p])
            .collect();
        if candidates.is_empty() {
            continue;
        }
        let seed = candidates[rng.random_range(0..candidates.len())];

        // --- Core ruin: 60% of chain budget ---
        let core_budget = (chain_budget * 3 / 5).max(1);
        let shaw_budget = chain_budget - core_budget;

        // Add seed
        if !in_removed[seed] {
            removed.push(seed);
            in_removed[seed] = true;
        }

        // Choose sub-method randomly: 0=worst-detour, 1=late-arrival, 2=cluster, 3=random-noisy
        let sub_method = rng.random_range(0..4u32);
        let core_p_rand: f64 = 4.0;

        // Score non-removed assigned pairs relative to seed based on sub-method
        let mut scored: Vec<(f64, usize)> = Vec::with_capacity(assigned.len());
        for &p in &assigned {
            if in_removed[p] {
                continue;
            }
            let score = match sub_method {
                0 => detour_cost[p], // worst-detour: higher detour → remove first
                1 => late_score[p],  // late-arrival: higher slack → remove first
                2 => {
                    // cluster around seed: smaller distance → higher priority (invert for descending)
                    let d_seed = inst.delivery_of(seed);
                    let d_p = inst.delivery_of(p);
                    -(inst.dist(seed, p) + inst.dist(d_seed, d_p))
                }
                _ => rng.random::<f64>(), // random-noisy
            };
            scored.push((score, p));
        }

        // Select core_budget - 1 pairs (seed already added) by score descending with randomization
        let mut core_added = 1usize; // seed already counted
        while core_added < core_budget && !scored.is_empty() {
            let idx = (rng.random::<f64>().powf(core_p_rand) * scored.len() as f64) as usize;
            let idx = idx.min(scored.len() - 1);
            scored.select_nth_unstable_by(idx, |a, b| b.0.partial_cmp(&a.0).unwrap());
            let pickup = scored[idx].1;
            if !in_removed[pickup] {
                removed.push(pickup);
                in_removed[pickup] = true;
                core_added += 1;
            }
            scored.swap_remove(idx);
        }

        // --- Shaw expansion: 40% of chain budget ---
        // Reference = only this chain's removals (locally coherent)
        if shaw_budget == 0 {
            continue;
        }
        let chain_start = removed.len() - core_added;
        let shaw_phi1: f64 = 9.0;
        let shaw_phi2: f64 = 3.0;
        let shaw_phi3: f64 = 2.0;
        let shaw_p_rand: f64 = 6.0;

        let mut similarities: Vec<(f64, usize)> = Vec::with_capacity(assigned.len());
        let mut shaw_added = 0usize;

        while shaw_added < shaw_budget {
            // Pick reference from this chain's removals only
            let ref_pickup = removed[rng.random_range(chain_start..removed.len())];

            similarities.clear();
            for &p in &assigned {
                if in_removed[p] {
                    continue;
                }
                let dist_sim = inst.dist(ref_pickup, p);
                let time_sim = (arrival_time[ref_pickup] - arrival_time[p]).abs();
                let demand_sim = f64::from((inst.demand[ref_pickup] - inst.demand[p]).abs());
                let sim = shaw_phi1 * dist_sim + shaw_phi2 * time_sim + shaw_phi3 * demand_sim;
                similarities.push((sim, p));
            }

            if similarities.is_empty() {
                break;
            }

            let idx = (rng.random::<f64>().powf(shaw_p_rand) * similarities.len() as f64) as usize;
            let idx = idx.min(similarities.len() - 1);
            similarities.select_nth_unstable_by(idx, |a, b| a.0.partial_cmp(&b.0).unwrap());
            let selected = similarities[idx].1;

            if !in_removed[selected] {
                removed.push(selected);
                in_removed[selected] = true;
                shaw_added += 1;
            } else {
                break; // avoid infinite loop
            }
        }
    }

    removed
}

/// Cluster removal: remove PD pairs that are geographically close to a random
/// seed pair. Uses pair-level proximity: dist(seed_p, p) + dist(seed_d, d).
/// Randomization exponent controls how tightly the cluster is selected.
#[allow(
    clippy::cast_possible_truncation,
    clippy::cast_sign_loss,
    clippy::cast_precision_loss
)]
fn cluster_removal(
    inst: &Instance,
    sol: &Solution,
    count: usize,
    rng: &mut impl Rng,
) -> Vec<usize> {
    let mut assigned: Vec<usize> = Vec::new();
    for route in &sol.routes {
        for &v in route {
            if v != 0 && inst.is_pickup(v) {
                assigned.push(v);
            }
        }
    }
    if assigned.is_empty() {
        return Vec::new();
    }

    let seed = assigned[rng.random_range(0..assigned.len())];
    let seed_d = inst.delivery_of(seed);

    // Compute pair-level proximity to seed and sort
    let mut dists: Vec<(f64, usize)> = assigned
        .iter()
        .filter(|&&p| p != seed)
        .map(|&p| {
            let d = inst.dist(seed, p) + inst.dist(seed_d, inst.delivery_of(p));
            (d, p)
        })
        .collect();

    let mut removed = Vec::with_capacity(count);
    removed.push(seed);

    // Select with randomization (exponent p=4 — tighter cluster than Shaw's p=6)
    let p_rand: f64 = 4.0;
    for _ in 1..count {
        if dists.is_empty() {
            break;
        }
        let idx = (rng.random::<f64>().powf(p_rand) * dists.len() as f64) as usize;
        let idx = idx.min(dists.len() - 1);
        dists.select_nth_unstable_by(idx, |a, b| a.0.partial_cmp(&b.0).unwrap());
        removed.push(dists[idx].1);
        dists.swap_remove(idx);
    }

    removed
}

/// Random removal: randomly select `count` PD pairs to remove.
fn random_removal(inst: &Instance, sol: &Solution, count: usize, rng: &mut impl Rng) -> Vec<usize> {
    let mut assigned: Vec<usize> = Vec::new();
    for route in &sol.routes {
        for &v in route {
            if v != 0 && inst.is_pickup(v) {
                assigned.push(v);
            }
        }
    }

    assigned.shuffle(rng);
    assigned.truncate(count);
    assigned
}

/// History removal (Masson et al. 2012, §3.2 heuristic 4).
/// Remove requests whose arcs rarely appear in the best-known solutions.
/// Low φ score = poorly placed relative to the search history.
#[allow(
    clippy::cast_possible_truncation,
    clippy::cast_sign_loss,
    clippy::cast_precision_loss
)]
fn history_removal(
    inst: &Instance,
    sol: &Solution,
    count: usize,
    rng: &mut impl Rng,
    history: &ArcHistory,
) -> Vec<usize> {
    let mut pos = vec![0usize; inst.n + 1];
    let mut node_route: Vec<u16> = vec![0; inst.n + 1];
    for (ri, route) in sol.routes.iter().enumerate() {
        for (i, &v) in route.iter().enumerate() {
            pos[v] = i;
            node_route[v] = ri as u16;
        }
    }

    // Compute φ_i for each assigned pickup: sum of arc frequencies
    // on the 4 arcs incident to its pickup and delivery nodes.
    let mut scores: Vec<(f64, usize)> = Vec::new();
    for route in &sol.routes {
        for &v in &route[1..route.len() - 1] {
            if v == 0 || !inst.is_pickup(v) {
                continue;
            }
            let d = inst.delivery_of(v);
            let ri = node_route[v] as usize;

            let prev_p = sol.routes[ri][pos[v] - 1];
            let next_p = sol.routes[ri][pos[v] + 1];
            let prev_d = sol.routes[ri][pos[d] - 1];
            let next_d = sol.routes[ri][pos[d] + 1];

            let phi = history.get(prev_p, v)
                + history.get(v, next_p)
                + history.get(prev_d, d)
                + history.get(d, next_d);

            scores.push((f64::from(phi), v));
        }
    }

    if scores.is_empty() {
        return Vec::new();
    }

    // Remove requests with LOWEST φ (poorly placed relative to best solutions)
    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];
    let p_rand: f64 = 3.0;

    while removed.len() < count && !scores.is_empty() {
        let idx = (rng.random::<f64>().powf(p_rand) * scores.len() as f64) as usize;
        let idx = idx.min(scores.len() - 1);

        // Ascending sort: lowest φ (worst-placed) first
        scores.select_nth_unstable_by(idx, |a, b| a.0.partial_cmp(&b.0).unwrap());
        let pickup = scores[idx].1;

        if !in_removed[pickup] {
            removed.push(pickup);
            in_removed[pickup] = true;
        }
        scores.swap_remove(idx);
    }

    removed
}

/// Worst-distance removal: remove pairs contributing most to total distance.
/// Cost(pair) = detour of pickup + detour of delivery in the route.
/// Uses position lookup (O(n) precompute) and select_nth for O(n) selection.
fn worst_removal(inst: &Instance, sol: &Solution, count: usize, rng: &mut impl Rng) -> Vec<usize> {
    // Build position lookup: pos[node] = index in its route
    let mut pos = vec![0usize; inst.n + 1];
    let mut node_route: Vec<u16> = vec![0; inst.n + 1];
    for (ri, route) in sol.routes.iter().enumerate() {
        for (i, &v) in route.iter().enumerate() {
            pos[v] = i;
            node_route[v] = ri as u16;
        }
    }

    // Compute detour cost for each assigned pair
    let mut costs: Vec<(f64, usize)> = Vec::new();
    for route in &sol.routes {
        for &v in &route[1..route.len() - 1] {
            if v == 0 || !inst.is_pickup(v) {
                continue;
            }
            let d = inst.delivery_of(v);
            let ri = node_route[v] as usize;

            // Pickup detour
            let prev_p = sol.routes[ri][pos[v] - 1];
            let next_p = sol.routes[ri][pos[v] + 1];
            let cost_p = inst.dist(prev_p, v) + inst.dist(v, next_p) - inst.dist(prev_p, next_p);

            // Delivery detour (skip adjustment if adjacent — minor approximation)
            let prev_d = sol.routes[ri][pos[d] - 1];
            let next_d = sol.routes[ri][pos[d] + 1];
            let cost_d = inst.dist(prev_d, d) + inst.dist(d, next_d) - inst.dist(prev_d, next_d);

            costs.push((cost_p + cost_d, v));
        }
    }

    if costs.is_empty() {
        return Vec::new();
    }

    // Select top-`count` by cost (descending) with randomization
    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];
    let p_rand: f64 = 3.0;

    while removed.len() < count && !costs.is_empty() {
        #[allow(
            clippy::cast_possible_truncation,
            clippy::cast_sign_loss,
            clippy::cast_precision_loss
        )]
        let idx = (rng.random::<f64>().powf(p_rand) * costs.len() as f64) as usize;
        let idx = idx.min(costs.len() - 1);

        // Partial sort: bring the idx-th largest to position idx
        costs.select_nth_unstable_by(idx, |a, b| b.0.partial_cmp(&a.0).unwrap());
        let pickup = costs[idx].1;

        if !in_removed[pickup] {
            removed.push(pickup);
            in_removed[pickup] = true;
        }
        costs.swap_remove(idx);
    }

    removed
}

/// Late Arrival Removal (Gasque & Munari 2025): remove pairs whose actual
/// service times deviate most from the start of their time windows.
/// Score = (arrival_p - early_p) + (arrival_d - early_d).
/// Targets pairs that waste the most TW slack, making them good candidates
/// for reinsertion at better positions.
fn late_arrival_removal(
    inst: &Instance,
    sol: &Solution,
    count: usize,
    rng: &mut impl Rng,
) -> Vec<usize> {
    use crate::solution::fmax;

    // Compute arrival times and score each pair
    let mut costs: Vec<(f64, usize)> = Vec::new();
    for route in &sol.routes {
        let n = route.len();
        if n < 3 {
            continue;
        }

        // Forward pass: compute arrival times
        let mut arrival = vec![0.0f64; n];
        arrival[0] = inst.early(0);
        for i in 1..n {
            arrival[i] = fmax(
                arrival[i - 1] + inst.svc(route[i - 1]) + inst.dist(route[i - 1], route[i]),
                inst.early(route[i]),
            );
        }

        // Collect pickup positions and compute LA score
        let mut pos_of = vec![0usize; inst.n + 1];
        for (i, &v) in route.iter().enumerate() {
            pos_of[v] = i;
        }

        for i in 1..n - 1 {
            let v = route[i];
            if !inst.is_pickup(v) {
                continue;
            }
            let d = inst.delivery_of(v);
            let slack_p = arrival[i] - inst.early(v);
            let slack_d = arrival[pos_of[d]] - inst.early(d);
            costs.push((slack_p + slack_d, v));
        }
    }

    if costs.is_empty() {
        return Vec::new();
    }

    // Select top-`count` by score (descending) with randomization (same as worst_removal)
    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];
    let p_rand: f64 = 3.0;

    while removed.len() < count && !costs.is_empty() {
        #[allow(
            clippy::cast_possible_truncation,
            clippy::cast_sign_loss,
            clippy::cast_precision_loss
        )]
        let idx = (rng.random::<f64>().powf(p_rand) * costs.len() as f64) as usize;
        let idx = idx.min(costs.len() - 1);

        costs.select_nth_unstable_by(idx, |a, b| b.0.partial_cmp(&a.0).unwrap());
        let pickup = costs[idx].1;

        if !in_removed[pickup] {
            removed.push(pickup);
            in_removed[pickup] = true;
        }
        costs.swap_remove(idx);
    }

    removed
}

/// Detour-misplacement removal: for each pair, compute the detour cost in its
/// current route and compare it against the best feasible insertion cost across
/// all other routes. The "misplacement score" = current_detour - best_alternative.
/// Pairs with the highest scores are poorly placed and would benefit most from
/// relocation. Falls back to detour-only scoring when no alternative exists.
fn detour_misplacement_removal(
    inst: &Instance,
    sol: &Solution,
    count: usize,
    rng: &mut impl Rng,
) -> Vec<usize> {
    let mut pos = vec![0usize; inst.n + 1];
    let mut node_route: Vec<u16> = vec![0; inst.n + 1];
    for (ri, route) in sol.routes.iter().enumerate() {
        for (i, &v) in route.iter().enumerate() {
            pos[v] = i;
            node_route[v] = ri as u16;
        }
    }

    // For each pair, compute current detour and best alternative insertion cost.
    let mut scores: Vec<(f64, usize)> = Vec::new();
    for route in &sol.routes {
        for &v in &route[1..route.len() - 1] {
            if v == 0 || !inst.is_pickup(v) {
                continue;
            }
            let d = inst.delivery_of(v);
            let ri = node_route[v] as usize;

            // Current detour of pickup
            let prev_p = sol.routes[ri][pos[v] - 1];
            let next_p = sol.routes[ri][pos[v] + 1];
            let detour_p = inst.dist(prev_p, v) + inst.dist(v, next_p) - inst.dist(prev_p, next_p);

            // Current detour of delivery
            let prev_d = sol.routes[ri][pos[d] - 1];
            let next_d = sol.routes[ri][pos[d] + 1];
            let detour_d = inst.dist(prev_d, d) + inst.dist(d, next_d) - inst.dist(prev_d, next_d);

            let current_detour = detour_p + detour_d;

            // Estimate best insertion cost in other routes (cheapest arc pair).
            // Only consider a sample of routes to keep cost bounded.
            let mut best_alt = f64::MAX;
            let max_routes_to_check = sol.routes.len().min(8);
            for (rj, other_route) in sol.routes.iter().enumerate().take(max_routes_to_check) {
                if rj == ri || other_route.len() < 2 {
                    continue;
                }
                // Find cheapest insertion edge for pickup in this route
                let mut best_p = f64::MAX;
                for w in other_route.windows(2) {
                    let cost = inst.dist(w[0], v) + inst.dist(v, w[1]) - inst.dist(w[0], w[1]);
                    if cost < best_p {
                        best_p = cost;
                    }
                }
                // Find cheapest insertion edge for delivery in this route
                let mut best_d_in_route = f64::MAX;
                for w in other_route.windows(2) {
                    let cost = inst.dist(w[0], d) + inst.dist(d, w[1]) - inst.dist(w[0], w[1]);
                    if cost < best_d_in_route {
                        best_d_in_route = cost;
                    }
                }
                let alt_cost = best_p + best_d_in_route;
                if alt_cost < best_alt {
                    best_alt = alt_cost;
                }
            }

            // Misplacement score: how much worse current placement is vs best alternative.
            // If no alternative found, just use current detour as score.
            let score = if best_alt < f64::MAX {
                current_detour - best_alt
            } else {
                current_detour
            };
            scores.push((score, v));
        }
    }

    if scores.is_empty() {
        return Vec::new();
    }

    // Select top-`count` by score (descending) with randomization
    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];
    let p_rand: f64 = 3.0;

    while removed.len() < count && !scores.is_empty() {
        #[allow(
            clippy::cast_possible_truncation,
            clippy::cast_sign_loss,
            clippy::cast_precision_loss
        )]
        let idx = (rng.random::<f64>().powf(p_rand) * scores.len() as f64) as usize;
        let idx = idx.min(scores.len() - 1);

        scores.select_nth_unstable_by(idx, |a, b| b.0.partial_cmp(&a.0).unwrap());
        let pickup = scores[idx].1;

        if !in_removed[pickup] {
            removed.push(pickup);
            in_removed[pickup] = true;
        }
        scores.swap_remove(idx);
    }

    removed
}

/// Reduced-cost guided removal: score each pair by the sum of LP reduced costs
/// on arcs incident to its pickup and delivery in the current solution.
/// Higher reduced costs indicate the pair is placed on "expensive" arcs from
/// the LP relaxation's perspective — promising removal candidates.
fn rc_removal(inst: &Instance, sol: &Solution, count: usize, rng: &mut impl Rng) -> Vec<usize> {
    let mut pos = vec![0usize; inst.n + 1];
    let mut node_route: Vec<u16> = vec![0; inst.n + 1];
    for (ri, route) in sol.routes.iter().enumerate() {
        for (i, &v) in route.iter().enumerate() {
            pos[v] = i;
            node_route[v] = ri as u16;
        }
    }

    let mut costs: Vec<(f64, usize)> = Vec::new();
    for route in &sol.routes {
        for &v in &route[1..route.len() - 1] {
            if v == 0 || !inst.is_pickup(v) {
                continue;
            }
            let d = inst.delivery_of(v);
            let ri = node_route[v] as usize;

            // Sum reduced costs on arcs adjacent to pickup and delivery
            let prev_p = sol.routes[ri][pos[v] - 1];
            let next_p = sol.routes[ri][pos[v] + 1];
            let rc_p = inst.rc(prev_p, v) + inst.rc(v, next_p);

            let prev_d = sol.routes[ri][pos[d] - 1];
            let next_d = sol.routes[ri][pos[d] + 1];
            let rc_d = inst.rc(prev_d, d) + inst.rc(d, next_d);

            // Clamp f64::MAX contributions (infeasible arcs shouldn't dominate)
            let rc_total = rc_p.min(1e12) + rc_d.min(1e12);
            costs.push((rc_total, v));
        }
    }

    if costs.is_empty() {
        return Vec::new();
    }

    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];
    let p_rand: f64 = 3.0;

    while removed.len() < count && !costs.is_empty() {
        #[allow(
            clippy::cast_possible_truncation,
            clippy::cast_sign_loss,
            clippy::cast_precision_loss
        )]
        let idx = (rng.random::<f64>().powf(p_rand) * costs.len() as f64) as usize;
        let idx = idx.min(costs.len() - 1);

        costs.select_nth_unstable_by(idx, |a, b| b.0.partial_cmp(&a.0).unwrap());
        let pickup = costs[idx].1;

        if !in_removed[pickup] {
            removed.push(pickup);
            in_removed[pickup] = true;
        }
        costs.swap_remove(idx);
    }

    removed
}

/// SISRS Adjacent String Removal (Christiaens & Vanden Berghe, Transportation Science, 2020).
/// Removes adjacent strings of consecutive customers from routes, inducing spatial slack.
/// Algorithm 2 from the paper, adapted for PDPTW (pickup-delivery pairs).
#[allow(
    clippy::cast_possible_truncation,
    clippy::cast_sign_loss,
    clippy::cast_precision_loss
)]
fn string_removal(inst: &Instance, sol: &Solution, count: usize, rng: &mut impl Rng) -> Vec<usize> {
    if sol.routes.is_empty() {
        return Vec::new();
    }

    // Average tour cardinality |t ∈ T| (customers per route, excluding depot)
    let total_customers: usize = sol.routes.iter().map(|r| r.len().saturating_sub(2)).sum();
    let num_routes = sol.routes.len();
    let avg_tour_len = if num_routes > 0 {
        total_customers as f64 / num_routes as f64
    } else {
        1.0
    };

    // Equations (5)-(7): compute max string cardinality and number of strings
    let ls_max = (SISRS_L_MAX as f64).min(avg_tour_len); // Eq. (5)
    let ks_max = ((4.0 * count as f64) / (1.0 + ls_max) - 1.0).max(1.0); // Eq. (6)
    let ks: usize = (rng.random::<f64>() * ks_max).floor().max(1.0) as usize; // Eq. (7)

    // Build adjacency list for each customer sorted by distance (precompute)
    // For PDPTW: use all non-depot nodes
    let all_nodes: Vec<usize> = (1..=inst.n).collect();

    // Pick random seed customer
    let seed_idx = rng.random_range(0..all_nodes.len());
    let c_seed = all_nodes[seed_idx];

    // Build sorted adjacency list for seed (by distance)
    let mut adj: Vec<(f64, usize)> = all_nodes
        .iter()
        .filter(|&&v| v != c_seed)
        .map(|&v| (inst.dist(c_seed, v), v))
        .collect();
    adj.sort_unstable_by(|a, b| a.0.partial_cmp(&b.0).unwrap());

    // Track which routes have been ruined
    let mut ruined_routes: Vec<bool> = vec![false; num_routes];
    // Map customer -> route index
    let mut node_route: Vec<usize> = vec![usize::MAX; inst.n + 1];
    for (ri, route) in sol.routes.iter().enumerate() {
        for &v in route {
            if v != 0 {
                node_route[v] = ri;
            }
        }
    }

    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];
    let mut strings_removed: usize = 0;

    // Iterate over adjacency list of seed
    for &(_, c) in &adj {
        if strings_removed >= ks || removed.len() >= count {
            break;
        }
        if in_removed[c] {
            continue;
        }
        let ri = node_route[c];
        if ri == usize::MAX || ruined_routes[ri] {
            continue;
        }

        let route = &sol.routes[ri];
        let route_cust_count = route.len() - 2; // excluding depots
        if route_cust_count == 0 {
            continue;
        }

        // Find position of c* (closest non-removed to seed) in this route
        // c* is c itself since we checked !in_removed[c]
        let c_star = c;
        let c_star_pos = route.iter().position(|&v| v == c_star);
        let Some(c_star_pos) = c_star_pos else {
            continue;
        };

        // Eq. (8)-(9): compute string cardinality
        let lt_max = ls_max.min(route_cust_count as f64);
        let lt: usize = (rng.random::<f64>() * lt_max).floor().max(1.0) as usize;
        let lt = lt.min(route_cust_count);

        // Use split string procedure with probability SISRS_ALPHA
        let use_split = rng.random::<f64>() < SISRS_ALPHA;

        if use_split {
            // Split string: remove string of length l+m centered on c_star,
            // preserving m intervening customers
            let l = lt;
            // m starts at 1, incremented with probability (1-SISRS_BETA)
            let m_max = route_cust_count.saturating_sub(l);
            let mut m: usize = 1;
            while m < m_max {
                if rng.random::<f64>() < SISRS_BETA {
                    break;
                }
                m += 1;
            }

            // Select a string of cardinality l+m that includes c_star_pos
            // Then remove l customers, keeping m consecutive preserved
            let total_len = l + m;
            let total_len = total_len.min(route_cust_count);
            // Start position: ensure c_star is included and within depot bounds [1..n-1]
            let min_start = if c_star_pos >= total_len {
                c_star_pos - total_len + 1
            } else {
                1
            };
            let max_start = (route.len() - 1 - total_len).max(1);
            let start = rng.random_range(min_start.max(1)..=max_start.max(min_start.max(1)));

            // Within [start..start+total_len), remove customers but preserve m consecutive ones
            // Choose which m positions to preserve (consecutive block)
            let preserve_start = if total_len > m {
                rng.random_range(0..=(total_len - m))
            } else {
                0
            };

            for j in 0..total_len {
                let pos = start + j;
                if pos >= route.len() - 1 {
                    break;
                }
                if j >= preserve_start && j < preserve_start + m {
                    continue; // preserved
                }
                let v = route[pos];
                if v != 0 && !in_removed[v] {
                    // For PDPTW: when removing a node, also mark its pair
                    let pickup = if inst.is_pickup(v) {
                        v
                    } else {
                        inst.pickup_of(v)
                    };
                    if !in_removed[pickup] {
                        in_removed[pickup] = true;
                        in_removed[inst.delivery_of(pickup)] = true;
                        removed.push(pickup);
                    }
                }
            }
        } else {
            // Regular string removal: remove lt consecutive customers around c_star
            // Pick a random contiguous substring of length lt that includes c_star_pos
            let min_start = if c_star_pos >= lt {
                c_star_pos - lt + 1
            } else {
                1
            };
            let max_start = (route.len() - 1 - lt).max(1);
            let start = rng.random_range(min_start.max(1)..=max_start.max(min_start.max(1)));

            for j in 0..lt {
                let pos = start + j;
                if pos >= route.len() - 1 {
                    break;
                }
                let v = route[pos];
                if v != 0 && !in_removed[v] {
                    let pickup = if inst.is_pickup(v) {
                        v
                    } else {
                        inst.pickup_of(v)
                    };
                    if !in_removed[pickup] {
                        in_removed[pickup] = true;
                        in_removed[inst.delivery_of(pickup)] = true;
                        removed.push(pickup);
                    }
                }
            }
        }

        ruined_routes[ri] = true;
        strings_removed += 1;
    }

    // If we haven't removed enough, fall back to random additional removal
    if removed.len() < count {
        let mut assigned: Vec<usize> = Vec::new();
        for route in &sol.routes {
            for &v in route {
                if v != 0 && inst.is_pickup(v) && !in_removed[v] {
                    assigned.push(v);
                }
            }
        }
        assigned.shuffle(rng);
        for p in assigned {
            if removed.len() >= count {
                break;
            }
            removed.push(p);
        }
    }

    removed.truncate(count);
    removed
}

/// Regret insertion heuristic: insert unassigned pairs using regret-k.
/// Section 3.3: k randomly selected from {2, 3, 4, 5, #R}.
/// `noise_level`: fraction of max insertion cost used as noise amplitude (0 = deterministic).
///
/// Uses a cached insertion matrix: only recomputes the column for the route
/// that changed, saving ~25x `best_pair_insertion` calls.
#[allow(clippy::similar_names, clippy::float_cmp)]
/// Regret-k insertion heuristic.
///
/// When `route_cap > 0`, no new routes beyond `route_cap` are created.
/// Pairs that cannot be feasibly inserted stay in `sol.unassigned`.
pub fn regret_insertion(
    inst: &Instance,
    sol: &mut Solution,
    k: usize,
    noise_level: f64,
    sparse: bool,
    rng: &mut impl Rng,
    infos: &mut Vec<RouteInfo>,
    route_cap: usize,
    insertion_difficulty: bool,
) -> f64 {
    compute_all_infos(inst, &sol.routes, infos);

    let nu = sol.unassigned.len();
    if nu == 0 {
        let nr = sol.routes.len();
        return infos.iter().take(nr).map(|info| info.distance).sum();
    }
    let mut nr = sol.routes.len();

    // Cache matrix: for each pair pi × route ri, store best insertion.
    // Layout: cache[pi * max_routes + ri], 16 bytes per entry (vs 32 for Option<(f64, usize, usize)>).
    // Sentinel: delta == f64::INFINITY means infeasible (replaces None).
    const INF: f64 = f64::INFINITY;
    let max_routes = nr + nu; // worst case: one new route per pair
    let mut cache: Vec<(f64, u32, u32)> = vec![(INF, 0, 0); nu * max_routes];
    // Lazy validity via route epochs: entry is valid iff cache_epoch == route_epoch[ri].
    let mut cache_epoch: Vec<u32> = vec![u32::MAX; nu * max_routes];
    let mut route_epoch: Vec<u32> = vec![0; max_routes];
    // Track which pairs are still unassigned (by original index)
    let mut pair_alive: Vec<bool> = vec![true; nu];
    // Reusable buffer for route construction (avoids per-insertion heap allocation)
    let max_route_len = sol.routes.iter().map(Vec::len).max().unwrap_or(4) + 2;
    let mut route_buf: Vec<usize> = Vec::with_capacity(max_route_len);

    // Pre-compute per-route request bitmasks for O(stride) conflict checks
    let stride = inst.conflict_stride;
    let mut route_masks: Vec<u64> = vec![0u64; max_routes * stride];
    let mut mask_buf: Vec<u64> = Vec::with_capacity(stride);
    for ri in 0..nr {
        inst.build_route_req_mask(&sol.routes[ri], &mut mask_buf);
        route_masks[ri * stride..(ri * stride + stride)].copy_from_slice(&mask_buf);
    }

    // Initial fill: compute all (pair, route) entries
    // Safety: pi < nu, ri < nr, cache.len() == nu * max_routes, nr <= max_routes.
    // So pi * max_routes + ri < nu * max_routes == cache.len().
    for pi in 0..nu {
        let pickup = sol.unassigned[pi];
        let delivery = inst.delivery_of(pickup);
        for ri in 0..nr {
            let idx = pi * max_routes + ri;
            // Skip expensive insertion if pair conflicts with any request in route
            if inst
                .has_conflict_with_mask(&route_masks[ri * stride..(ri * stride + stride)], pickup)
            {
                unsafe {
                    *cache.get_unchecked_mut(idx) = (INF, 0, 0);
                    *cache_epoch.get_unchecked_mut(idx) = route_epoch[ri];
                }
                continue;
            }
            let result = if sparse {
                best_pair_insertion_with_info::<true>(
                    inst,
                    &sol.routes[ri],
                    &infos[ri],
                    pickup,
                    delivery,
                )
            } else {
                best_pair_insertion_with_info::<false>(
                    inst,
                    &sol.routes[ri],
                    &infos[ri],
                    pickup,
                    delivery,
                )
            };
            unsafe {
                *cache.get_unchecked_mut(idx) = match result {
                    Some((dist, ep, ed)) => (dist - infos[ri].distance, ep as u32, ed as u32),
                    None => (INF, 0, 0),
                };
                *cache_epoch.get_unchecked_mut(idx) = route_epoch[ri];
            }
        }
    }

    // Top-k buffer for regret computation (avoids sort — O(nr) instead of O(nr log nr))
    const TOP_K_CAP: usize = 6;

    let mut remaining = nu;
    while remaining > 0 {
        let mut best_regret_pi: Option<usize> = None;
        let mut best_regret_value: f64 = f64::MIN;
        let mut best_route_idx: Option<usize> = None;
        let mut best_ep: u32 = 0;
        let mut best_ed: u32 = 0;
        let mut best_cost: f64 = f64::MAX;

        // Effective k for top-k: cap at TOP_K_CAP for the sorted tracker
        let ek = k.min(TOP_K_CAP);

        for (pi, alive) in pair_alive.iter().enumerate() {
            if !alive {
                continue;
            }

            // Lazily recompute any invalid entries for this pair
            let pickup = sol.unassigned[pi];
            let delivery = inst.delivery_of(pickup);
            for ri in 0..nr {
                let idx = pi * max_routes + ri;
                unsafe {
                    if *cache_epoch.get_unchecked(idx) != route_epoch[ri] {
                        // Skip expensive insertion if pair conflicts with route
                        let result = if inst.has_conflict_with_mask(
                            &route_masks[ri * stride..(ri * stride + stride)],
                            pickup,
                        ) {
                            None
                        } else if sparse {
                            best_pair_insertion_with_info::<true>(
                                inst,
                                &sol.routes[ri],
                                &infos[ri],
                                pickup,
                                delivery,
                            )
                        } else {
                            best_pair_insertion_with_info::<false>(
                                inst,
                                &sol.routes[ri],
                                &infos[ri],
                                pickup,
                                delivery,
                            )
                        };
                        *cache.get_unchecked_mut(idx) = match result {
                            Some((dist, ep, ed)) => {
                                (dist - infos[ri].distance, ep as u32, ed as u32)
                            }
                            None => (INF, 0, 0),
                        };
                        *cache_epoch.get_unchecked_mut(idx) = route_epoch[ri];
                    }
                }
            }

            // Track the best (minimum-cost) insertion and the sum of the k
            // smallest noisy deltas for regret computation.  Uses an unsorted
            // buffer with explicit max-tracking to avoid insertion-sort swaps.
            let mut min_delta: f64 = f64::MAX;
            let mut min_ri: usize = 0;
            let mut min_ep: u32 = 0;
            let mut min_ed: u32 = 0;

            let mut top_costs: [f64; TOP_K_CAP] = [f64::MAX; TOP_K_CAP];
            let mut top_count: usize = 0;
            let mut top_sum: f64 = 0.0;
            let mut top_max_val: f64 = f64::MIN;
            let mut top_max_idx: usize = 0;

            let mut feasible_count: usize = 0;
            let mut total_sum: f64 = 0.0;

            // Compute max delta for noise scaling (only when noise is active)
            let max_delta = if noise_level > 0.0 {
                let mut md: f64 = 0.0;
                for ri in 0..nr {
                    let d = unsafe { cache.get_unchecked(pi * max_routes + ri).0 };
                    if d < INF && d > md {
                        md = d;
                    }
                }
                md
            } else {
                0.0
            };
            let noise_amp = noise_level * max_delta.max(0.0);

            for ri in 0..nr {
                let (delta, ep, ed) = unsafe { *cache.get_unchecked(pi * max_routes + ri) };
                if delta < INF {
                    feasible_count += 1;
                    total_sum += delta;

                    // Add noise to delta for top-k ranking (Ropke & Pisinger 2006)
                    let noisy_delta = if noise_amp > 0.0 {
                        delta + (rng.random::<f64>() * 2.0 - 1.0) * noise_amp
                    } else {
                        delta
                    };

                    // Track overall minimum for insertion
                    if noisy_delta < min_delta {
                        min_delta = noisy_delta;
                        min_ri = ri;
                        min_ep = ep;
                        min_ed = ed;
                    }

                    // Maintain unsorted top-k buffer with max-tracking
                    if top_count < ek {
                        top_costs[top_count] = noisy_delta;
                        top_sum += noisy_delta;
                        if top_count == 0 || noisy_delta > top_max_val {
                            top_max_val = noisy_delta;
                            top_max_idx = top_count;
                        }
                        top_count += 1;
                    } else if noisy_delta < top_max_val {
                        top_sum = top_sum - top_max_val + noisy_delta;
                        top_costs[top_max_idx] = noisy_delta;
                        // Rescan for new max
                        top_max_val = f64::MIN;
                        for (j, &cost) in top_costs.iter().enumerate().take(ek) {
                            if cost > top_max_val {
                                top_max_val = cost;
                                top_max_idx = j;
                            }
                        }
                    }
                }
            }

            if feasible_count == 0 {
                continue;
            }

            #[allow(clippy::cast_precision_loss)]
            let regret = if insertion_difficulty {
                // Insertion Difficulty (Definition 4, Kim et al.):
                // ID = Δ_2nd / Δ_best  if Δ_best > 0  (≥1, high when 2nd much worse)
                // ID = 0               if Δ_best ≈ 0
                // ID = MAX             if Δ_best < 0, Δ_2nd ≥ 0  (best saves, 2nd costs)
                // ID = Δ_best / Δ_2nd  if Δ_best < 0, Δ_2nd < 0  (≥1, more neg / less neg)
                // Higher ID → harder to insert → insert first.
                if feasible_count == 1 {
                    f64::MAX
                } else {
                    // Find 2nd-best delta from top-k buffer
                    let second_best = if top_count >= 2 {
                        let mut s = f64::MAX;
                        for &cost in &top_costs[..top_count] {
                            if cost > min_delta && cost < s {
                                s = cost;
                            }
                        }
                        // If all top entries equal min_delta, scan remaining routes
                        if s == f64::MAX { min_delta } else { s }
                    } else {
                        min_delta
                    };
                    if min_delta > 0.0 {
                        if min_delta < 1e-12 {
                            f64::MAX
                        } else {
                            second_best / min_delta
                        }
                    } else if min_delta.abs() < 1e-12 {
                        0.0
                    } else {
                        // min_delta < 0: beneficial insertion
                        if second_best > -1e-12 {
                            f64::MAX
                        } else {
                            min_delta / second_best
                        }
                    }
                }
            } else if feasible_count == 1 {
                f64::MAX
            } else if k >= feasible_count {
                // All feasible routes contribute to regret
                total_sum - feasible_count as f64 * min_delta
            } else {
                // regret = sum_{h=1..top_k}(top[h] - top[0])
                //        = top_sum - top_count * min_delta
                top_sum - top_count as f64 * min_delta
            };

            if regret > best_regret_value || (regret == best_regret_value && min_delta < best_cost)
            {
                best_regret_value = regret;
                best_regret_pi = Some(pi);
                best_route_idx = Some(min_ri);
                best_ep = min_ep;
                best_ed = min_ed;
                best_cost = min_delta;
            }
        }

        if let Some(pi) = best_regret_pi {
            let pickup = sol.unassigned[pi];
            let delivery = inst.delivery_of(pickup);
            let ri = best_route_idx.unwrap();
            // Fused: build route with pair + compute light info in one pass
            // (saves one O(n) traversal vs. separate build + compute)
            compute_info_with_pair(
                inst,
                &sol.routes[ri],
                pickup,
                best_ep as usize,
                delivery,
                best_ed as usize,
                &mut route_buf,
                &mut infos[ri],
            );
            std::mem::swap(&mut sol.routes[ri], &mut route_buf);
            // Upgrade to full info (cap_excess/tw_excess needed by penalty-based caller)
            infos[ri].compute_tw_sector(inst, &sol.routes[ri], sol.routes[ri].len());

            pair_alive[pi] = false;
            remaining -= 1;

            // O(1) incremental mask update instead of O(n) rebuild
            let req = inst.pickup_to_req[pickup];
            if req != usize::MAX {
                route_masks[ri * stride + (req >> 6)] |= 1u64 << (req & 63);
            }

            // Invalidate cache column ri in O(1): bump route epoch.
            route_epoch[ri] = route_epoch[ri].wrapping_add(1);
        } else if route_cap > 0 && nr >= route_cap {
            // Route cap active and reached — leave remaining pairs unassigned.
            break;
        } else {
            // No feasible insertion found — create a new route for the first alive pair
            let pi = (0..nu).find(|&i| pair_alive[i]).unwrap();
            let pickup = sol.unassigned[pi];
            let delivery = inst.delivery_of(pickup);
            sol.routes.push(vec![0, pickup, delivery, 0]);

            // Grow infos and compute new route
            while infos.len() < sol.routes.len() {
                infos.push(RouteInfo::new());
            }
            let new_ri = sol.routes.len() - 1;
            infos[new_ri].compute(inst, &sol.routes[new_ri]);
            nr = sol.routes.len();

            // Grow and compute mask for the new route
            route_masks.resize(nr * stride, 0u64);
            inst.build_route_req_mask(&sol.routes[new_ri], &mut mask_buf);
            route_masks[new_ri * stride..(new_ri * stride + stride)].copy_from_slice(&mask_buf);

            pair_alive[pi] = false;
            remaining -= 1;

            // New route cache entries are already stale (epoch mismatch from u32::MAX).
        }
    }

    // Update sol.unassigned: keep only pairs that were not inserted.
    let old_unassigned = std::mem::take(&mut sol.unassigned);
    sol.unassigned = old_unassigned
        .into_iter()
        .enumerate()
        .filter(|&(pi, _)| pair_alive[pi])
        .map(|(_, p)| p)
        .collect();

    // Return total distance from cached infos (avoids recomputing in caller)
    let nr = sol.routes.len();
    infos.iter().take(nr).map(|info| info.distance).sum()
}

/// Route elimination removal: remove ALL pairs from the shortest route,
/// then fill up to `count` with Shaw-similar pairs from neighboring routes.
/// Gives LNS a chance to eliminate a vehicle by absorbing a short route.
#[allow(
    clippy::cast_possible_truncation,
    clippy::cast_sign_loss,
    clippy::cast_precision_loss
)]
fn route_elimination_removal(
    inst: &Instance,
    sol: &Solution,
    count: usize,
    rng: &mut impl Rng,
) -> Vec<usize> {
    if sol.routes.is_empty() {
        return Vec::new();
    }

    // Find shortest route by customer count
    let min_idx = sol
        .routes
        .iter()
        .enumerate()
        .min_by_key(|(_, r)| r.len())
        .unwrap()
        .0;

    let mut removed = Vec::new();
    let mut in_removed = vec![false; inst.n + 1];

    // Remove ALL pickups from shortest route
    for &v in &sol.routes[min_idx] {
        if v != 0 && inst.is_pickup(v) {
            removed.push(v);
            in_removed[v] = true;
        }
    }

    // Fill up to count with Shaw-similar pairs from other routes
    if removed.len() < count {
        let phi1: f64 = 9.0;
        let phi2: f64 = 3.0;
        let phi3: f64 = 2.0;

        // Compute arrival times for Shaw similarity
        let mut arrival_time: Vec<f64> = vec![0.0; inst.n + 1];
        for route in &sol.routes {
            let mut time = inst.early(0);
            for i in 1..route.len() {
                let prev = route[i - 1];
                let curr = route[i];
                time = crate::solution::fmax(
                    time + inst.svc(prev) + inst.dist(prev, curr),
                    inst.early(curr),
                );
                arrival_time[curr] = time;
            }
        }

        // Collect non-removed assigned pickups
        let mut assigned: Vec<usize> = Vec::new();
        for route in &sol.routes {
            for &v in route {
                if v != 0 && inst.is_pickup(v) && !in_removed[v] {
                    assigned.push(v);
                }
            }
        }

        let p_rand: f64 = 6.0;
        let mut similarities: Vec<(f64, usize)> = Vec::with_capacity(assigned.len());
        while removed.len() < count && !assigned.is_empty() {
            let ref_pickup = removed[rng.random_range(0..removed.len())];

            similarities.clear();
            similarities.extend(assigned.iter().filter(|&&p| !in_removed[p]).map(|&p| {
                let sim = phi1 * inst.dist(ref_pickup, p)
                    + phi2 * (arrival_time[ref_pickup] - arrival_time[p]).abs()
                    + phi3 * f64::from((inst.demand[ref_pickup] - inst.demand[p]).abs());
                (sim, p)
            }));

            if similarities.is_empty() {
                break;
            }

            let idx = (rng.random::<f64>().powf(p_rand) * similarities.len() as f64) as usize;
            let idx = idx.min(similarities.len() - 1);
            similarities.select_nth_unstable_by(idx, |a, b| a.0.partial_cmp(&b.0).unwrap());
            let selected = similarities[idx].1;

            removed.push(selected);
            in_removed[selected] = true;
        }
    }

    removed
}

#[cfg(test)]
mod tests {
    use super::{TabuList, branch_and_bound_repair, is_tabu_admissible, tabu_move_signature};
    use crate::instance::tests::make_test_instance;
    use crate::solution::Solution;
    use std::time::{Duration, Instant};

    #[test]
    fn tabu_signature_ignores_removed_order() {
        let a = tabu_move_signature(&[3, 7, 11]);
        let b = tabu_move_signature(&[11, 3, 7]);
        assert_eq!(a, b);
    }

    #[test]
    fn tabu_entry_expires_after_tenure() {
        let mut tabu = TabuList::new(3);
        let sig = tabu_move_signature(&[5, 9]);
        tabu.add_with_tenure(sig, 10, 3);

        assert!(tabu.contains(sig, 10));
        assert!(tabu.contains(sig, 12));
        assert!(!tabu.contains(sig, 13));
    }

    #[test]
    fn aspiration_overrides_tabu_for_global_improvement() {
        let mut tabu = TabuList::new(5);
        let sig = tabu_move_signature(&[1, 2]);
        tabu.add_with_tenure(sig, 20, 5);

        assert!(!is_tabu_admissible(&tabu, sig, 16, 100.0, 90.0));
        assert!(is_tabu_admissible(&tabu, sig, 16, 80.0, 90.0));
    }

    #[test]
    fn bnb_repair_inserts_removed_pairs_on_tiny_instance() {
        let inst = make_test_instance(
            2,
            2,
            &[
                0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 2.0, 1.0, 2.0, 1.0, 2.0, 0.0, 2.0, 1.0, 1.0,
                1.0, 2.0, 0.0, 2.0, 1.0, 2.0, 1.0, 2.0, 0.0,
            ],
            &[0, 1, 1, -1, -1],
            &[0.0, 0.0, 0.0, 0.0, 0.0],
            &[1000.0, 1000.0, 1000.0, 1000.0, 1000.0],
            &[0.0, 0.0, 0.0, 0.0, 0.0],
        );

        let mut sol = Solution {
            routes: vec![vec![0, 0]],
            unassigned: vec![1, 2],
        };

        let deadline = Instant::now() + Duration::from_millis(100);
        assert!(branch_and_bound_repair(&inst, &mut sol, deadline));
        assert!(sol.unassigned.is_empty());
        assert!(sol.is_feasible(&inst));
    }
}