pdp_lns 0.1.1

use crate::four_opt;
use crate::instance::Instance;
use crate::solution::{
    CircleSector, PenaltyState, RouteInfo, Solution, best_pair_insertion_penalized,
    best_pair_insertion_with_info, build_route_with_pair, build_route_with_pair_buf,
    compute_all_infos, compute_info_with_pair, compute_info_without_pair, route_without_pair,
    route_without_pair_buf, single_pair_feasible_distance,
};
use crate::special_opt;
use crate::two_k_opt;
use std::fmt::Write;
use std::time::Instant;

const OP_NAMES: [&str; 10] = ["M1", "M2", "M3", "M4", "M5", "M6", "M7", "M8", "M9", "M10"];

/// Compute a fast signature for a route. Used to detect which routes changed
/// between operations, enabling dirty route hints for local search.
#[inline]
fn route_sig(route: &[usize]) -> u64 {
    let mut h: u64 = route.len() as u64;
    for &v in route {
        h = h
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(v as u64);
    }
    h
}

/// Snapshot route signatures for a solution. Used with `dirty_route_indices`
/// to detect which routes changed after an operation.
pub fn snapshot_routes(routes: &[Vec<usize>]) -> Vec<u64> {
    routes.iter().map(|r| route_sig(r)).collect()
}

/// Compare current routes against a snapshot and return indices of changed routes.
/// Also includes any newly added routes (indices beyond the snapshot length).
pub fn dirty_route_indices(routes: &[Vec<usize>], snapshot: &[u64]) -> Vec<usize> {
    let mut dirty = Vec::new();
    for (i, r) in routes.iter().enumerate() {
        if i >= snapshot.len() || route_sig(r) != snapshot[i] {
            dirty.push(i);
        }
    }
    dirty
}

#[derive(Default, Clone, Debug)]
pub struct LsStats {
    pub moves: [u64; 10],
    pub delta_sum: [f64; 10],
    pub time_us: [u64; 10],
    pub scans: [u64; 10], // number of invocations (including 0-move scans)
    pub outer_loops: u64,
    pub total_time_us: u64,
    pub calls: u64, // number of local_search(_full) invocations
}

impl LsStats {
    #[allow(dead_code)]
    pub fn merge_from(&mut self, other: &Self) {
        for i in 0..10 {
            self.moves[i] += other.moves[i];
            self.delta_sum[i] += other.delta_sum[i];
            self.time_us[i] += other.time_us[i];
            self.scans[i] += other.scans[i];
        }
        self.outer_loops += other.outer_loops;
        self.total_time_us += other.total_time_us;
        self.calls += other.calls;
    }

    pub fn total_moves(&self) -> u64 {
        self.moves.iter().sum()
    }

    pub fn summary_line(&self) -> String {
        let total_ms = self.total_time_us / 1000;
        let mut parts = String::new();
        for (i, &name) in OP_NAMES.iter().enumerate() {
            let ms = self.time_us[i] / 1000;
            if self.moves[i] > 0 {
                if !parts.is_empty() {
                    parts.push(' ');
                }
                let _ = write!(
                    parts,
                    "{}:{}/{:.1}/{}ms",
                    name, self.moves[i], self.delta_sum[i], ms,
                );
            } else if ms > 0 {
                // Show operators that consumed time but found nothing
                if !parts.is_empty() {
                    parts.push(' ');
                }
                let _ = write!(parts, "{name}:0/{ms}ms");
            }
        }
        format!(
            "{} moves, {}ms, {} calls, {} loops | {}",
            self.total_moves(),
            total_ms,
            self.calls,
            self.outer_loops,
            parts,
        )
    }
}

/// Apply the full local search: cycle M1->M2->M3->M4 until no improvement.
/// Per the paper (Section 3.1):
/// - M1, M2: best improvement
/// - M3, M4: first improvement
/// - Each operator is applied exhaustively before moving to the next.
/// - When M4 is exhausted, return to M1. Repeat until no improvement from any operator.
///
/// If `full` is true, also applies 2-Opt, Or-Opt, and 2k-Opt (Pacheco et al. 2022)
/// for deeper intra-route optimization. These are slower and best used during
/// distance refinement rather than vehicle minimization phases.
pub fn local_search(
    inst: &Instance,
    sol: &mut Solution,
    penalty: &PenaltyState,
    stats: Option<&mut LsStats>,
) {
    local_search_inner(inst, sol, false, penalty, stats, None);
}

/// Full local search including expensive intra-route operators (2-Opt, Or-Opt, 2k-Opt).
pub fn local_search_full(
    inst: &Instance,
    sol: &mut Solution,
    penalty: &PenaltyState,
    stats: Option<&mut LsStats>,
) {
    local_search_inner(inst, sol, true, penalty, stats, None);
}

/// Local search with dirty route hints: only routes in `dirty_routes` are marked
/// as modified initially, allowing operators to skip pairs involving clean routes.
pub fn local_search_dirty(
    inst: &Instance,
    sol: &mut Solution,
    penalty: &PenaltyState,
    stats: Option<&mut LsStats>,
    dirty_routes: &[usize],
) {
    local_search_inner(inst, sol, false, penalty, stats, Some(dirty_routes));
}

/// Full local search with dirty route hints.
pub fn local_search_full_dirty(
    inst: &Instance,
    sol: &mut Solution,
    penalty: &PenaltyState,
    stats: Option<&mut LsStats>,
    dirty_routes: &[usize],
) {
    local_search_inner(inst, sol, true, penalty, stats, Some(dirty_routes));
}

/// Fast intra-route optimization only: 2-Opt + Or-Opt.
/// Does NOT move pairs between routes — only improves within each route.
/// Typical cost: 1-3ms for 24 routes, ~50 customers each.
pub fn local_search_intra(inst: &Instance, sol: &mut Solution) {
    // 2-opt and or-opt never READ from infos during search — they compute
    // deltas from inst.dist() directly.  They only WRITE to infos[ri] after
    // a successful move (via infos[ri].compute()).  So skip the expensive
    // upfront compute_all_infos; just ensure the Vec has the right length
    // with empty RouteInfos.  This saves ~125μs per call (full route-info
    // recomputation for all ~25 routes), which adds up over the hundreds
    // of calls from the LNS inner loop in distance mode.
    let nr = sol.routes.len();
    let mut infos: Vec<RouteInfo> = Vec::with_capacity(nr);
    infos.resize_with(nr, RouteInfo::new);
    loop {
        let mut any = false;
        while two_k_opt::apply_two_opt(inst, &mut sol.routes, &mut infos) {
            any = true;
        }
        while two_k_opt::apply_or_opt(inst, &mut sol.routes, &mut infos) {
            any = true;
        }
        if !any {
            break;
        }
    }
}

fn local_search_inner(
    inst: &Instance,
    sol: &mut Solution,
    full: bool,
    penalty: &PenaltyState,
    mut stats: Option<&mut LsStats>,
    dirty_routes: Option<&[usize]>,
) {
    let t_total = Instant::now();
    let tracking = stats.is_some();
    if let Some(s) = stats.as_deref_mut() {
        s.calls += 1;
    }

    // Reusable pool of RouteInfo — grows once, reused across all operators
    let mut infos: Vec<RouteInfo> = Vec::new();
    let mut info_buf = RouteInfo::new();
    let mut info_buf2 = RouteInfo::new();
    let mut dirty = true;

    // Route modification tracking for inter-route operator skip logic.
    // route_epoch[ri] = counter value when route ri was last modified.
    // Operators skip route pair (ri, rj) if both epochs <= their clean_epoch.
    let max_nr = sol.routes.len() + sol.unassigned.len() + 1;
    let mut move_ctr: u32 = 1;
    let mut route_epoch: Vec<u32>;
    let mut m2_clean: u32;
    let mut m3_clean: u32;
    let mut m4_clean: u32;
    let mut m5_clean: u32;
    let mut m6_clean: u32;
    let mut full_clean: u32;
    let mut m10_clean: u32;
    let mut m5_route_clean: Vec<u32>;
    let mut m6_route_clean: Vec<u32>;

    if let Some(dr) = dirty_routes {
        // Selective initialization: only specified routes are dirty.
        // Clean routes keep epoch=0, dirty routes get epoch=1.
        // Operator clean epochs start at 0, so operators will skip
        // pairs where both routes are clean.
        route_epoch = vec![0; max_nr];
        for &ri in dr {
            if ri < max_nr {
                route_epoch[ri] = 1;
            }
        }
        // If no routes are dirty, skip entirely
        if dr.is_empty() {
            if let Some(s) = stats.as_deref_mut() {
                s.total_time_us += t_total.elapsed().as_micros() as u64;
            }
            return;
        }
        m2_clean = 0;
        m3_clean = 0;
        m4_clean = 0;
        m5_clean = 0;
        m6_clean = 0;
        full_clean = 0;
        m10_clean = 0;
        // Pre-mark clean routes for M5/M6 per-route tracking:
        // clean routes get route_clean[ri] = 1 (matching their epoch of 0,
        // so route_epoch[ri] <= route_clean[ri] and they'll be skipped).
        m5_route_clean = vec![1; max_nr];
        m6_route_clean = vec![1; max_nr];
        for &ri in dr {
            if ri < max_nr {
                m5_route_clean[ri] = 0;
                m6_route_clean[ri] = 0;
            }
        }
    } else {
        // Default: all routes dirty (original behavior)
        route_epoch = vec![1; max_nr];
        m2_clean = 0;
        m3_clean = 0;
        m4_clean = 0;
        m5_clean = 0;
        m6_clean = 0;
        full_clean = 0;
        m10_clean = 0;
        m5_route_clean = vec![0; max_nr];
        m6_route_clean = vec![0; max_nr];
    }

    loop {
        if let Some(s) = stats.as_deref_mut() {
            s.outer_loops += 1;
        }
        let mut any_improved = false;

        // M1: Insert unassigned pairs
        {
            if let Some(s) = stats.as_deref_mut() {
                s.scans[0] += 1;
            }
            let t = Instant::now();
            let dist_before = if tracking {
                sol.total_distance(inst)
            } else {
                0.0
            };
            while apply_m1(
                inst,
                sol,
                &mut infos,
                &mut dirty,
                penalty,
                &mut move_ctr,
                &mut route_epoch,
            ) {
                if let Some(s) = stats.as_deref_mut() {
                    s.moves[0] += 1;
                }
                any_improved = true;
            }
            if let Some(s) = stats.as_deref_mut() {
                s.time_us[0] += t.elapsed().as_micros() as u64;
                s.delta_sum[0] += sol.total_distance(inst) - dist_before;
            }
        }

        // M2: Relocate assigned pair to different route
        // Skip if no routes changed since M2 last exhausted.
        if move_ctr > m2_clean {
            if let Some(s) = stats.as_deref_mut() {
                s.scans[1] += 1;
            }
            let t = Instant::now();
            let dist_before = if tracking {
                sol.total_distance(inst)
            } else {
                0.0
            };
            while apply_m2(
                inst,
                sol,
                &mut infos,
                &mut dirty,
                penalty,
                &mut move_ctr,
                &mut route_epoch,
            ) {
                if let Some(s) = stats.as_deref_mut() {
                    s.moves[1] += 1;
                }
                any_improved = true;
            }
            m2_clean = move_ctr;
            if let Some(s) = stats.as_deref_mut() {
                s.time_us[1] += t.elapsed().as_micros() as u64;
                s.delta_sum[1] += sol.total_distance(inst) - dist_before;
            }
        }

        // M3: Swap pairs between two routes
        // Skip if no routes changed since M3 last exhausted.
        if move_ctr > m3_clean {
            if let Some(s) = stats.as_deref_mut() {
                s.scans[2] += 1;
            }
            let t = Instant::now();
            let dist_before = if tracking {
                sol.total_distance(inst)
            } else {
                0.0
            };
            while apply_m3(
                inst,
                sol,
                &mut infos,
                &mut info_buf,
                &mut info_buf2,
                &mut dirty,
                &mut route_epoch,
                m3_clean,
                &mut move_ctr,
                penalty,
            ) {
                if let Some(s) = stats.as_deref_mut() {
                    s.moves[2] += 1;
                }
                any_improved = true;
            }
            m3_clean = move_ctr;
            if let Some(s) = stats.as_deref_mut() {
                s.time_us[2] += t.elapsed().as_micros() as u64;
                s.delta_sum[2] += sol.total_distance(inst) - dist_before;
            }
        }

        // M4: Transfer pairs between routes
        // Skip if no routes changed since M4 last exhausted.
        if move_ctr > m4_clean {
            if let Some(s) = stats.as_deref_mut() {
                s.scans[3] += 1;
            }
            let t = Instant::now();
            let dist_before = if tracking {
                sol.total_distance(inst)
            } else {
                0.0
            };
            while apply_m4(
                inst,
                sol,
                &mut infos,
                &mut info_buf,
                &mut dirty,
                penalty,
                &mut move_ctr,
                &mut route_epoch,
                m4_clean,
            ) {
                if let Some(s) = stats.as_deref_mut() {
                    s.moves[3] += 1;
                }
                any_improved = true;
            }
            m4_clean = move_ctr;
            if let Some(s) = stats.as_deref_mut() {
                s.time_us[3] += t.elapsed().as_micros() as u64;
                s.delta_sum[3] += sol.total_distance(inst) - dist_before;
            }
        }

        // M5: Intra-route pair relocation (improve positions within each route)
        // Skip if no routes changed since M5 last exhausted.
        if move_ctr > m5_clean {
            if let Some(s) = stats.as_deref_mut() {
                s.scans[4] += 1;
            }
            let t = Instant::now();
            let dist_before = if tracking {
                sol.total_distance(inst)
            } else {
                0.0
            };
            while apply_m5(
                inst,
                sol,
                &mut infos,
                &mut info_buf,
                &mut dirty,
                &mut move_ctr,
                &mut route_epoch,
                penalty,
                &mut m5_route_clean,
            ) {
                if let Some(s) = stats.as_deref_mut() {
                    s.moves[4] += 1;
                }
                any_improved = true;
            }
            m5_clean = move_ctr;
            if let Some(s) = stats.as_deref_mut() {
                s.time_us[4] += t.elapsed().as_micros() as u64;
                s.delta_sum[4] += sol.total_distance(inst) - dist_before;
            }
        }

        // M6: Intra-route 2-pair relocation (escape saddle points M5 can't)
        // Skip if no routes changed since M6 last exhausted.
        if move_ctr > m6_clean {
            if let Some(s) = stats.as_deref_mut() {
                s.scans[5] += 1;
            }
            let t = Instant::now();
            let dist_before = if tracking {
                sol.total_distance(inst)
            } else {
                0.0
            };
            while apply_m6(
                inst,
                sol,
                &mut infos,
                &mut info_buf,
                &mut info_buf2,
                &mut dirty,
                &mut move_ctr,
                &mut route_epoch,
                penalty,
                &mut m6_route_clean,
            ) {
                if let Some(s) = stats.as_deref_mut() {
                    s.moves[5] += 1;
                }
                any_improved = true;
            }
            m6_clean = move_ctr;
            if let Some(s) = stats.as_deref_mut() {
                s.time_us[5] += t.elapsed().as_micros() as u64;
                s.delta_sum[5] += sol.total_distance(inst) - dist_before;
            }
        }

        if full {
            // Ensure infos are up-to-date before M7-M10
            if dirty {
                compute_all_infos(inst, &sol.routes, &mut infos);
                dirty = false;
            }

            // M7-M9: intra-route optimization block.
            // Skip if no routes changed since this block last exhausted.
            if move_ctr > full_clean {
                // Snapshot distances before M7-M9 for targeted epoch updates
                let nr_full = sol.routes.len();
                let pre_full_dist: Vec<f64> =
                    infos.iter().take(nr_full).map(|i| i.distance).collect();

                // M7: Intra-route 2-Opt (Pacheco et al. 2022)
                {
                    if let Some(s) = stats.as_deref_mut() {
                        s.scans[6] += 1;
                    }
                    let t = Instant::now();
                    let dist_before = if tracking {
                        infos.iter().take(nr_full).map(|i| i.distance).sum::<f64>()
                    } else {
                        0.0
                    };
                    while two_k_opt::apply_two_opt(inst, &mut sol.routes, &mut infos) {
                        if let Some(s) = stats.as_deref_mut() {
                            s.moves[6] += 1;
                        }
                        any_improved = true;
                    }
                    if let Some(s) = stats.as_deref_mut() {
                        s.time_us[6] += t.elapsed().as_micros() as u64;
                        let dist_after: f64 = infos.iter().take(nr_full).map(|i| i.distance).sum();
                        s.delta_sum[6] += dist_after - dist_before;
                    }
                }

                // M8: Intra-route Or-Opt — segment relocation (Pacheco et al. 2022)
                {
                    if let Some(s) = stats.as_deref_mut() {
                        s.scans[7] += 1;
                    }
                    let t = Instant::now();
                    let dist_before = if tracking {
                        infos.iter().take(nr_full).map(|i| i.distance).sum::<f64>()
                    } else {
                        0.0
                    };
                    while two_k_opt::apply_or_opt(inst, &mut sol.routes, &mut infos) {
                        if let Some(s) = stats.as_deref_mut() {
                            s.moves[7] += 1;
                        }
                        any_improved = true;
                    }
                    if let Some(s) = stats.as_deref_mut() {
                        s.time_us[7] += t.elapsed().as_micros() as u64;
                        let dist_after: f64 = infos.iter().take(nr_full).map(|i| i.distance).sum();
                        s.delta_sum[7] += dist_after - dist_before;
                    }
                }

                // 4-Opt: intra-route 4-edge exchange (Pacheco et al. 2022)
                while four_opt::apply_four_opt(inst, &mut sol.routes, &mut infos) {
                    any_improved = true;
                }

                // M9: Inter-route segment exchange
                {
                    if let Some(s) = stats.as_deref_mut() {
                        s.scans[8] += 1;
                    }
                    let t = Instant::now();
                    let dist_before = if tracking {
                        infos.iter().take(nr_full).map(|i| i.distance).sum::<f64>()
                    } else {
                        0.0
                    };
                    while special_opt::apply_segment_exchange(inst, &mut sol.routes, &mut infos) {
                        if let Some(s) = stats.as_deref_mut() {
                            s.moves[8] += 1;
                        }
                        any_improved = true;
                    }
                    if let Some(s) = stats.as_deref_mut() {
                        s.time_us[8] += t.elapsed().as_micros() as u64;
                        let dist_after: f64 = infos.iter().take(nr_full).map(|i| i.distance).sum();
                        s.delta_sum[8] += dist_after - dist_before;
                    }
                }

                // Update epochs only for routes whose distance changed during M7-M9
                for ri in 0..nr_full {
                    if (infos[ri].distance - pre_full_dist[ri]).abs() > 1e-10 {
                        move_ctr += 1;
                        route_epoch[ri] = move_ctr;
                    }
                }
                full_clean = move_ctr;
            }

            // M10: SWAP* — inter-route pair exchange with free reinsertion (Vidal 2022)
            {
                if let Some(s) = stats.as_deref_mut() {
                    s.scans[9] += 1;
                }
                let t = Instant::now();
                let dist_before = if tracking {
                    sol.total_distance(inst)
                } else {
                    0.0
                };
                while apply_swap_star(
                    inst,
                    sol,
                    &mut infos,
                    &mut dirty,
                    &mut route_epoch,
                    m10_clean,
                    &mut move_ctr,
                    penalty,
                ) {
                    if let Some(s) = stats.as_deref_mut() {
                        s.moves[9] += 1;
                    }
                    any_improved = true;
                }
                m10_clean = move_ctr;
                if let Some(s) = stats.as_deref_mut() {
                    s.time_us[9] += t.elapsed().as_micros() as u64;
                    s.delta_sum[9] += sol.total_distance(inst) - dist_before;
                }
            }
        }

        if !any_improved {
            break;
        }
    }

    if let Some(s) = stats {
        s.total_time_us += t_total.elapsed().as_micros() as u64;
    }
}

/// Ensure infos are up-to-date. Skips recomputation if not dirty.
#[inline]
fn ensure_infos(
    inst: &Instance,
    routes: &[Vec<usize>],
    infos: &mut Vec<RouteInfo>,
    dirty: &mut bool,
) {
    if *dirty {
        compute_all_infos(inst, routes, infos);
        *dirty = false;
    }
}

/// M1: Insert unassigned PD pairs into existing routes or create new routes.
/// Uses best improvement strategy with a cached insertion matrix: precomputes
/// best insertion for each (pair, route) upfront, then greedily picks the global
/// best, inserts it, and only recomputes the changed route's column.
/// This reduces insertion evaluations from O(N² × R) to O(N × R + N²/2).
#[allow(clippy::similar_names)]
fn apply_m1(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    dirty: &mut bool,
    penalty: &PenaltyState,
    move_ctr: &mut u32,
    route_epoch: &mut Vec<u32>,
) -> bool {
    if sol.unassigned.is_empty() {
        return false;
    }

    ensure_infos(inst, &sol.routes, infos, dirty);

    let nr = sol.routes.len();
    let unassigned: Vec<usize> = sol.unassigned.clone();
    let np = unassigned.len();

    if nr == 0 {
        // No existing routes — go directly to new route creation
        return m1_create_new_route(inst, sol, infos, move_ctr, route_epoch);
    }

    // Precompute insertion cache: ins_cache[pi * nr + ri] = Option<(delta, ep, ed)>
    let mut ins_cache: Vec<Option<(f64, usize, usize)>> = vec![None; np * nr];
    for pi in 0..np {
        let pickup = unassigned[pi];
        let delivery = inst.delivery_of(pickup);
        for ri in 0..nr {
            ins_cache[pi * nr + ri] =
                m1_eval_insertion(inst, &sol.routes[ri], &infos[ri], pickup, delivery, penalty);
        }
    }

    let mut inserted = vec![false; np];
    let mut any_improved = false;

    loop {
        // Find global best (pi, ri) with minimum delta
        let mut best_pi = usize::MAX;
        let mut best_ri = usize::MAX;
        let mut best_ep = 0;
        let mut best_ed = 0;
        let mut best_delta = f64::MAX;

        for pi in 0..np {
            if inserted[pi] {
                continue;
            }
            for ri in 0..nr {
                if let Some((delta, ep, ed)) = ins_cache[pi * nr + ri]
                    && delta < best_delta
                {
                    best_delta = delta;
                    best_pi = pi;
                    best_ri = ri;
                    best_ep = ep;
                    best_ed = ed;
                }
            }
        }

        if best_pi == usize::MAX {
            break; // No feasible insertion found
        }

        // Insert the pair
        let pickup = unassigned[best_pi];
        let delivery = inst.delivery_of(pickup);
        sol.routes[best_ri] =
            build_route_with_pair(&sol.routes[best_ri], pickup, best_ep, delivery, best_ed);
        infos[best_ri].compute(inst, &sol.routes[best_ri]);
        inserted[best_pi] = true;
        *move_ctr += 1;
        route_epoch[best_ri] = *move_ctr;
        any_improved = true;

        // Recompute only the changed route's column for remaining pairs
        for pj in 0..np {
            if inserted[pj] {
                continue;
            }
            let p = unassigned[pj];
            let d = inst.delivery_of(p);
            ins_cache[pj * nr + best_ri] =
                m1_eval_insertion(inst, &sol.routes[best_ri], &infos[best_ri], p, d, penalty);
        }
    }

    // Update sol.unassigned: keep only non-inserted pairs
    sol.unassigned.clear();
    for pi in 0..np {
        if !inserted[pi] {
            sol.unassigned.push(unassigned[pi]);
        }
    }

    // If there are still unassigned pairs, create a new route for the cheapest
    if !sol.unassigned.is_empty() && m1_create_new_route(inst, sol, infos, move_ctr, route_epoch) {
        any_improved = true;
    }

    any_improved
}

/// Evaluate inserting a PD pair into a route, returning delta and positions.
#[inline]
fn m1_eval_insertion(
    inst: &Instance,
    route: &[usize],
    info: &RouteInfo,
    pickup: usize,
    delivery: usize,
    penalty: &PenaltyState,
) -> Option<(f64, usize, usize)> {
    if penalty.is_active() {
        best_pair_insertion_penalized(
            inst,
            route,
            info,
            pickup,
            delivery,
            penalty.penalty_cap,
            penalty.penalty_tw,
        )
        .map(|(pc, ep, ed)| {
            let old_pc = penalty.penalized_cost(info);
            (pc - old_pc, ep, ed)
        })
    } else {
        best_pair_insertion_with_info::<false>(inst, route, info, pickup, delivery)
            .map(|(dist, ep, ed)| (dist - info.distance, ep, ed))
    }
}

/// Create a new route for the cheapest feasible unassigned pair.
fn m1_create_new_route(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    move_ctr: &mut u32,
    route_epoch: &mut Vec<u32>,
) -> bool {
    let mut best_new_pickup: Option<usize> = None;
    let mut best_new_dist = f64::MAX;
    for &pickup in &sol.unassigned {
        let delivery = inst.delivery_of(pickup);
        if let Some(dist) = single_pair_feasible_distance(inst, pickup, delivery)
            && dist < best_new_dist
        {
            best_new_dist = dist;
            best_new_pickup = Some(pickup);
        }
    }
    if let Some(pickup) = best_new_pickup {
        let delivery = inst.delivery_of(pickup);
        sol.routes.push(vec![0, pickup, delivery, 0]);
        if let Some(pos) = sol.unassigned.iter().position(|&p| p == pickup) {
            sol.unassigned.remove(pos);
        }
        while infos.len() < sol.routes.len() {
            infos.push(RouteInfo::new());
        }
        let last = sol.routes.len() - 1;
        infos[last].compute(inst, &sol.routes[last]);
        *move_ctr += 1;
        if route_epoch.len() < sol.routes.len() {
            route_epoch.resize(sol.routes.len(), 0);
        }
        route_epoch[last] = *move_ctr;
        return true;
    }
    false
}

/// M2: Un-assign an assigned PD pair and insert into a different route (or new route).
/// Uses best improvement with hierarchical objective, cached `RouteInfo`, and
/// a cached insertion matrix that avoids recomputing unchanged routes after each move.
#[allow(clippy::similar_names, clippy::cast_precision_loss)]
fn apply_m2(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    dirty: &mut bool,
    penalty: &PenaltyState,
    move_ctr: &mut u32,
    route_epoch: &mut Vec<u32>,
) -> bool {
    ensure_infos(inst, &sol.routes, infos, dirty);

    let nr = sol.routes.len();
    if nr < 2 {
        return false;
    }

    // Collect all pickups with their route indices
    let mut pickups: Vec<(usize, usize)> = Vec::new(); // (route_idx, pickup_id)
    for ri in 0..nr {
        for &v in &sol.routes[ri] {
            if v != 0 && inst.is_pickup(v) {
                pickups.push((ri, v));
            }
        }
    }
    let np = pickups.len();
    if np == 0 {
        return false;
    }

    // Build node-to-route mapping for sparse pre-filtering
    let mut node_to_route: Vec<usize> = vec![usize::MAX; inst.n + 1];
    for ri in 0..nr {
        for &v in &sol.routes[ri] {
            if v != 0 {
                node_to_route[v] = ri;
            }
        }
    }
    let mut candidate_mask: Vec<bool> = vec![false; nr];

    // Build node position map for O(1) removal delta computation
    let mut node_pos: Vec<usize> = vec![0; inst.n + 1];
    for ri in 0..nr {
        for (pos, &v) in sol.routes[ri].iter().enumerate() {
            if v != 0 {
                node_pos[v] = pos;
            }
        }
    }

    // Precompute conflict masks for all full routes (O(8) lookup vs O(route_len))
    let stride = inst.conflict_stride;
    let mut route_masks: Vec<u64> = vec![0u64; nr * stride];
    let mut mask_buf: Vec<u64> = Vec::new();
    for ri in 0..nr {
        inst.build_route_req_mask(&sol.routes[ri], &mut mask_buf);
        route_masks[ri * stride..(ri + 1) * stride].copy_from_slice(&mask_buf);
    }

    // Insertion cache: for each (pickup_idx, route_idx), store delta + positions
    let mut ins_cache: Vec<Option<(f64, usize, usize)>> = vec![None; np * nr];
    let mut ins_valid: Vec<bool> = vec![false; np * nr];

    // Initial fill with sparse pre-filtering
    for pi in 0..np {
        let (_, pickup) = pickups[pi];
        let delivery = inst.delivery_of(pickup);

        // Check if depot-related arcs make all routes candidates
        let all_routes = inst.is_arc_sparse(pickup, delivery)
            || inst.is_arc_sparse(0, pickup)
            || inst.is_arc_sparse(pickup, 0)
            || inst.is_arc_sparse(delivery, 0);

        if all_routes {
            candidate_mask.fill(true);
        } else {
            candidate_mask.fill(false);
            for &u in &inst.sparse_in[pickup] {
                let r = node_to_route[u];
                if r < nr {
                    candidate_mask[r] = true;
                }
            }
            for &u in &inst.sparse_out[pickup] {
                let r = node_to_route[u];
                if r < nr {
                    candidate_mask[r] = true;
                }
            }
            for &u in &inst.sparse_in[delivery] {
                let r = node_to_route[u];
                if r < nr {
                    candidate_mask[r] = true;
                }
            }
            for &u in &inst.sparse_out[delivery] {
                let r = node_to_route[u];
                if r < nr {
                    candidate_mask[r] = true;
                }
            }
        }

        for ri in 0..nr {
            if !candidate_mask[ri]
                || inst.has_conflict_with_mask(&route_masks[ri * stride..(ri + 1) * stride], pickup)
            {
                ins_valid[pi * nr + ri] = true;
                continue;
            }
            if penalty.is_active() {
                let result = best_pair_insertion_penalized(
                    inst,
                    &sol.routes[ri],
                    &infos[ri],
                    pickup,
                    delivery,
                    penalty.penalty_cap,
                    penalty.penalty_tw,
                );
                let old_pc = penalty.penalized_cost(&infos[ri]);
                ins_cache[pi * nr + ri] = result.map(|(pc, ep, ed)| (pc - old_pc, ep, ed));
            } else {
                let result = best_pair_insertion_with_info::<true>(
                    inst,
                    &sol.routes[ri],
                    &infos[ri],
                    pickup,
                    delivery,
                );
                ins_cache[pi * nr + ri] =
                    result.map(|(dist, ep, ed)| (dist - infos[ri].distance, ep, ed));
            }
            ins_valid[pi * nr + ri] = true;
        }
    }

    // Precompute removal costs using O(1) edge distance delta.
    // Avoids O(route_len) route_without_pair_buf + route_distance per pair per scan.
    let mut rem_new_dist: Vec<f64> = vec![0.0; np];
    let mut rem_veh_delta: Vec<i64> = vec![0; np];
    for pi in 0..np {
        let (ri, pickup) = pickups[pi];
        let delivery = inst.delivery_of(pickup);
        let route = &sol.routes[ri];
        if route.len() <= 4 {
            rem_veh_delta[pi] = -1;
        } else {
            let pos_p = node_pos[pickup];
            let pos_d = node_pos[delivery];
            let delta = if pos_d == pos_p + 1 {
                inst.dist(route[pos_p - 1], route[pos_d + 1])
                    - infos[ri].edge_dists[pos_p - 1]
                    - infos[ri].edge_dists[pos_p]
                    - infos[ri].edge_dists[pos_d]
            } else {
                (inst.dist(route[pos_p - 1], route[pos_p + 1])
                    - infos[ri].edge_dists[pos_p - 1]
                    - infos[ri].edge_dists[pos_p])
                    + (inst.dist(route[pos_d - 1], route[pos_d + 1])
                        - infos[ri].edge_dists[pos_d - 1]
                        - infos[ri].edge_dists[pos_d])
            };
            rem_new_dist[pi] = infos[ri].distance + delta;
        }
    }

    let mut any_improved = false;

    // Compute current cost once; update incrementally after moves (O(1) vs O(nr))
    #[allow(clippy::cast_precision_loss)]
    let mut current_cost = sol.routes.len() as f64 * 1e9
        + infos
            .iter()
            .take(sol.routes.len())
            .map(|info| penalty.penalized_cost(info))
            .sum::<f64>();

    loop {
        let mut best_move: Option<(usize, usize, usize, usize)> = None; // (pi, rj, ep, ed)
        let mut best_move_cost = current_cost;

        for pi in 0..np {
            let (src_ri, pickup) = pickups[pi];
            let delivery = inst.delivery_of(pickup);

            // Pre-check if depot arcs make all routes candidates for this pair
            let all_routes = inst.is_arc_sparse(pickup, delivery)
                || inst.is_arc_sparse(0, pickup)
                || inst.is_arc_sparse(pickup, 0)
                || inst.is_arc_sparse(delivery, 0);

            // O(1) removal cost from precomputed edge-delta cache
            let new_pc_ri = rem_new_dist[pi];
            let old_pc_ri = penalty.penalized_cost(&infos[src_ri]);
            let vehicles_delta = rem_veh_delta[pi];

            for rj in 0..sol.routes.len() {
                if rj == src_ri {
                    continue;
                }

                // Lazily recompute invalid cache entries with sparse pre-filtering
                if !ins_valid[pi * nr + rj] {
                    let is_candidate = !inst.has_conflict_with_mask(
                        &route_masks[rj * stride..(rj + 1) * stride],
                        pickup,
                    ) && (all_routes
                        || inst.sparse_in[pickup]
                            .iter()
                            .any(|&u| node_to_route[u] == rj)
                        || inst.sparse_out[pickup]
                            .iter()
                            .any(|&u| node_to_route[u] == rj)
                        || inst.sparse_in[delivery]
                            .iter()
                            .any(|&u| node_to_route[u] == rj)
                        || inst.sparse_out[delivery]
                            .iter()
                            .any(|&u| node_to_route[u] == rj));
                    if is_candidate {
                        if penalty.is_active() {
                            let result = best_pair_insertion_penalized(
                                inst,
                                &sol.routes[rj],
                                &infos[rj],
                                pickup,
                                delivery,
                                penalty.penalty_cap,
                                penalty.penalty_tw,
                            );
                            let old_pc = penalty.penalized_cost(&infos[rj]);
                            ins_cache[pi * nr + rj] =
                                result.map(|(pc, ep, ed)| (pc - old_pc, ep, ed));
                        } else {
                            let result = best_pair_insertion_with_info::<true>(
                                inst,
                                &sol.routes[rj],
                                &infos[rj],
                                pickup,
                                delivery,
                            );
                            ins_cache[pi * nr + rj] =
                                result.map(|(dist, ep, ed)| (dist - infos[rj].distance, ep, ed));
                        }
                    } else {
                        ins_cache[pi * nr + rj] = None;
                    }
                    ins_valid[pi * nr + rj] = true;
                }

                if let Some((delta_rj, ep, ed)) = ins_cache[pi * nr + rj] {
                    let total_delta = (new_pc_ri - old_pc_ri) + delta_rj;
                    let new_total_cost = current_cost + total_delta + vehicles_delta as f64 * 1e9;

                    if new_total_cost < best_move_cost {
                        best_move_cost = new_total_cost;
                        best_move = Some((pi, rj, ep, ed));
                    }
                }
            }
        }

        if let Some((pi, dst, ep, ed)) = best_move {
            let (src, pickup) = pickups[pi];
            let delivery = inst.delivery_of(pickup);

            let new_dst = build_route_with_pair(&sol.routes[dst], pickup, ep, delivery, ed);
            let new_src = route_without_pair(&sol.routes[src], pickup, delivery);

            sol.routes[dst] = new_dst;
            if new_src.len() <= 2 {
                sol.routes.remove(src);
                route_epoch.remove(src);
                // Route removed — indices shifted, cache is invalid
                *dirty = true;
                *move_ctr += 1;
                for e in route_epoch.iter_mut().take(sol.routes.len()) {
                    *e = *move_ctr;
                }
                any_improved = true;
                break;
            }

            sol.routes[src] = new_src;
            let old_pc_src = penalty.penalized_cost(&infos[src]);
            let old_pc_dst = penalty.penalized_cost(&infos[dst]);
            infos[src].compute(inst, &sol.routes[src]);
            infos[dst].compute(inst, &sol.routes[dst]);
            // Incremental cost update: O(1) instead of O(nr)
            current_cost += (penalty.penalized_cost(&infos[src]) - old_pc_src)
                + (penalty.penalized_cost(&infos[dst]) - old_pc_dst);

            // Update pickup's route index and node-to-route mapping
            pickups[pi].0 = dst;
            node_to_route[pickup] = dst;
            node_to_route[delivery] = dst;

            // Invalidate insertion cache columns for changed routes (src and dst)
            for pj in 0..np {
                ins_valid[pj * nr + src] = false;
                ins_valid[pj * nr + dst] = false;
            }

            *move_ctr += 1;
            route_epoch[src] = *move_ctr;
            route_epoch[dst] = *move_ctr;

            // Rebuild conflict masks for changed routes
            inst.build_route_req_mask(&sol.routes[src], &mut mask_buf);
            route_masks[src * stride..(src + 1) * stride].copy_from_slice(&mask_buf);
            inst.build_route_req_mask(&sol.routes[dst], &mut mask_buf);
            route_masks[dst * stride..(dst + 1) * stride].copy_from_slice(&mask_buf);

            // Update node_pos and removal costs for pairs in changed routes
            for (pos, &v) in sol.routes[src].iter().enumerate() {
                if v != 0 {
                    node_pos[v] = pos;
                }
            }
            for (pos, &v) in sol.routes[dst].iter().enumerate() {
                if v != 0 {
                    node_pos[v] = pos;
                }
            }
            for pj in 0..np {
                let (rr, pp) = pickups[pj];
                if rr != src && rr != dst {
                    continue;
                }
                let dd = inst.delivery_of(pp);
                let route = &sol.routes[rr];
                if route.len() <= 4 {
                    rem_new_dist[pj] = 0.0;
                    rem_veh_delta[pj] = -1;
                } else {
                    let pos_p = node_pos[pp];
                    let pos_d = node_pos[dd];
                    let delta = if pos_d == pos_p + 1 {
                        inst.dist(route[pos_p - 1], route[pos_d + 1])
                            - infos[rr].edge_dists[pos_p - 1]
                            - infos[rr].edge_dists[pos_p]
                            - infos[rr].edge_dists[pos_d]
                    } else {
                        (inst.dist(route[pos_p - 1], route[pos_p + 1])
                            - infos[rr].edge_dists[pos_p - 1]
                            - infos[rr].edge_dists[pos_p])
                            + (inst.dist(route[pos_d - 1], route[pos_d + 1])
                                - infos[rr].edge_dists[pos_d - 1]
                                - infos[rr].edge_dists[pos_d])
                    };
                    rem_new_dist[pj] = infos[rr].distance + delta;
                    rem_veh_delta[pj] = 0;
                }
            }

            any_improved = true;
        } else {
            break;
        }
    }

    any_improved
}

/// M3: Swap PD pairs between two routes. First improvement with cached `RouteInfo`.
/// Pre-computes r2_without RouteInfos for all p2 in routes[rj] to avoid redundant
/// recomputation across different p1 values (saves ~93% of RouteInfo::compute calls).
#[allow(clippy::similar_names)]
fn apply_m3(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    _info_buf: &mut RouteInfo,
    info_buf2: &mut RouteInfo,
    dirty: &mut bool,
    route_epoch: &mut [u32],
    clean_epoch: u32,
    move_ctr: &mut u32,
    penalty: &PenaltyState,
) -> bool {
    ensure_infos(inst, &sol.routes, infos, dirty);

    let nr = sol.routes.len();

    // Build node-to-route mapping for sparse pre-filtering
    let mut node_to_route: Vec<usize> = vec![usize::MAX; inst.n + 1];
    for ri in 0..nr {
        for &v in &sol.routes[ri] {
            if v != 0 {
                node_to_route[v] = ri;
            }
        }
    }

    // Precompute conflict masks for all full routes (O(8) lookup vs O(route_len))
    let stride = inst.conflict_stride;
    let mut route_masks: Vec<u64> = vec![0u64; nr * stride];
    let mut mask_buf: Vec<u64> = Vec::new();
    for ri in 0..nr {
        inst.build_route_req_mask(&sol.routes[ri], &mut mask_buf);
        route_masks[ri * stride..(ri + 1) * stride].copy_from_slice(&mask_buf);
    }

    // Reusable buffer for r1_without_pair
    let mut r1_without_buf: Vec<usize> = Vec::new();

    // Per-rj cache for r2_without precomputation.
    // r2_without only depends on sol.routes[rj], not on ri, so we cache it
    // across all ri values for the same rj. Since M3 is first-improvement
    // and returns immediately on success, the cache stays valid.
    struct M3RjCache {
        routes: Vec<Vec<usize>>,
        infos: Vec<RouteInfo>,
        pickups: Vec<usize>,
        count: usize,
        valid: bool,
    }
    let mut rj_cache: Vec<M3RjCache> = (0..nr)
        .map(|_| M3RjCache {
            routes: Vec::new(),
            infos: Vec::new(),
            pickups: Vec::new(),
            count: 0,
            valid: false,
        })
        .collect();

    for ri in 0..nr {
        for rj in (ri + 1)..nr {
            // Skip if neither route changed since last full scan
            if route_epoch[ri] <= clean_epoch && route_epoch[rj] <= clean_epoch {
                continue;
            }

            // Precompute penalized/distance costs of original routes (for cost bound)
            let old_pc_ri = penalty.penalized_cost(&infos[ri]);
            let old_pc_rj = penalty.penalized_cost(&infos[rj]);

            // Lazily compute and cache r2_without for all pickups in routes[rj].
            // Same rj is visited by multiple ri values — cache avoids redundancy.
            if !rj_cache[rj].valid {
                let cache = &mut rj_cache[rj];
                cache.pickups.clear();
                cache.count = 0;
                for idx_j in 0..sol.routes[rj].len() {
                    let p2 = sol.routes[rj][idx_j];
                    if p2 == 0 || !inst.is_pickup(p2) {
                        continue;
                    }
                    let d2 = inst.delivery_of(p2);

                    // Grow cache as needed
                    if cache.count >= cache.routes.len() {
                        cache.routes.push(Vec::new());
                        cache.infos.push(RouteInfo::new());
                    }

                    compute_info_without_pair(
                        inst,
                        &sol.routes[rj],
                        p2,
                        d2,
                        &mut cache.routes[cache.count],
                        &mut cache.infos[cache.count],
                    );
                    if cache.routes[cache.count].len() < 2 {
                        cache.pickups.push(0); // invalid
                        cache.count += 1;
                        continue;
                    }
                    if penalty.is_active() {
                        cache.infos[cache.count].compute_tw_sector(
                            inst,
                            &cache.routes[cache.count],
                            cache.routes[cache.count].len(),
                        );
                    }
                    cache.pickups.push(p2);
                    cache.count += 1;
                }
                cache.valid = true;
            }

            for idx_i in 0..sol.routes[ri].len() {
                let p1 = sol.routes[ri][idx_i];
                if p1 == 0 || !inst.is_pickup(p1) {
                    continue;
                }
                let d1 = inst.delivery_of(p1);

                // Sparse pre-filter: skip if p1/d1 have no sparse arcs to route rj
                let p1_has_sparse_to_rj = inst.is_arc_sparse(p1, d1)
                    || inst.sparse_out[p1].iter().any(|&u| node_to_route[u] == rj)
                    || inst.sparse_in[p1].iter().any(|&u| node_to_route[u] == rj)
                    || inst.sparse_out[d1].iter().any(|&u| node_to_route[u] == rj)
                    || inst.sparse_in[d1].iter().any(|&u| node_to_route[u] == rj);
                if !p1_has_sparse_to_rj {
                    continue;
                }

                compute_info_without_pair(
                    inst,
                    &sol.routes[ri],
                    p1,
                    d1,
                    &mut r1_without_buf,
                    info_buf2,
                );

                if r1_without_buf.len() < 2 {
                    continue;
                }

                if penalty.is_active() {
                    info_buf2.compute_tw_sector(inst, &r1_without_buf, r1_without_buf.len());
                }

                // Lower bound on r1 cost after inserting any p2: the without-pair cost.
                // Inserting a pair can only increase distance, cap_excess, and tw_excess.
                let r1_without_cost = penalty.penalized_cost(info_buf2);

                let r2_count = rj_cache[rj].count;
                for k in 0..r2_count {
                    let p2 = rj_cache[rj].pickups[k];
                    if p2 == 0 {
                        continue;
                    }
                    let d2 = inst.delivery_of(p2);

                    // Sparse pre-filter: skip if p2/d2 have no sparse arcs to route ri
                    let p2_has_sparse_to_ri = inst.is_arc_sparse(p2, d2)
                        || inst.sparse_out[p2].iter().any(|&u| node_to_route[u] == ri)
                        || inst.sparse_in[p2].iter().any(|&u| node_to_route[u] == ri)
                        || inst.sparse_out[d2].iter().any(|&u| node_to_route[u] == ri)
                        || inst.sparse_in[d2].iter().any(|&u| node_to_route[u] == ri);
                    if !p2_has_sparse_to_ri {
                        continue;
                    }

                    // Conflict pre-filters: check BOTH directions before expensive
                    // insertion calls. O(8) bitwise mask checks are cheap vs O(n²) insertion.
                    // p1 must be compatible with rj\p2, and p2 with ri\p1.
                    if inst.has_conflict_with_mask_excluding(
                        &route_masks[rj * stride..(rj + 1) * stride],
                        p1,
                        p2,
                    ) {
                        continue;
                    }
                    if inst.has_conflict_with_mask_excluding(
                        &route_masks[ri * stride..(ri + 1) * stride],
                        p2,
                        p1,
                    ) {
                        continue;
                    }
                    let (pc_new_r2, ep1, ed1) = if penalty.is_active() {
                        match best_pair_insertion_penalized(
                            inst,
                            &rj_cache[rj].routes[k],
                            &rj_cache[rj].infos[k],
                            p1,
                            d1,
                            penalty.penalty_cap,
                            penalty.penalty_tw,
                        ) {
                            Some((pc, e1, e2)) => (pc, e1, e2),
                            None => continue,
                        }
                    } else {
                        match best_pair_insertion_with_info::<true>(
                            inst,
                            &rj_cache[rj].routes[k],
                            &rj_cache[rj].infos[k],
                            p1,
                            d1,
                        ) {
                            Some((dist, e1, e2)) => (dist, e1, e2),
                            None => continue,
                        }
                    };
                    // Cost-based bound: r1_without_cost is a lower bound on r1 after
                    // inserting p2 (adding a pair can only increase dist/cap/tw).
                    // If partial delta with this lower bound is non-negative, skip.
                    if (r1_without_cost - old_pc_ri) + (pc_new_r2 - old_pc_rj) >= -1e-10 {
                        continue;
                    }
                    let ins_r1 = if penalty.is_active() {
                        best_pair_insertion_penalized(
                            inst,
                            &r1_without_buf,
                            info_buf2,
                            p2,
                            d2,
                            penalty.penalty_cap,
                            penalty.penalty_tw,
                        )
                    } else {
                        best_pair_insertion_with_info::<true>(
                            inst,
                            &r1_without_buf,
                            info_buf2,
                            p2,
                            d2,
                        )
                    };
                    if let Some((pc_new_r1, ep2, ed2)) = ins_r1 {
                        let delta = (pc_new_r1 - old_pc_ri) + (pc_new_r2 - old_pc_rj);

                        if delta < -1e-10 {
                            sol.routes[ri] =
                                build_route_with_pair(&r1_without_buf, p2, ep2, d2, ed2);
                            sol.routes[rj] =
                                build_route_with_pair(&rj_cache[rj].routes[k], p1, ep1, d1, ed1);
                            // Incremental update
                            infos[ri].compute(inst, &sol.routes[ri]);
                            infos[rj].compute(inst, &sol.routes[rj]);
                            *move_ctr += 1;
                            route_epoch[ri] = *move_ctr;
                            *move_ctr += 1;
                            route_epoch[rj] = *move_ctr;
                            return true;
                        }
                    }
                }
            }
        }
    }

    false
}

/// M4: Transfer PD pairs: pd1 from r1 -> r2, pd2 from r2 -> r3.
/// First improvement with hierarchical objective and cached `RouteInfo`.
/// Pre-computes rj_without RouteInfos for all p2 in routes[rj] to avoid
/// redundant recomputation across different p1 values.
///
/// Optimization: precomputes the insertion of every pair into every route once,
/// eliminating redundant `best_pair_insertion_with_info` calls in the inner rk loop.
/// This reduces O(nr² × pairs² × nr) insertion calls to O(pairs × nr) upfront.
#[allow(clippy::similar_names, clippy::cast_precision_loss)]
fn apply_m4(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    _info_buf: &mut RouteInfo,
    dirty: &mut bool,
    penalty: &PenaltyState,
    move_ctr: &mut u32,
    route_epoch: &mut Vec<u32>,
    clean_epoch: u32,
) -> bool {
    ensure_infos(inst, &sol.routes, infos, dirty);

    let nr = sol.routes.len();
    if nr < 3 {
        return false; // M4 needs 3 distinct routes (ri, rj, rk)
    }

    // Quick check: if no route was modified since last M4 scan, skip entirely
    // (moved before expensive precomputations)
    if route_epoch.iter().take(nr).all(|&e| e <= clean_epoch) {
        return false;
    }

    // Collect all assigned pickups with their route indices
    let mut all_pickups: Vec<usize> = Vec::new();
    let mut pickup_route: Vec<usize> = vec![usize::MAX; inst.n + 1];
    for ri in 0..nr {
        for &v in &sol.routes[ri] {
            if v != 0 && inst.is_pickup(v) {
                all_pickups.push(v);
                pickup_route[v] = ri;
            }
        }
    }
    let np = all_pickups.len();
    if np == 0 {
        return false;
    }

    // Build pickup -> global index mapping
    let mut pickup_to_idx: Vec<usize> = vec![usize::MAX; inst.n + 1];
    for (idx, &p) in all_pickups.iter().enumerate() {
        pickup_to_idx[p] = idx;
    }

    // Build node-to-route and node-to-position mappings (needed for O(1) removal delta)
    let mut node_to_route: Vec<usize> = vec![usize::MAX; inst.n + 1];
    let mut node_pos: Vec<usize> = vec![0; inst.n + 1];
    for ri in 0..nr {
        for (pos, &v) in sol.routes[ri].iter().enumerate() {
            if v != 0 {
                node_to_route[v] = ri;
                node_pos[v] = pos;
            }
        }
    }

    // Precompute conflict masks for all full routes (O(8) lookup vs O(route_len))
    let stride = inst.conflict_stride;
    let mut route_masks: Vec<u64> = vec![0u64; nr * stride];
    let mut mask_buf: Vec<u64> = Vec::new();
    for ri in 0..nr {
        inst.build_route_req_mask(&sol.routes[ri], &mut mask_buf);
        route_masks[ri * stride..(ri + 1) * stride].copy_from_slice(&mask_buf);
    }

    // Precompute removal costs for ALL pickups using O(1) edge-distance delta.
    let mut removal_cost: Vec<f64> = vec![0.0; np];
    let mut removal_veh_delta: Vec<i64> = vec![0; np];
    let mut removal_ri_delta: Vec<f64> = vec![0.0; np];
    for pi in 0..np {
        let p = all_pickups[pi];
        let d = inst.delivery_of(p);
        let ri = pickup_route[p];
        let route = &sol.routes[ri];
        let old_cost = if penalty.is_active() {
            penalty.penalized_cost(&infos[ri])
        } else {
            infos[ri].distance
        };
        if route.len() <= 4 {
            removal_veh_delta[pi] = -1;
            removal_ri_delta[pi] = -old_cost;
        } else {
            let pos_p = node_pos[p];
            let pos_d = node_pos[d];
            let delta = if pos_d == pos_p + 1 {
                inst.dist(route[pos_p - 1], route[pos_d + 1])
                    - infos[ri].edge_dists[pos_p - 1]
                    - infos[ri].edge_dists[pos_p]
                    - infos[ri].edge_dists[pos_d]
            } else {
                (inst.dist(route[pos_p - 1], route[pos_p + 1])
                    - infos[ri].edge_dists[pos_p - 1]
                    - infos[ri].edge_dists[pos_p])
                    + (inst.dist(route[pos_d - 1], route[pos_d + 1])
                        - infos[ri].edge_dists[pos_d - 1]
                        - infos[ri].edge_dists[pos_d])
            };
            let nc = infos[ri].distance + delta;
            removal_cost[pi] = nc;
            removal_ri_delta[pi] = nc - old_cost;
        }
    }

    // Per-route aggregates for (ri, rj) pair-level pruning
    let mut min_ri_delta_per_route: Vec<f64> = vec![f64::MAX; nr];
    let mut has_veh_reduction: Vec<bool> = vec![false; nr];
    for pi in 0..np {
        let ri = pickup_route[all_pickups[pi]];
        if removal_ri_delta[pi] < min_ri_delta_per_route[ri] {
            min_ri_delta_per_route[ri] = removal_ri_delta[pi];
        }
        if removal_veh_delta[pi] == -1 {
            has_veh_reduction[ri] = true;
        }
    }

    let mut candidate_mask: Vec<bool> = vec![false; nr];

    // Precompute insertion cache: for each pair (pi) and each route (rk),
    // the best sparse insertion delta and positions. Sparse pre-filtering
    // skips (pair, route) combinations with no sparse arc interaction,
    // since best_pair_insertion_with_info::<true> would return None anyway.
    let mut ins_cache: Vec<Option<(f64, usize, usize)>> = vec![None; np * nr];
    let mut best_ins_delta: Vec<f64> = vec![f64::MAX; np]; // global min delta per pair

    for pi in 0..np {
        let p = all_pickups[pi];
        let d = inst.delivery_of(p);
        let own_ri = pickup_route[p];

        // Sparse pre-filter: identify candidate routes for this pair
        let all_routes = inst.is_arc_sparse(p, d)
            || inst.is_arc_sparse(0, p)
            || inst.is_arc_sparse(p, 0)
            || inst.is_arc_sparse(d, 0);

        if all_routes {
            candidate_mask.fill(true);
        } else {
            candidate_mask.fill(false);
            for &u in &inst.sparse_in[p] {
                let r = node_to_route[u];
                if r < nr {
                    candidate_mask[r] = true;
                }
            }
            for &u in &inst.sparse_out[p] {
                let r = node_to_route[u];
                if r < nr {
                    candidate_mask[r] = true;
                }
            }
            for &u in &inst.sparse_in[d] {
                let r = node_to_route[u];
                if r < nr {
                    candidate_mask[r] = true;
                }
            }
            for &u in &inst.sparse_out[d] {
                let r = node_to_route[u];
                if r < nr {
                    candidate_mask[r] = true;
                }
            }
        }

        for rk in 0..nr {
            if rk == own_ri {
                continue;
            }
            if !candidate_mask[rk]
                || inst.has_conflict_with_mask(&route_masks[rk * stride..(rk + 1) * stride], p)
            {
                // No sparse arc interaction or request conflict — insertion infeasible
                continue;
            }
            let cached = if penalty.is_active() {
                let result = best_pair_insertion_penalized(
                    inst,
                    &sol.routes[rk],
                    &infos[rk],
                    p,
                    d,
                    penalty.penalty_cap,
                    penalty.penalty_tw,
                );
                let old_pc = penalty.penalized_cost(&infos[rk]);
                result.map(|(pc, ep, ed)| (pc - old_pc, ep, ed))
            } else {
                let result =
                    best_pair_insertion_with_info::<true>(inst, &sol.routes[rk], &infos[rk], p, d);
                result.map(|(dist, ep, ed)| (dist - infos[rk].distance, ep, ed))
            };
            ins_cache[pi * nr + rk] = cached;
            if let Some((delta, _, _)) = cached
                && delta < best_ins_delta[pi]
            {
                best_ins_delta[pi] = delta;
            }
        }
    }

    // Per-route min insertion delta (for pairs IN that route) for pair-level pruning
    let mut min_ins_delta_per_route: Vec<f64> = vec![f64::MAX; nr];
    for pi in 0..np {
        let ri = pickup_route[all_pickups[pi]];
        if best_ins_delta[pi] < min_ins_delta_per_route[ri] {
            min_ins_delta_per_route[ri] = best_ins_delta[pi];
        }
    }

    // Reusable buffer
    let mut r1_without_buf: Vec<usize> = Vec::new();

    // Per-rj cache for rj_without precomputation.
    // rj_without only depends on sol.routes[rj], not on ri, so we cache it
    // across all ri values for the same rj. Saves ~(nr-1) redundant
    // compute_info_without_pair calls per rj per scan.
    struct RjWoCache {
        routes: Vec<Vec<usize>>,
        infos: Vec<RouteInfo>,
        pickups: Vec<usize>, // 0 = invalid (route too short after removal)
        count: usize,
        valid: bool,
    }
    let mut rj_wo: Vec<RjWoCache> = (0..nr)
        .map(|_| RjWoCache {
            routes: Vec::new(),
            infos: Vec::new(),
            pickups: Vec::new(),
            count: 0,
            valid: false,
        })
        .collect();

    let mut any_improved = false;

    // Loop: scan for improving move, apply, incrementally update cache, repeat.
    // Avoids expensive full cache rebuild between consecutive M4 moves.
    'restart: loop {
        for ri in 0..nr {
            for rj in 0..nr {
                if rj == ri {
                    continue;
                }

                // Route-pair level pruning: if min removal from ri + min insertion of
                // any pair from rj can't improve, skip this entire (ri, rj) pair
                // (avoids expensive rj_without precomputation and inner loop scan)
                if !has_veh_reduction[ri]
                    && min_ri_delta_per_route[ri] + min_ins_delta_per_route[rj] >= -1e-10
                {
                    continue;
                }

                // Lazily compute and cache rj_without for this rj value.
                // Same rj may be visited by multiple ri values — cache avoids redundancy.
                if !rj_wo[rj].valid {
                    let cache = &mut rj_wo[rj];
                    cache.pickups.clear();
                    cache.count = 0;
                    for idx_j in 0..sol.routes[rj].len() {
                        let p2 = sol.routes[rj][idx_j];
                        if p2 == 0 || !inst.is_pickup(p2) {
                            continue;
                        }
                        let d2 = inst.delivery_of(p2);

                        if cache.count >= cache.routes.len() {
                            cache.routes.push(Vec::new());
                            cache.infos.push(RouteInfo::new());
                        }

                        compute_info_without_pair(
                            inst,
                            &sol.routes[rj],
                            p2,
                            d2,
                            &mut cache.routes[cache.count],
                            &mut cache.infos[cache.count],
                        );
                        if cache.routes[cache.count].len() < 2 {
                            cache.pickups.push(0);
                            cache.count += 1;
                            continue;
                        }
                        if penalty.is_active() {
                            cache.infos[cache.count].compute_tw_sector(
                                inst,
                                &cache.routes[cache.count],
                                cache.routes[cache.count].len(),
                            );
                        }
                        cache.pickups.push(p2);
                        cache.count += 1;
                    }
                    cache.valid = true;
                }

                let rj_count = rj_wo[rj].count;

                for idx_i in 0..sol.routes[ri].len() {
                    let p1 = sol.routes[ri][idx_i];
                    if p1 == 0 || !inst.is_pickup(p1) {
                        continue;
                    }
                    let d1 = inst.delivery_of(p1);
                    let p1_idx = pickup_to_idx[p1];

                    // Use precomputed removal costs (O(1) lookup instead of O(route_len))
                    let vehicles_delta = removal_veh_delta[p1_idx];
                    let ri_delta = removal_ri_delta[p1_idx];

                    for k in 0..rj_count {
                        let p2 = rj_wo[rj].pickups[k];
                        if p2 == 0 {
                            continue;
                        }

                        let p2_idx = pickup_to_idx[p2];

                        // Tighter pruning: if partial_delta from ri alone + best possible
                        // p1 insertion (>= 0) + best possible p2 insertion can't improve,
                        // skip the expensive p1 insertion call
                        if vehicles_delta == 0 && ri_delta + best_ins_delta[p2_idx] >= -1e-10 {
                            continue;
                        }

                        // Conflict pre-filter: skip if p1 conflicts with remaining nodes in rj
                        // 2-level: O(8) full mask check, then O(8) derived check excluding p2
                        if inst.has_conflict_with_mask_excluding(
                            &route_masks[rj * stride..(rj + 1) * stride],
                            p1,
                            p2,
                        ) {
                            continue;
                        }

                        let (cost_new_rj, ep1, ed1) = if penalty.is_active() {
                            match best_pair_insertion_penalized(
                                inst,
                                &rj_wo[rj].routes[k],
                                &rj_wo[rj].infos[k],
                                p1,
                                d1,
                                penalty.penalty_cap,
                                penalty.penalty_tw,
                            ) {
                                Some((pc, e1, e2)) => (pc, e1, e2),
                                None => continue,
                            }
                        } else {
                            match best_pair_insertion_with_info::<true>(
                                inst,
                                &rj_wo[rj].routes[k],
                                &rj_wo[rj].infos[k],
                                p1,
                                d1,
                            ) {
                                Some((d, e1, e2)) => (d, e1, e2),
                                None => continue,
                            }
                        };

                        let old_rj_cost = if penalty.is_active() {
                            penalty.penalized_cost(&infos[rj])
                        } else {
                            infos[rj].distance
                        };
                        let partial_delta = ri_delta + (cost_new_rj - old_rj_cost);
                        if vehicles_delta == 0 && partial_delta >= -1e-10 {
                            continue;
                        }

                        // Even tighter: partial_delta + best p2 insertion can't improve
                        if vehicles_delta == 0 && partial_delta + best_ins_delta[p2_idx] >= -1e-10 {
                            continue;
                        }

                        // Use precomputed insertion cache for rk loop (O(1) per entry)
                        for rk in 0..nr {
                            if rk == ri || rk == rj {
                                continue;
                            }

                            if let Some((delta_rk, ep2, ed2)) = ins_cache[p2_idx * nr + rk] {
                                let dist_delta = partial_delta + delta_rk;
                                let cost_delta = dist_delta + vehicles_delta as f64 * 1e9;

                                if cost_delta < -1e-10 {
                                    let d2 = inst.delivery_of(p2);
                                    sol.routes[rj] = build_route_with_pair(
                                        &rj_wo[rj].routes[k],
                                        p1,
                                        ep1,
                                        d1,
                                        ed1,
                                    );
                                    sol.routes[rk] =
                                        build_route_with_pair(&sol.routes[rk], p2, ep2, d2, ed2);
                                    if vehicles_delta == -1 {
                                        sol.routes.remove(ri);
                                        route_epoch.remove(ri);
                                        // Route removed — indices shifted, bail out
                                        *dirty = true;
                                        *move_ctr += 1;
                                        for e in route_epoch.iter_mut().take(sol.routes.len()) {
                                            *e = *move_ctr;
                                        }
                                        return true;
                                    }

                                    // Non-removal move: apply, incrementally update cache, restart scan
                                    route_without_pair_buf(
                                        &sol.routes[ri],
                                        p1,
                                        d1,
                                        &mut r1_without_buf,
                                    );
                                    std::mem::swap(&mut sol.routes[ri], &mut r1_without_buf);
                                    infos[ri].compute(inst, &sol.routes[ri]);
                                    infos[rj].compute(inst, &sol.routes[rj]);
                                    infos[rk].compute(inst, &sol.routes[rk]);
                                    *move_ctr += 1;
                                    route_epoch[ri] = *move_ctr;
                                    route_epoch[rj] = *move_ctr;
                                    route_epoch[rk] = *move_ctr;
                                    any_improved = true;

                                    // --- Incremental cache update ---

                                    // Invalidate rj_wo cache for changed routes
                                    rj_wo[ri].valid = false;
                                    rj_wo[rj].valid = false;
                                    rj_wo[rk].valid = false;

                                    // Update node_to_route for moved nodes
                                    node_to_route[p1] = rj;
                                    node_to_route[d1] = rj;
                                    node_to_route[p2] = rk;
                                    node_to_route[d2] = rk;

                                    // Update pickup_route for moved pairs
                                    pickup_route[p1] = rj;
                                    pickup_route[p2] = rk;

                                    // Rebuild conflict masks for changed routes
                                    let changed = [ri, rj, rk];
                                    for &cr in &changed {
                                        inst.build_route_req_mask(&sol.routes[cr], &mut mask_buf);
                                        route_masks[cr * stride..(cr + 1) * stride]
                                            .copy_from_slice(&mask_buf);
                                    }

                                    // Recompute removal costs for pairs in changed routes
                                    // using O(1) edge-distance delta
                                    for &cr in &changed {
                                        for (pos, &v) in sol.routes[cr].iter().enumerate() {
                                            if v != 0 {
                                                node_pos[v] = pos;
                                            }
                                        }
                                    }
                                    for &cr in &changed {
                                        let route = &sol.routes[cr];
                                        let oc = if penalty.is_active() {
                                            penalty.penalized_cost(&infos[cr])
                                        } else {
                                            infos[cr].distance
                                        };
                                        for &v in route {
                                            if v != 0 && inst.is_pickup(v) {
                                                let pidx = pickup_to_idx[v];
                                                let dv = inst.delivery_of(v);
                                                if route.len() <= 4 {
                                                    removal_cost[pidx] = 0.0;
                                                    removal_veh_delta[pidx] = -1;
                                                    removal_ri_delta[pidx] = -oc;
                                                } else {
                                                    let pp = node_pos[v];
                                                    let pd = node_pos[dv];
                                                    let delta = if pd == pp + 1 {
                                                        inst.dist(route[pp - 1], route[pd + 1])
                                                            - infos[cr].edge_dists[pp - 1]
                                                            - infos[cr].edge_dists[pp]
                                                            - infos[cr].edge_dists[pd]
                                                    } else {
                                                        (inst.dist(route[pp - 1], route[pp + 1])
                                                            - infos[cr].edge_dists[pp - 1]
                                                            - infos[cr].edge_dists[pp])
                                                            + (inst
                                                                .dist(route[pd - 1], route[pd + 1])
                                                                - infos[cr].edge_dists[pd - 1]
                                                                - infos[cr].edge_dists[pd])
                                                    };
                                                    let nc = infos[cr].distance + delta;
                                                    removal_cost[pidx] = nc;
                                                    removal_veh_delta[pidx] = 0;
                                                    removal_ri_delta[pidx] = nc - oc;
                                                }
                                            }
                                        }
                                    }

                                    // Recompute per-route aggregates
                                    min_ri_delta_per_route.fill(f64::MAX);
                                    has_veh_reduction.fill(false);
                                    for pi in 0..np {
                                        let r = pickup_route[all_pickups[pi]];
                                        if removal_ri_delta[pi] < min_ri_delta_per_route[r] {
                                            min_ri_delta_per_route[r] = removal_ri_delta[pi];
                                        }
                                        if removal_veh_delta[pi] == -1 {
                                            has_veh_reduction[r] = true;
                                        }
                                    }

                                    // Recompute ins_cache columns for changed routes
                                    for &cr in &changed {
                                        for pi in 0..np {
                                            let p = all_pickups[pi];
                                            let dp = inst.delivery_of(p);
                                            let own = pickup_route[p];
                                            if cr == own {
                                                ins_cache[pi * nr + cr] = None;
                                                continue;
                                            }
                                            let is_cand = !inst.has_conflict_with_mask(
                                                &route_masks[cr * stride..(cr + 1) * stride],
                                                p,
                                            ) && (inst.is_arc_sparse(p, dp)
                                                || inst.is_arc_sparse(0, p)
                                                || inst.is_arc_sparse(p, 0)
                                                || inst.is_arc_sparse(dp, 0)
                                                || inst.sparse_in[p]
                                                    .iter()
                                                    .any(|&u| node_to_route[u] == cr)
                                                || inst.sparse_out[p]
                                                    .iter()
                                                    .any(|&u| node_to_route[u] == cr)
                                                || inst.sparse_in[dp]
                                                    .iter()
                                                    .any(|&u| node_to_route[u] == cr)
                                                || inst.sparse_out[dp]
                                                    .iter()
                                                    .any(|&u| node_to_route[u] == cr));
                                            if !is_cand {
                                                ins_cache[pi * nr + cr] = None;
                                                continue;
                                            }
                                            ins_cache[pi * nr + cr] = if penalty.is_active() {
                                                let res = best_pair_insertion_penalized(
                                                    inst,
                                                    &sol.routes[cr],
                                                    &infos[cr],
                                                    p,
                                                    dp,
                                                    penalty.penalty_cap,
                                                    penalty.penalty_tw,
                                                );
                                                let opc = penalty.penalized_cost(&infos[cr]);
                                                res.map(|(pc, ep, ed)| (pc - opc, ep, ed))
                                            } else {
                                                let res = best_pair_insertion_with_info::<true>(
                                                    inst,
                                                    &sol.routes[cr],
                                                    &infos[cr],
                                                    p,
                                                    dp,
                                                );
                                                res.map(|(dist, ep, ed)| {
                                                    (dist - infos[cr].distance, ep, ed)
                                                })
                                            };
                                        }
                                    }

                                    // Recompute best_ins_delta for all pairs
                                    best_ins_delta.fill(f64::MAX);
                                    for pi in 0..np {
                                        let own = pickup_route[all_pickups[pi]];
                                        for rk_i in 0..nr {
                                            if rk_i == own {
                                                continue;
                                            }
                                            if let Some((delta, _, _)) = ins_cache[pi * nr + rk_i]
                                                && delta < best_ins_delta[pi]
                                            {
                                                best_ins_delta[pi] = delta;
                                            }
                                        }
                                    }

                                    // Recompute min_ins_delta_per_route
                                    min_ins_delta_per_route.fill(f64::MAX);
                                    for pi in 0..np {
                                        let r = pickup_route[all_pickups[pi]];
                                        if best_ins_delta[pi] < min_ins_delta_per_route[r] {
                                            min_ins_delta_per_route[r] = best_ins_delta[pi];
                                        }
                                    }

                                    continue 'restart;
                                }
                            }
                        }
                    }
                }
            }
        }
        break; // Full scan completed without improvement
    } // 'restart loop

    any_improved
}

/// M5: Intra-route pair relocation — move a PD pair to a better position within its route.
/// Best improvement strategy with per-route caching: after each move, only rescan the
/// modified route (other routes' best intra-route moves are unchanged).
#[allow(clippy::similar_names)]
fn apply_m5(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    info_buf: &mut RouteInfo,
    dirty: &mut bool,
    move_ctr: &mut u32,
    route_epoch: &mut [u32],
    penalty: &PenaltyState,
    route_clean: &mut [u32],
) -> bool {
    ensure_infos(inst, &sol.routes, infos, dirty);

    let nr = sol.routes.len();
    let mut r_without_buf: Vec<usize> = Vec::new();

    // Per-route cache of best intra-route move
    let mut rt_delta: Vec<f64> = vec![f64::MAX; nr];
    let mut rt_pickup: Vec<usize> = vec![0; nr];
    let mut rt_ep: Vec<usize> = vec![0; nr];
    let mut rt_ed: Vec<usize> = vec![0; nr];
    let mut rt_valid: Vec<bool> = vec![false; nr];

    // Pre-mark routes unchanged since last M5 scan as already valid (skip scan)
    for ri in 0..nr {
        if route_epoch[ri] <= route_clean[ri] {
            rt_valid[ri] = true;
        }
    }

    let mut any_improved = false;

    loop {
        // Compute best move for invalid (unscanned) routes
        for ri in 0..nr {
            if rt_valid[ri] {
                continue;
            }
            rt_delta[ri] = f64::MAX;
            rt_valid[ri] = true;

            let route = &sol.routes[ri];
            if route.len() < 6 {
                continue;
            }

            let orig_cost = if penalty.is_active() {
                penalty.penalized_cost(&infos[ri])
            } else {
                infos[ri].distance
            };
            let mut best_delta = -1e-10;
            for idx in 1..route.len() - 1 {
                let p = route[idx];
                if !inst.is_pickup(p) {
                    continue;
                }
                let d = inst.delivery_of(p);

                compute_info_without_pair(inst, route, p, d, &mut r_without_buf, info_buf);
                if penalty.is_active() {
                    info_buf.compute_tw_sector(inst, &r_without_buf, r_without_buf.len());
                }
                let result = if penalty.is_active() {
                    best_pair_insertion_penalized(
                        inst,
                        &r_without_buf,
                        info_buf,
                        p,
                        d,
                        penalty.penalty_cap,
                        penalty.penalty_tw,
                    )
                } else {
                    best_pair_insertion_with_info::<true>(inst, &r_without_buf, info_buf, p, d)
                };
                if let Some((new_cost, ep, ed)) = result {
                    let delta = new_cost - orig_cost;
                    if delta < best_delta {
                        best_delta = delta;
                        rt_delta[ri] = delta;
                        rt_pickup[ri] = p;
                        rt_ep[ri] = ep;
                        rt_ed[ri] = ed;
                    }
                }
            }
        }

        // Find global best across all routes
        let mut best_ri = usize::MAX;
        let mut best_delta = -1e-10;
        for (ri, &delta) in rt_delta.iter().enumerate().take(nr) {
            if delta < best_delta {
                best_delta = delta;
                best_ri = ri;
            }
        }

        if best_ri == usize::MAX {
            break;
        }

        // Apply the best move
        let p = rt_pickup[best_ri];
        let d = inst.delivery_of(p);
        route_without_pair_buf(&sol.routes[best_ri], p, d, &mut r_without_buf);
        sol.routes[best_ri] =
            build_route_with_pair(&r_without_buf, p, rt_ep[best_ri], d, rt_ed[best_ri]);
        infos[best_ri].compute(inst, &sol.routes[best_ri]);
        *move_ctr += 1;
        route_epoch[best_ri] = *move_ctr;

        // Only the modified route needs rescan
        rt_valid[best_ri] = false;
        any_improved = true;
    }

    // Update route_clean for all scanned (valid) routes
    for ri in 0..nr {
        if rt_valid[ri] {
            route_clean[ri] = *move_ctr;
        }
    }

    any_improved
}

/// M6: Intra-route 2-pair relocation — remove 2 PD pairs, reinsert both optimally trying
/// both insertion orders. Best improvement with per-route caching (like M5).
/// Escapes saddle points that single-pair M5 cannot reach.
#[allow(clippy::similar_names, clippy::too_many_lines)]
fn apply_m6(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    info_buf: &mut RouteInfo,
    info_buf2: &mut RouteInfo,
    dirty: &mut bool,
    move_ctr: &mut u32,
    route_epoch: &mut [u32],
    penalty: &PenaltyState,
    route_clean: &mut [u32],
) -> bool {
    ensure_infos(inst, &sol.routes, infos, dirty);

    let nr = sol.routes.len();

    // Work buffers (reused across all iterations)
    let mut buf_without_both: Vec<usize> = Vec::new();
    let mut buf_intermediate: Vec<usize> = Vec::new();
    let mut pickups_buf: Vec<(f64, usize)> = Vec::new(); // (removal_delta, pickup)
    let mut pos_in_route: Vec<usize> = vec![0; inst.n + 1]; // node → position in route

    // Per-route cache of best 2-pair intra-route move
    let mut rt_delta: Vec<f64> = vec![f64::MAX; nr];
    let mut rt_p1: Vec<usize> = vec![0; nr];
    let mut rt_p2: Vec<usize> = vec![0; nr];
    let mut rt_first_is_p1: Vec<bool> = vec![false; nr];
    let mut rt_ep_first: Vec<usize> = vec![0; nr];
    let mut rt_ed_first: Vec<usize> = vec![0; nr];
    let mut rt_ep_second: Vec<usize> = vec![0; nr];
    let mut rt_ed_second: Vec<usize> = vec![0; nr];
    let mut rt_valid: Vec<bool> = vec![false; nr];

    // Pre-mark routes unchanged since last M6 scan as already valid (skip scan)
    for ri in 0..nr {
        if route_epoch[ri] <= route_clean[ri] {
            rt_valid[ri] = true;
        }
    }

    let mut any_improved = false;

    loop {
        // Scan invalid (unscanned) routes
        for ri in 0..nr {
            if rt_valid[ri] {
                continue;
            }
            rt_delta[ri] = f64::MAX;
            rt_valid[ri] = true;

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

            // Build position map for this route
            for (pos, &v) in route.iter().enumerate() {
                pos_in_route[v] = pos;
            }

            // Collect pickups with precomputed O(1) removal deltas, sorted by delta
            pickups_buf.clear();
            let edge_dists_ri = &infos[ri].edge_dists;
            for &v in &route[1..route.len() - 1] {
                if !inst.is_pickup(v) {
                    continue;
                }
                let dv = inst.delivery_of(v);
                let pos_p = pos_in_route[v];
                let pos_d = pos_in_route[dv];
                let removal_delta = if pos_d == pos_p + 1 {
                    inst.dist(route[pos_p - 1], route[pos_d + 1])
                        - edge_dists_ri[pos_p - 1]
                        - edge_dists_ri[pos_p]
                        - edge_dists_ri[pos_d]
                } else {
                    (inst.dist(route[pos_p - 1], route[pos_p + 1])
                        - edge_dists_ri[pos_p - 1]
                        - edge_dists_ri[pos_p])
                        + (inst.dist(route[pos_d - 1], route[pos_d + 1])
                            - edge_dists_ri[pos_d - 1]
                            - edge_dists_ri[pos_d])
                };
                pickups_buf.push((removal_delta, v));
            }
            if pickups_buf.len() < 2 {
                continue;
            }
            // Sort by removal delta (most negative first) for better pruning
            pickups_buf.sort_unstable_by(|a, b| a.0.partial_cmp(&b.0).unwrap());

            let orig_cost = if penalty.is_active() {
                penalty.penalized_cost(&infos[ri])
            } else {
                infos[ri].distance
            };
            let mut best_delta = -1e-10;

            for i in 0..pickups_buf.len() {
                let (_delta_i, p1) = pickups_buf[i];
                let d1 = inst.delivery_of(p1);

                for &(_delta_j, p2) in &pickups_buf[(i + 1)..] {
                    let d2 = inst.delivery_of(p2);

                    // O(1) exact combined distance removal delta for removing 4 nodes.
                    // Groups consecutive positions into runs, handles all adjacency patterns.
                    // Since reinsertion can only add distance (Euclidean triangle inequality),
                    // this is a safe lower bound on the final delta in non-penalty mode.
                    if !penalty.is_active() {
                        let mut positions = [
                            pos_in_route[p1],
                            pos_in_route[d1],
                            pos_in_route[p2],
                            pos_in_route[d2],
                        ];
                        positions.sort_unstable();
                        let mut removal_delta = 0.0;
                        let mut pi = 0;
                        while pi < 4 {
                            let run_start = positions[pi];
                            let mut run_end = run_start;
                            while pi + 1 < 4 && positions[pi + 1] == run_end + 1 {
                                pi += 1;
                                run_end = positions[pi];
                            }
                            removal_delta += inst.dist(route[run_start - 1], route[run_end + 1]);
                            for &ed in &edge_dists_ri[(run_start - 1)..=run_end] {
                                removal_delta -= ed;
                            }
                            pi += 1;
                        }
                        if removal_delta >= best_delta {
                            continue;
                        }
                    }

                    // Remove both pairs from route
                    buf_without_both.clear();
                    for &v in route {
                        if v != p1 && v != d1 && v != p2 && v != d2 {
                            buf_without_both.push(v);
                        }
                    }

                    if penalty.is_active() {
                        info_buf.compute(inst, &buf_without_both);
                    } else {
                        info_buf.compute_light(inst, &buf_without_both);
                    }

                    // Post-compute bound: reduced cost is a lower bound on final cost
                    // (reinserting pairs can only increase dist, cap_excess, tw_excess).
                    let reduced_cost = if penalty.is_active() {
                        penalty.penalized_cost(info_buf)
                    } else {
                        info_buf.distance
                    };
                    if reduced_cost - orig_cost >= best_delta {
                        continue;
                    }

                    // Order A: insert (p1,d1) first, then (p2,d2)
                    let ins_a1 = if penalty.is_active() {
                        best_pair_insertion_penalized(
                            inst,
                            &buf_without_both,
                            info_buf,
                            p1,
                            d1,
                            penalty.penalty_cap,
                            penalty.penalty_tw,
                        )
                    } else {
                        best_pair_insertion_with_info::<false>(
                            inst,
                            &buf_without_both,
                            info_buf,
                            p1,
                            d1,
                        )
                    };
                    if let Some((_, ep1, ed1)) = ins_a1 {
                        compute_info_with_pair(
                            inst,
                            &buf_without_both,
                            p1,
                            ep1,
                            d1,
                            ed1,
                            &mut buf_intermediate,
                            info_buf2,
                        );
                        if penalty.is_active() {
                            info_buf2.compute_tw_sector(
                                inst,
                                &buf_intermediate,
                                buf_intermediate.len(),
                            );
                        }
                        let ins_a2 = if penalty.is_active() {
                            best_pair_insertion_penalized(
                                inst,
                                &buf_intermediate,
                                info_buf2,
                                p2,
                                d2,
                                penalty.penalty_cap,
                                penalty.penalty_tw,
                            )
                        } else {
                            best_pair_insertion_with_info::<false>(
                                inst,
                                &buf_intermediate,
                                info_buf2,
                                p2,
                                d2,
                            )
                        };
                        if let Some((cost_final, ep2, ed2)) = ins_a2 {
                            let delta = cost_final - orig_cost;
                            if delta < best_delta {
                                best_delta = delta;
                                rt_delta[ri] = delta;
                                rt_p1[ri] = p1;
                                rt_p2[ri] = p2;
                                rt_first_is_p1[ri] = true;
                                rt_ep_first[ri] = ep1;
                                rt_ed_first[ri] = ed1;
                                rt_ep_second[ri] = ep2;
                                rt_ed_second[ri] = ed2;
                            }
                        }
                    }

                    // Order B: insert (p2,d2) first, then (p1,d1)
                    // info_buf still holds RouteInfo for buf_without_both (unchanged)
                    let ins_b1 = if penalty.is_active() {
                        best_pair_insertion_penalized(
                            inst,
                            &buf_without_both,
                            info_buf,
                            p2,
                            d2,
                            penalty.penalty_cap,
                            penalty.penalty_tw,
                        )
                    } else {
                        best_pair_insertion_with_info::<false>(
                            inst,
                            &buf_without_both,
                            info_buf,
                            p2,
                            d2,
                        )
                    };
                    if let Some((_, ep2, ed2)) = ins_b1 {
                        compute_info_with_pair(
                            inst,
                            &buf_without_both,
                            p2,
                            ep2,
                            d2,
                            ed2,
                            &mut buf_intermediate,
                            info_buf2,
                        );
                        if penalty.is_active() {
                            info_buf2.compute_tw_sector(
                                inst,
                                &buf_intermediate,
                                buf_intermediate.len(),
                            );
                        }
                        let ins_b2 = if penalty.is_active() {
                            best_pair_insertion_penalized(
                                inst,
                                &buf_intermediate,
                                info_buf2,
                                p1,
                                d1,
                                penalty.penalty_cap,
                                penalty.penalty_tw,
                            )
                        } else {
                            best_pair_insertion_with_info::<false>(
                                inst,
                                &buf_intermediate,
                                info_buf2,
                                p1,
                                d1,
                            )
                        };
                        if let Some((cost_final, ep1, ed1)) = ins_b2 {
                            let delta = cost_final - orig_cost;
                            if delta < best_delta {
                                best_delta = delta;
                                rt_delta[ri] = delta;
                                rt_p1[ri] = p1;
                                rt_p2[ri] = p2;
                                rt_first_is_p1[ri] = false;
                                rt_ep_first[ri] = ep2;
                                rt_ed_first[ri] = ed2;
                                rt_ep_second[ri] = ep1;
                                rt_ed_second[ri] = ed1;
                            }
                        }
                    }
                }
            }
        }

        // Find global best across all routes
        let mut best_ri = usize::MAX;
        let mut best_delta = -1e-10;
        for (ri, &delta) in rt_delta.iter().enumerate().take(nr) {
            if delta < best_delta {
                best_delta = delta;
                best_ri = ri;
            }
        }

        if best_ri == usize::MAX {
            break;
        }

        // Apply the best move
        let p1 = rt_p1[best_ri];
        let d1 = inst.delivery_of(p1);
        let p2 = rt_p2[best_ri];
        let d2 = inst.delivery_of(p2);

        let (first_p, first_d, second_p, second_d) = if rt_first_is_p1[best_ri] {
            (p1, d1, p2, d2)
        } else {
            (p2, d2, p1, d1)
        };

        // Rebuild: remove both pairs
        buf_without_both.clear();
        for &v in &sol.routes[best_ri] {
            if v != p1 && v != d1 && v != p2 && v != d2 {
                buf_without_both.push(v);
            }
        }

        // Insert first pair
        build_route_with_pair_buf(
            &buf_without_both,
            first_p,
            rt_ep_first[best_ri],
            first_d,
            rt_ed_first[best_ri],
            &mut buf_intermediate,
        );

        // Insert second pair → final route
        sol.routes[best_ri] = build_route_with_pair(
            &buf_intermediate,
            second_p,
            rt_ep_second[best_ri],
            second_d,
            rt_ed_second[best_ri],
        );
        infos[best_ri].compute(inst, &sol.routes[best_ri]);
        *move_ctr += 1;
        route_epoch[best_ri] = *move_ctr;

        // Only the modified route needs rescan
        rt_valid[best_ri] = false;
        any_improved = true;
    }

    // Update route_clean for all scanned (valid) routes
    for ri in 0..nr {
        if rt_valid[ri] {
            route_clean[ri] = *move_ctr;
        }
    }

    any_improved
}

/// SWAP* (Vidal 2022): Inter-route pair exchange with free reinsertion.
/// Best improvement with precomputed screening. Each call finds and applies
/// the single best swap across all route pairs. Unlike M3 (first improvement),
/// SWAP* evaluates all candidates with screening to find moves M3 misses.
#[allow(clippy::similar_names, clippy::too_many_lines)]
fn apply_swap_star(
    inst: &Instance,
    sol: &mut Solution,
    infos: &mut Vec<RouteInfo>,
    dirty: &mut bool,
    route_epoch: &mut [u32],
    clean_epoch: u32,
    move_ctr: &mut u32,
    penalty: &PenaltyState,
) -> bool {
    ensure_infos(inst, &sol.routes, infos, dirty);

    let nr = sol.routes.len();
    if nr < 2 {
        return false;
    }

    // Skip if no route has been modified since last full scan
    if !route_epoch.iter().take(nr).any(|&e| e > clean_epoch) {
        return false;
    }

    // --- Phase 1: Collect pairs and precompute removal deltas (O(1) per pair) ---
    let mut pair_pickups: Vec<usize> = Vec::new();
    let mut pair_routes: Vec<usize> = Vec::new();
    let mut removal_deltas: Vec<f64> = Vec::new();

    for (ri, (route, info)) in sol.routes.iter().zip(infos.iter()).enumerate().take(nr) {
        for pos in 1..route.len() - 1 {
            let p = route[pos];
            if !inst.is_pickup(p) {
                continue;
            }
            let d = inst.delivery_of(p);

            // Find delivery position (delivery always after pickup in a valid route)
            let mut pos_d = 0;
            for (i, &node) in route.iter().enumerate().take(route.len() - 1).skip(pos + 1) {
                if node == d {
                    pos_d = i;
                    break;
                }
            }

            // Compute removal delta using edge_dists (O(1))
            let removal_delta = if pos_d == pos + 1 {
                // Adjacent: a->p->d->e becomes a->e
                inst.dist(route[pos - 1], route[pos_d + 1])
                    - info.edge_dists[pos - 1]
                    - info.edge_dists[pos]
                    - info.edge_dists[pos_d]
            } else {
                // Non-adjacent: a->p->b ... c->d->e becomes a->b ... c->e
                (inst.dist(route[pos - 1], route[pos + 1])
                    - info.edge_dists[pos - 1]
                    - info.edge_dists[pos])
                    + (inst.dist(route[pos_d - 1], route[pos_d + 1])
                        - info.edge_dists[pos_d - 1]
                        - info.edge_dists[pos_d])
            };

            pair_pickups.push(p);
            pair_routes.push(ri);
            removal_deltas.push(removal_delta);
        }
    }

    let np = pair_pickups.len();
    if np == 0 {
        return false;
    }

    // Precompute conflict masks for all full routes (O(8) lookup vs O(route_len))
    let stride = inst.conflict_stride;
    let mut route_masks: Vec<u64> = vec![0u64; nr * stride];
    let mut mask_buf: Vec<u64> = Vec::new();
    for ri in 0..nr {
        inst.build_route_req_mask(&sol.routes[ri], &mut mask_buf);
        route_masks[ri * stride..(ri + 1) * stride].copy_from_slice(&mask_buf);
    }

    // --- Phase 2: Precompute screening insertions (each pair into each full route) ---
    let mut ins_screen: Vec<f64> = vec![f64::MAX; np * nr];

    for pi in 0..np {
        let p = pair_pickups[pi];
        let d = inst.delivery_of(p);
        let own_ri = pair_routes[pi];
        for ri in 0..nr {
            if ri == own_ri {
                continue;
            }
            // Conflict pre-filter: skip if p conflicts with any request in target route
            if inst.has_conflict_with_mask(&route_masks[ri * stride..(ri + 1) * stride], p) {
                continue;
            }
            if penalty.is_active() {
                if let Some((pc, _, _)) = best_pair_insertion_penalized(
                    inst,
                    &sol.routes[ri],
                    &infos[ri],
                    p,
                    d,
                    penalty.penalty_cap,
                    penalty.penalty_tw,
                ) {
                    let old_pc = penalty.penalized_cost(&infos[ri]);
                    ins_screen[pi * nr + ri] = pc - old_pc;
                }
            } else if let Some((dist, _, _)) =
                best_pair_insertion_with_info::<true>(inst, &sol.routes[ri], &infos[ri], p, d)
            {
                ins_screen[pi * nr + ri] = dist - infos[ri].distance;
            }
        }
    }

    // --- Phase 3: Best-improvement search across all route pairs ---
    let mut best_delta = -1e-10;
    let mut best_is_swap = true;
    let mut best_p1_idx: usize = usize::MAX;
    let mut best_p2_idx: usize = usize::MAX;
    let mut best_ep1: usize = 0;
    let mut best_ed1: usize = 0;
    let mut best_ep2: usize = 0;
    let mut best_ed2: usize = 0;
    // Relocate move tracking (competes with swaps via combined delta)
    let mut best_reloc_pi: usize = usize::MAX;
    let mut best_reloc_dst: usize = usize::MAX;

    // Reusable buffers for route-without-pair computation
    let mut r1_bufs: Vec<Vec<usize>> = Vec::new();
    let mut r1_info_bufs: Vec<RouteInfo> = Vec::new();
    let mut r2_bufs: Vec<Vec<usize>> = Vec::new();
    let mut r2_info_bufs: Vec<RouteInfo> = Vec::new();

    // Pair indices per route for fast iteration
    let mut route_pair_indices: Vec<Vec<usize>> = vec![Vec::new(); nr];
    for pi in 0..np {
        route_pair_indices[pair_routes[pi]].push(pi);
    }

    for ri in 0..nr {
        if route_pair_indices[ri].is_empty() {
            continue;
        }
        for rj in (ri + 1)..nr {
            if route_pair_indices[rj].is_empty() {
                continue;
            }
            // Skip if neither route changed since last full scan
            if route_epoch[ri] <= clean_epoch && route_epoch[rj] <= clean_epoch {
                continue;
            }

            // Skip route pairs with non-overlapping circle sectors
            if inst.polar_angle.is_some()
                && !CircleSector::overlap(&infos[ri].sector, &infos[rj].sector)
            {
                continue;
            }

            let ri_pairs = &route_pair_indices[ri];
            let rj_pairs = &route_pair_indices[rj];

            // Precompute ri_without_p1 + RouteInfo for all p1 in ri
            while r1_bufs.len() < ri_pairs.len() {
                r1_bufs.push(Vec::new());
                r1_info_bufs.push(RouteInfo::new());
            }
            for (k, &pi) in ri_pairs.iter().enumerate() {
                let p = pair_pickups[pi];
                let d = inst.delivery_of(p);
                compute_info_without_pair(
                    inst,
                    &sol.routes[ri],
                    p,
                    d,
                    &mut r1_bufs[k],
                    &mut r1_info_bufs[k],
                );
                if r1_bufs[k].len() >= 2 && penalty.is_active() {
                    r1_info_bufs[k].compute_tw_sector(inst, &r1_bufs[k], r1_bufs[k].len());
                }
            }

            // Precompute rj_without_p2 + RouteInfo for all p2 in rj
            while r2_bufs.len() < rj_pairs.len() {
                r2_bufs.push(Vec::new());
                r2_info_bufs.push(RouteInfo::new());
            }
            for (k, &pi) in rj_pairs.iter().enumerate() {
                let p = pair_pickups[pi];
                let d = inst.delivery_of(p);
                compute_info_without_pair(
                    inst,
                    &sol.routes[rj],
                    p,
                    d,
                    &mut r2_bufs[k],
                    &mut r2_info_bufs[k],
                );
                if r2_bufs[k].len() >= 2 && penalty.is_active() {
                    r2_info_bufs[k].compute_tw_sector(inst, &r2_bufs[k], r2_bufs[k].len());
                }
            }

            // Try all (p1 from ri, p2 from rj) combinations
            for (k1, &pi1) in ri_pairs.iter().enumerate() {
                if r1_bufs[k1].len() < 2 {
                    continue;
                }

                for (k2, &pi2) in rj_pairs.iter().enumerate() {
                    if r2_bufs[k2].len() < 2 {
                        continue;
                    }

                    // Phase 3a: Approximate screen
                    let approx = removal_deltas[pi1]
                        + removal_deltas[pi2]
                        + ins_screen[pi1 * nr + rj]
                        + ins_screen[pi2 * nr + ri];
                    if approx >= best_delta {
                        continue;
                    }

                    // Phase 3b: Exact insertion of p1 into rj_without_p2
                    let p1 = pair_pickups[pi1];
                    let d1 = inst.delivery_of(p1);
                    // Conflict pre-filter (2-level: full mask, then excluding p2)
                    let p2_for_exclude = pair_pickups[pi2];
                    if inst.has_conflict_with_mask_excluding(
                        &route_masks[rj * stride..(rj + 1) * stride],
                        p1,
                        p2_for_exclude,
                    ) {
                        continue;
                    }
                    let (pc_rj_new, ep1, ed1) = if penalty.is_active() {
                        match best_pair_insertion_penalized(
                            inst,
                            &r2_bufs[k2],
                            &r2_info_bufs[k2],
                            p1,
                            d1,
                            penalty.penalty_cap,
                            penalty.penalty_tw,
                        ) {
                            Some((pc, e1, e2)) => (pc, e1, e2),
                            None => continue,
                        }
                    } else {
                        match best_pair_insertion_with_info::<true>(
                            inst,
                            &r2_bufs[k2],
                            &r2_info_bufs[k2],
                            p1,
                            d1,
                        ) {
                            Some((d, e1, e2)) => (d, e1, e2),
                            None => continue,
                        }
                    };
                    let old_pc_rj = penalty.penalized_cost(&infos[rj]);
                    let delta_rj = pc_rj_new - old_pc_rj;

                    // Early termination: delta_rj + approximate ri side
                    if delta_rj + removal_deltas[pi1] + ins_screen[pi2 * nr + ri] >= best_delta {
                        continue;
                    }

                    // Phase 3c: Exact insertion of p2 into ri_without_p1
                    let p2 = pair_pickups[pi2];
                    let d2 = inst.delivery_of(p2);
                    // Conflict pre-filter (2-level: full mask, then excluding p1)
                    if inst.has_conflict_with_mask_excluding(
                        &route_masks[ri * stride..(ri + 1) * stride],
                        p2,
                        p1,
                    ) {
                        continue;
                    }
                    let (pc_ri_new, ep2, ed2) = if penalty.is_active() {
                        match best_pair_insertion_penalized(
                            inst,
                            &r1_bufs[k1],
                            &r1_info_bufs[k1],
                            p2,
                            d2,
                            penalty.penalty_cap,
                            penalty.penalty_tw,
                        ) {
                            Some((pc, e1, e2)) => (pc, e1, e2),
                            None => continue,
                        }
                    } else {
                        match best_pair_insertion_with_info::<true>(
                            inst,
                            &r1_bufs[k1],
                            &r1_info_bufs[k1],
                            p2,
                            d2,
                        ) {
                            Some((d, e1, e2)) => (d, e1, e2),
                            None => continue,
                        }
                    };
                    let old_pc_ri = penalty.penalized_cost(&infos[ri]);
                    let delta_ri = pc_ri_new - old_pc_ri;

                    // Phase 3d: Track global best
                    let total_delta = delta_ri + delta_rj;
                    if total_delta < best_delta {
                        best_delta = total_delta;
                        best_is_swap = true;
                        best_p1_idx = pi1;
                        best_p2_idx = pi2;
                        best_ep1 = ep1;
                        best_ed1 = ed1;
                        best_ep2 = ep2;
                        best_ed2 = ed2;
                    }
                }
            }

            // --- Relocate: pairs from ri → full rj ---
            for (k1, &pi1) in ri_pairs.iter().enumerate() {
                let ins_delta = ins_screen[pi1 * nr + rj];
                if ins_delta >= f64::MAX / 2.0 {
                    continue;
                }

                let empties = r1_bufs[k1].len() <= 2;
                let old_pc_ri = penalty.penalized_cost(&infos[ri]);
                let exact_removal = if empties {
                    -old_pc_ri
                } else {
                    let new_pc_ri = penalty.penalized_cost(&r1_info_bufs[k1]);
                    new_pc_ri - old_pc_ri
                };

                let vehicle_bonus = if empties { -1e9 } else { 0.0 };
                let total = exact_removal + ins_delta + vehicle_bonus;
                if total < best_delta {
                    best_delta = total;
                    best_is_swap = false;
                    best_reloc_pi = pi1;
                    best_reloc_dst = rj;
                }
            }

            // --- Relocate: pairs from rj → full ri ---
            for (k2, &pi2) in rj_pairs.iter().enumerate() {
                let ins_delta = ins_screen[pi2 * nr + ri];
                if ins_delta >= f64::MAX / 2.0 {
                    continue;
                }

                let empties = r2_bufs[k2].len() <= 2;
                let old_pc_rj = penalty.penalized_cost(&infos[rj]);
                let exact_removal = if empties {
                    -old_pc_rj
                } else {
                    let new_pc_rj = penalty.penalized_cost(&r2_info_bufs[k2]);
                    new_pc_rj - old_pc_rj
                };

                let vehicle_bonus = if empties { -1e9 } else { 0.0 };
                let total = exact_removal + ins_delta + vehicle_bonus;
                if total < best_delta {
                    best_delta = total;
                    best_is_swap = false;
                    best_reloc_pi = pi2;
                    best_reloc_dst = ri;
                }
            }
        }
    }

    // --- Phase 4: Apply the best move (swap or relocate) ---
    if best_is_swap {
        if best_p1_idx == usize::MAX {
            return false;
        }

        let ri = pair_routes[best_p1_idx];
        let rj = pair_routes[best_p2_idx];
        let p1 = pair_pickups[best_p1_idx];
        let d1 = inst.delivery_of(p1);
        let p2 = pair_pickups[best_p2_idx];
        let d2 = inst.delivery_of(p2);

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

        sol.routes[ri] = build_route_with_pair(&r1_without, p2, best_ep2, d2, best_ed2);
        sol.routes[rj] = build_route_with_pair(&r2_without, p1, best_ep1, d1, best_ed1);

        infos[ri].compute(inst, &sol.routes[ri]);
        infos[rj].compute(inst, &sol.routes[rj]);
        *move_ctr += 1;
        route_epoch[ri] = *move_ctr;
        *move_ctr += 1;
        route_epoch[rj] = *move_ctr;
    } else {
        if best_reloc_pi == usize::MAX {
            return false;
        }

        let pi = best_reloc_pi;
        let dst = best_reloc_dst;
        let p = pair_pickups[pi];
        let d = inst.delivery_of(p);
        let src = pair_routes[pi];

        // Recompute exact insertion into dst (need positions ep, ed)
        let (ep, ed) = if penalty.is_active() {
            let (_, ep, ed) = best_pair_insertion_penalized(
                inst,
                &sol.routes[dst],
                &infos[dst],
                p,
                d,
                penalty.penalty_cap,
                penalty.penalty_tw,
            )
            .expect("relocate screened as feasible");
            (ep, ed)
        } else {
            let (_, ep, ed) =
                best_pair_insertion_with_info::<true>(inst, &sol.routes[dst], &infos[dst], p, d)
                    .expect("relocate screened as feasible");
            (ep, ed)
        };

        let new_dst = build_route_with_pair(&sol.routes[dst], p, ep, d, ed);
        let src_without = route_without_pair(&sol.routes[src], p, d);

        sol.routes[dst] = new_dst;
        if src_without.len() <= 2 {
            sol.routes.remove(src);
            *dirty = true;
        } else {
            sol.routes[src] = src_without;
            infos[src].compute(inst, &sol.routes[src]);
            infos[dst].compute(inst, &sol.routes[dst]);
            *move_ctr += 1;
            route_epoch[src] = *move_ctr;
            *move_ctr += 1;
            route_epoch[dst] = *move_ctr;
        }
    }

    true
}