pdp_lns 0.1.1

use crate::instance::Instance;
use crate::lns;
use crate::solution::{RouteInfo, Solution};
use rand::prelude::*;
use std::collections::{HashMap, HashSet};

/// Number of random subset-pairs generated before ranking by LCS.
const CREX_NTOTAL: usize = 40;
/// Number of top-ranked subset-pairs used to construct offspring.
const CREX_NCROSS: usize = 20;
/// Max number of local-improvement iterations in MAKE_NEIGHBORHOOD loop.
const CREX_MAX_NEIGHBOR_ITERS: usize = 16;
/// Max random (a_cand, b_cand) attempts per grow step before giving up.
const GROW_ATTEMPTS_PER_STEP: usize = 5;
/// Regret-k parameter for offspring repair.
const CREX_REPAIR_K: usize = 4;

#[derive(Clone, Copy)]
struct CrossoverParams {
    ntotal: usize,
    ncross: usize,
    max_neighbor_iters: usize,
}

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

fn load_crossover_params() -> CrossoverParams {
    let ntotal = env_parse_usize("HP_CREX_NTOTAL", CREX_NTOTAL).max(1);
    let ncross = env_parse_usize("HP_CREX_NCROSS", CREX_NCROSS)
        .max(1)
        .min(ntotal);
    let max_neighbor_iters =
        env_parse_usize("HP_CREX_MAX_NEIGHBOR_ITERS", CREX_MAX_NEIGHBOR_ITERS).max(1);

    CrossoverParams {
        ntotal,
        ncross,
        max_neighbor_iters,
    }
}

#[derive(Clone)]
struct SubsetPair {
    sa: Vec<usize>,
    sb: Vec<usize>,
    lcs: usize,
}

#[inline]
fn random_subset_indices(len: usize, k: usize, rng: &mut impl Rng) -> Vec<usize> {
    let mut idx: Vec<usize> = (0..len).collect();
    idx.shuffle(rng);
    idx.truncate(k);
    idx.sort_unstable();
    idx
}

fn route_pickups(inst: &Instance, route: &[usize]) -> Vec<usize> {
    route
        .iter()
        .copied()
        .filter(|&v| v != 0 && inst.is_pickup(v))
        .collect()
}

fn subset_arc_sequence(routes: &[Vec<usize>], subset: &[usize], node_base: usize) -> Vec<u64> {
    let mut seq = Vec::new();
    for &ri in subset {
        let route = &routes[ri];
        for w in route.windows(2) {
            let a = w[0] as u64;
            let b = w[1] as u64;
            seq.push(a * node_base as u64 + b);
        }
    }
    seq
}

fn lcs_len(a: &[u64], b: &[u64]) -> usize {
    if a.is_empty() || b.is_empty() {
        return 0;
    }

    let mut prev = vec![0usize; b.len() + 1];
    let mut curr = vec![0usize; b.len() + 1];
    for &x in a {
        curr[0] = 0;
        for (j, &y) in b.iter().enumerate() {
            curr[j + 1] = if x == y {
                prev[j] + 1
            } else {
                prev[j + 1].max(curr[j])
            };
        }
        std::mem::swap(&mut prev, &mut curr);
    }
    prev[b.len()]
}

fn subset_pair_lcs(
    parent_a: &Solution,
    parent_b: &Solution,
    sa: &[usize],
    sb: &[usize],
    node_base: usize,
) -> usize {
    let xa = subset_arc_sequence(&parent_a.routes, sa, node_base);
    let xb = subset_arc_sequence(&parent_b.routes, sb, node_base);
    lcs_len(&xa, &xb)
}

fn improve_subset_pair(
    parent_a: &Solution,
    parent_b: &Solution,
    sa: &mut Vec<usize>,
    sb: &mut Vec<usize>,
    node_base: usize,
    max_neighbor_iters: usize,
    rng: &mut impl Rng,
) -> usize {
    let mut lcs = subset_pair_lcs(parent_a, parent_b, sa, sb, node_base);
    let ra = parent_a.routes.len();
    let rb = parent_b.routes.len();

    // Phase 1: Grow-based improvement (Blocho & Nalepa 2017 MakeNeighborhood).
    // Add a route to both S_A and S_B simultaneously while LCS improves.
    for _ in 0..max_neighbor_iters {
        if sa.len() >= ra || sb.len() >= rb {
            break;
        }

        let mut in_sa = vec![false; ra];
        for &ri in sa.iter() {
            in_sa[ri] = true;
        }
        let mut in_sb = vec![false; rb];
        for &ri in sb.iter() {
            in_sb[ri] = true;
        }

        let mut a_cands: Vec<usize> = (0..ra).filter(|&i| !in_sa[i]).collect();
        let mut b_cands: Vec<usize> = (0..rb).filter(|&i| !in_sb[i]).collect();
        if a_cands.is_empty() || b_cands.is_empty() {
            break;
        }
        a_cands.shuffle(rng);
        b_cands.shuffle(rng);

        let mut grew = false;
        let mut attempts = 0;
        'grow: for &ac in &a_cands {
            for &bc in &b_cands {
                attempts += 1;
                let mut sa2 = sa.clone();
                sa2.push(ac);
                sa2.sort_unstable();
                let mut sb2 = sb.clone();
                sb2.push(bc);
                sb2.sort_unstable();
                let l = subset_pair_lcs(parent_a, parent_b, &sa2, &sb2, node_base);
                if l > lcs {
                    *sa = sa2;
                    *sb = sb2;
                    lcs = l;
                    grew = true;
                    break 'grow;
                }
                if attempts >= GROW_ATTEMPTS_PER_STEP {
                    break 'grow;
                }
            }
        }
        if !grew {
            break;
        }
    }

    // Phase 2: Swap-based refinement — try replacing individual routes
    // within fixed-size subsets to maximize LCS.
    for _ in 0..max_neighbor_iters {
        let mut best_lcs = lcs;
        let mut best_sa = sa.clone();
        let mut best_sb = sb.clone();

        let mut in_sa = vec![false; ra];
        for &ri in sa.iter() {
            in_sa[ri] = true;
        }

        for pos in 0..sa.len() {
            for (cand, &is_in_sa) in in_sa.iter().enumerate() {
                if is_in_sa {
                    continue;
                }
                let mut sa2 = sa.clone();
                sa2[pos] = cand;
                sa2.sort_unstable();
                let l = subset_pair_lcs(parent_a, parent_b, &sa2, sb, node_base);
                if l > best_lcs {
                    best_lcs = l;
                    best_sa = sa2;
                    best_sb.clone_from(sb);
                }
            }
        }

        let mut in_sb = vec![false; rb];
        for &ri in sb.iter() {
            in_sb[ri] = true;
        }

        for pos in 0..sb.len() {
            for (cand, &is_in_sb) in in_sb.iter().enumerate() {
                if is_in_sb {
                    continue;
                }
                let mut sb2 = sb.clone();
                sb2[pos] = cand;
                sb2.sort_unstable();
                let l = subset_pair_lcs(parent_a, parent_b, sa, &sb2, node_base);
                if l > best_lcs {
                    best_lcs = l;
                    best_sa.clone_from(sa);
                    best_sb = sb2;
                }
            }
        }

        if best_lcs > lcs {
            *sa = best_sa;
            *sb = best_sb;
            lcs = best_lcs;
        } else {
            break;
        }
    }

    lcs
}

fn filter_route_with_mask(
    inst: &Instance,
    route: &[usize],
    keep_pickup: &[bool],
    served_pickup: &mut [bool],
) -> Vec<usize> {
    let mut out = Vec::with_capacity(route.len());
    let mut keep_local = vec![false; inst.n + 1];
    out.push(0);

    for &v in &route[1..route.len() - 1] {
        if inst.is_pickup(v) {
            if keep_pickup[v] && !served_pickup[v] {
                served_pickup[v] = true;
                keep_local[v] = true;
                out.push(v);
            }
        } else {
            let p = inst.pickup_of(v);
            if p != 0 && keep_local[p] {
                out.push(v);
                keep_local[p] = false;
            }
        }
    }

    out.push(0);
    out
}

fn build_offspring(
    inst: &Instance,
    parent_a: &Solution,
    parent_b: &Solution,
    a_route_pickups: &[Vec<usize>],
    b_route_pickups: &[Vec<usize>],
    sa: &[usize],
    sb: &[usize],
    include_all_sb: bool,
    rng: &mut impl Rng,
) -> Solution {
    let mut a_removed = vec![false; inst.n + 1];
    for &ri in sa {
        for &p in &a_route_pickups[ri] {
            a_removed[p] = true;
        }
    }

    let mut b_subset = vec![false; inst.n + 1];
    for &ri in sb {
        for &p in &b_route_pickups[ri] {
            b_subset[p] = true;
        }
    }

    let mut ejected = vec![false; inst.n + 1];
    for &p in &inst.pickups {
        ejected[p] = b_subset[p] && !a_removed[p];
    }

    let mut keep_from_a = vec![true; inst.n + 1];
    for &p in &inst.pickups {
        if ejected[p] {
            keep_from_a[p] = false;
        }
    }

    let mut in_sa = vec![false; parent_a.routes.len()];
    for &ri in sa {
        in_sa[ri] = true;
    }

    let mut served = vec![false; inst.n + 1];
    let mut child_routes: Vec<Vec<usize>> = Vec::new();

    for (ri, route) in parent_a.routes.iter().enumerate() {
        if in_sa[ri] {
            continue;
        }
        let r = filter_route_with_mask(inst, route, &keep_from_a, &mut served);
        if r.len() > 2 {
            child_routes.push(r);
        }
    }

    let keep_all = vec![true; inst.n + 1];
    for &ri in sb {
        if !include_all_sb && !b_route_pickups[ri].iter().any(|&p| ejected[p]) {
            continue;
        }
        let r = filter_route_with_mask(inst, &parent_b.routes[ri], &keep_all, &mut served);
        if r.len() > 2 {
            child_routes.push(r);
        }
    }

    let mut child = Solution {
        routes: child_routes,
        unassigned: inst
            .pickups
            .iter()
            .copied()
            .filter(|&p| !served[p])
            .collect(),
    };

    let mut infos: Vec<RouteInfo> = Vec::new();
    lns::regret_insertion(
        inst,
        &mut child,
        CREX_REPAIR_K,
        0.0,
        true,
        rng,
        &mut infos,
        0,
        false,
    );

    child
}

/// Full LCS-SREX crossover (Algorithm 2 style) adapted for this solver.
///
/// Returns the best child among `Ncross` offspring candidates, each constructed
/// from one of the top-ranked subset pairs by LCS value.
pub fn lcs_srex_crossover(
    inst: &Instance,
    parent_a: &Solution,
    parent_b: &Solution,
    rng: &mut impl Rng,
) -> Option<Solution> {
    if parent_a.routes.is_empty() || parent_b.routes.is_empty() {
        return None;
    }

    let max_k = parent_a.routes.len().min(parent_b.routes.len());
    if max_k == 0 {
        return None;
    }

    let params = load_crossover_params();

    let node_base = inst.n + 1;
    let a_route_pickups: Vec<Vec<usize>> = parent_a
        .routes
        .iter()
        .map(|r| route_pickups(inst, r))
        .collect();
    let b_route_pickups: Vec<Vec<usize>> = parent_b
        .routes
        .iter()
        .map(|r| route_pickups(inst, r))
        .collect();

    let mut subset_pairs: Vec<SubsetPair> = Vec::with_capacity(params.ntotal);
    for _ in 0..params.ntotal {
        let k = rng.random_range(1..=max_k);
        let mut sa = random_subset_indices(parent_a.routes.len(), k, rng);
        let mut sb = random_subset_indices(parent_b.routes.len(), k, rng);
        let lcs = improve_subset_pair(
            parent_a,
            parent_b,
            &mut sa,
            &mut sb,
            node_base,
            params.max_neighbor_iters,
            rng,
        );
        subset_pairs.push(SubsetPair { sa, sb, lcs });
    }

    subset_pairs.sort_by_key(|b| std::cmp::Reverse(b.lcs));

    let mut seen: HashSet<(Vec<usize>, Vec<usize>)> = HashSet::new();
    let mut best_children: Vec<SubsetPair> = Vec::new();
    for cand in subset_pairs {
        let key = (cand.sa.clone(), cand.sb.clone());
        if seen.insert(key) {
            best_children.push(cand);
            if best_children.len() >= params.ncross {
                break;
            }
        }
    }

    let mut best_child: Option<Solution> = None;
    let mut best_cost = f64::MAX;

    for cand in best_children {
        let child1 = build_offspring(
            inst,
            parent_a,
            parent_b,
            &a_route_pickups,
            &b_route_pickups,
            &cand.sa,
            &cand.sb,
            true,
            rng,
        );
        let child2 = build_offspring(
            inst,
            parent_a,
            parent_b,
            &a_route_pickups,
            &b_route_pickups,
            &cand.sa,
            &cand.sb,
            false,
            rng,
        );

        let (candidate, cand_cost) = if child1.cost(inst) <= child2.cost(inst) {
            let c = child1.cost(inst);
            (child1, c)
        } else {
            let c = child2.cost(inst);
            (child2, c)
        };

        if candidate.is_feasible(inst) && cand_cost < best_cost {
            best_cost = cand_cost;
            best_child = Some(candidate);
        }
    }

    best_child
}

/// SREX-like perturbation: swap a small random subset of routes from `parent_b`
/// into `parent_a`, then repair unassigned requests with regret insertion.
///
/// `swap_routes` is clamped to `[1, parent_b.routes.len()]`.
pub fn srex_route_exchange(
    inst: &Instance,
    parent_a: &Solution,
    parent_b: &Solution,
    swap_routes: usize,
    rng: &mut impl Rng,
) -> Option<Solution> {
    if parent_a.routes.is_empty() || parent_b.routes.is_empty() {
        return None;
    }

    let k = swap_routes.clamp(1, parent_b.routes.len());
    let sb = random_subset_indices(parent_b.routes.len(), k, rng);
    let sa: Vec<usize> = Vec::new();

    let a_route_pickups: Vec<Vec<usize>> = parent_a
        .routes
        .iter()
        .map(|r| route_pickups(inst, r))
        .collect();
    let b_route_pickups: Vec<Vec<usize>> = parent_b
        .routes
        .iter()
        .map(|r| route_pickups(inst, r))
        .collect();

    let child1 = build_offspring(
        inst,
        parent_a,
        parent_b,
        &a_route_pickups,
        &b_route_pickups,
        &sa,
        &sb,
        true,
        rng,
    );
    let child2 = build_offspring(
        inst,
        parent_a,
        parent_b,
        &a_route_pickups,
        &b_route_pickups,
        &sa,
        &sb,
        false,
        rng,
    );

    let candidate = if child1.cost(inst) <= child2.cost(inst) {
        child1
    } else {
        child2
    };

    if candidate.is_feasible(inst) {
        Some(candidate)
    } else {
        None
    }
}

/// Number of offspring to generate per EAX crossover call.
const EAX_NUM_OFFSPRING: usize = 8;

/// Score a tour direction by PD violations, arc mismatches, and distance.
/// `parent_arcs` encodes directed arcs as `u * stride + v`.
/// Returns (num_pd_violations, num_arc_mismatches, total_distance).
fn tour_direction_score(
    inst: &Instance,
    inner: &[usize],
    parent_arcs: &HashSet<u64>,
    stride: u64,
) -> (usize, usize, f64) {
    if inner.is_empty() {
        return (0, 0, 0.0);
    }
    let mut dist = inst.dist(0, inner[0]) + inst.dist(inner[inner.len() - 1], 0);
    let mut mismatches = 0usize;
    // Check depot→first arc
    if !parent_arcs.contains(&(inner[0] as u64)) {
        mismatches += 1;
    }
    // Check last→depot arc
    if !parent_arcs.contains(&(inner[inner.len() - 1] as u64 * stride)) {
        mismatches += 1;
    }
    for w in inner.windows(2) {
        dist += inst.dist(w[0], w[1]);
        if !parent_arcs.contains(&(w[0] as u64 * stride + w[1] as u64)) {
            mismatches += 1;
        }
    }
    let mut pickup_seen = vec![false; inst.n + 1];
    let mut violations = 0usize;
    for &v in inner {
        if inst.is_pickup(v) {
            pickup_seen[v] = true;
        } else {
            let p = inst.pickup_of(v);
            if p != 0 && !pickup_seen[p] {
                violations += 1;
            }
        }
    }
    (violations, mismatches, dist)
}

/// EAX crossover for PDPTW (Edge Assembly Crossover).
///
/// Adapted from Nagata & Bräysy (2009) for pickup-and-delivery:
/// 1. Build undirected edge sets from both parents
/// 2. Find AB-cycles in the symmetric difference
/// 3. Generate offspring by selecting random AB-cycles
/// 4. Extract and orient tours, repair PD constraints
pub fn eax_crossover(
    inst: &Instance,
    parent_a: &Solution,
    parent_b: &Solution,
    rng: &mut impl Rng,
) -> Option<Solution> {
    if parent_a.routes.is_empty() || parent_b.routes.is_empty() {
        return None;
    }
    let n = inst.n;

    // --- Step 1: Build undirected edge multisets ---
    // Key: (min(u,v), max(u,v)), Value: [count_a, count_b]
    let mut edge_counts: HashMap<(usize, usize), [u32; 2]> = HashMap::new();
    for route in &parent_a.routes {
        for w in route.windows(2) {
            let key = (w[0].min(w[1]), w[0].max(w[1]));
            edge_counts.entry(key).or_default()[0] += 1;
        }
    }
    for route in &parent_b.routes {
        for w in route.windows(2) {
            let key = (w[0].min(w[1]), w[0].max(w[1]));
            edge_counts.entry(key).or_default()[1] += 1;
        }
    }

    // --- Step 2: Classify edges into common and symmetric difference ---
    let mut common_edges: Vec<(usize, usize)> = Vec::new();
    // Sym-diff edges stored as parallel arrays for cache efficiency
    let mut sd_u: Vec<usize> = Vec::new();
    let mut sd_v: Vec<usize> = Vec::new();
    let mut sd_is_a: Vec<bool> = Vec::new();

    for (&(u, v), &[ca, cb]) in &edge_counts {
        let common = ca.min(cb);
        for _ in 0..common {
            common_edges.push((u, v));
        }
        for _ in 0..(ca - common) {
            sd_u.push(u);
            sd_v.push(v);
            sd_is_a.push(true);
        }
        for _ in 0..(cb - common) {
            sd_u.push(u);
            sd_v.push(v);
            sd_is_a.push(false);
        }
    }
    if sd_u.is_empty() {
        return None; // parents have identical edge sets
    }

    // --- Step 3: Balance depot degree in symmetric difference ---
    // Each parent's depot degree = 2 * num_routes. If route counts differ,
    // depot has unequal A/B degree in sym-diff. Add self-loops to balance.
    let mut depot_deg_a = 0u32;
    let mut depot_deg_b = 0u32;
    for i in 0..sd_u.len() {
        let deg = u32::from(sd_u[i] == 0) + u32::from(sd_v[i] == 0);
        if sd_is_a[i] {
            depot_deg_a += deg;
        } else {
            depot_deg_b += deg;
        }
    }
    // Each self-loop adds 2 to its parent's depot degree
    if depot_deg_a < depot_deg_b {
        for _ in 0..(depot_deg_b - depot_deg_a) / 2 {
            sd_u.push(0);
            sd_v.push(0);
            sd_is_a.push(true);
        }
    } else if depot_deg_b < depot_deg_a {
        for _ in 0..(depot_deg_a - depot_deg_b) / 2 {
            sd_u.push(0);
            sd_v.push(0);
            sd_is_a.push(false);
        }
    }
    let ne = sd_u.len();

    // --- Step 4: Build adjacency and trace AB-cycles ---
    let mut used = vec![false; ne];
    // adj[node] = [(other_endpoint, edge_index)]
    let mut adj: Vec<Vec<(usize, usize)>> = vec![Vec::new(); n + 1];
    for i in 0..ne {
        let u = sd_u[i];
        let v = sd_v[i];
        adj[u].push((v, i));
        if u == v {
            // Self-loop: appears twice in adjacency for correct degree
            adj[u].push((u, i));
        } else {
            adj[v].push((u, i));
        }
    }

    let mut cycles: Vec<Vec<usize>> = Vec::new();
    for start_ei in 0..ne {
        if used[start_ei] {
            continue;
        }

        used[start_ei] = true;
        let start_node = sd_u[start_ei];
        let mut current = if sd_u[start_ei] == sd_v[start_ei] {
            sd_u[start_ei]
        } else {
            sd_v[start_ei]
        };
        let mut cycle = vec![start_ei];
        let mut next_is_a = !sd_is_a[start_ei];

        let mut closed = false;
        loop {
            if current == start_node && cycle.len() >= 2 && cycle.len() % 2 == 0 {
                closed = true;
                break;
            }
            // Limit cycle length to prevent infinite loops on malformed graphs
            if cycle.len() >= ne {
                break;
            }

            let mut found = false;
            for &(neighbor, ei) in &adj[current] {
                if !used[ei] && sd_is_a[ei] == next_is_a {
                    used[ei] = true;
                    cycle.push(ei);
                    current = neighbor;
                    next_is_a = !next_is_a;
                    found = true;
                    break;
                }
            }
            if !found {
                break;
            }
        }

        if closed && cycle.len() >= 2 {
            cycles.push(cycle);
        } else {
            // Failed to close — unmark edges so they may be used elsewhere
            for &ei in &cycle {
                used[ei] = false;
            }
        }
    }

    if cycles.is_empty() {
        return None;
    }

    // Build directed arc set from both parents for tour orientation scoring.
    let stride = (n + 1) as u64;
    let mut parent_arcs: HashSet<u64> = HashSet::new();
    for route in &parent_a.routes {
        for w in route.windows(2) {
            parent_arcs.insert(w[0] as u64 * stride + w[1] as u64);
        }
    }
    for route in &parent_b.routes {
        for w in route.windows(2) {
            parent_arcs.insert(w[0] as u64 * stride + w[1] as u64);
        }
    }

    // --- Step 5: Generate offspring ---
    let num_offspring = (cycles.len() * 2).min(EAX_NUM_OFFSPRING);
    let mut best_child: Option<Solution> = None;
    let mut best_cost = f64::MAX;

    for oi in 0..num_offspring {
        // Select cycles: each offspring uses a different primary cycle.
        // First batch: single-cycle (small perturbation).
        // Second batch: multi-cycle (aggressive recombination).
        let primary = oi % cycles.len();
        let mut in_selected = vec![false; ne];
        for &ei in &cycles[primary] {
            in_selected[ei] = true;
        }
        if oi >= cycles.len() {
            for (ci, cycle) in cycles.iter().enumerate() {
                if ci != primary && rng.random_bool(0.3) {
                    for &ei in cycle {
                        in_selected[ei] = true;
                    }
                }
            }
        }

        // Build offspring undirected adjacency.
        // Include edge i if: it's a B-edge in a selected cycle, or an A-edge NOT in a selected cycle.
        // Equivalently: include if in_selected[i] XOR sd_is_a[i] (i.e., they differ).
        let mut off_adj: Vec<Vec<usize>> = vec![Vec::new(); n + 1];

        // Common edges (always included)
        for &(u, v) in &common_edges {
            off_adj[u].push(v);
            if u != v {
                off_adj[v].push(u);
            }
        }

        // Sym-diff edges based on selection
        for i in 0..ne {
            if sd_u[i] == 0 && sd_v[i] == 0 {
                continue; // skip self-loops
            }
            if in_selected[i] != sd_is_a[i] {
                let u = sd_u[i];
                let v = sd_v[i];
                off_adj[u].push(v);
                if u != v {
                    off_adj[v].push(u);
                }
            }
        }

        // --- Step 6: Extract tours from depot ---
        let mut off_used: Vec<Vec<bool>> = off_adj
            .iter()
            .map(|neighbors| vec![false; neighbors.len()])
            .collect();
        let mut routes: Vec<Vec<usize>> = Vec::new();
        let mut visited = vec![false; n + 1];

        loop {
            // Find unused edge at depot
            let start_pos = off_adj[0]
                .iter()
                .enumerate()
                .position(|(i, _)| !off_used[0][i]);
            let Some(start_pos) = start_pos else { break };

            let mut tour = vec![0usize];
            let mut cur = 0usize;
            let mut cur_edge_pos = start_pos;

            loop {
                off_used[cur][cur_edge_pos] = true;
                let next = off_adj[cur][cur_edge_pos];

                // Mark reverse edge as used
                if let Some(rev) = off_adj[next]
                    .iter()
                    .enumerate()
                    .position(|(i, &v)| v == cur && !off_used[next][i])
                {
                    off_used[next][rev] = true;
                }

                cur = next;
                if cur == 0 {
                    tour.push(0);
                    break;
                }
                visited[cur] = true;
                tour.push(cur);

                // Find next unused edge at cur
                if let Some(pos) = off_adj[cur]
                    .iter()
                    .enumerate()
                    .position(|(i, _)| !off_used[cur][i])
                {
                    cur_edge_pos = pos;
                } else {
                    tour.push(0); // dead end, close at depot
                    break;
                }
            }

            if tour.len() > 2 {
                routes.push(tour);
            }
        }

        // --- Step 7: Orient tours and fix PD constraints ---
        let mut child_routes: Vec<Vec<usize>> = Vec::new();
        for tour in &routes {
            let inner = &tour[1..tour.len() - 1];
            if inner.is_empty() {
                continue;
            }

            // Choose better orientation (fewer PD violations, then fewer arc mismatches, then shorter distance)
            let fwd_score = tour_direction_score(inst, inner, &parent_arcs, stride);
            let rev: Vec<usize> = inner.iter().rev().copied().collect();
            let rev_score = tour_direction_score(inst, &rev, &parent_arcs, stride);

            let chosen = if fwd_score <= rev_score {
                inner.to_vec()
            } else {
                rev
            };

            let mut route = Vec::with_capacity(chosen.len() + 2);
            route.push(0);
            route.extend_from_slice(&chosen);
            route.push(0);
            child_routes.push(route);
        }

        // Map nodes to their route and position
        let mut node_route: Vec<usize> = vec![usize::MAX; n + 1];
        let mut node_pos: Vec<usize> = vec![0; n + 1];
        for (ri, route) in child_routes.iter().enumerate() {
            for (pos, &v) in route.iter().enumerate() {
                if v != 0 {
                    node_route[v] = ri;
                    node_pos[v] = pos;
                }
            }
        }

        // Find PD pairs that violate constraints
        let mut bad_pickups: Vec<usize> = Vec::new();
        for &p in &inst.pickups {
            let d = inst.delivery_of(p);
            let rp = node_route[p];
            let rd = node_route[d];
            if rp == usize::MAX || rd == usize::MAX || rp != rd || node_pos[p] > node_pos[d] {
                bad_pickups.push(p);
            }
        }

        // Remove bad pairs from routes
        if !bad_pickups.is_empty() {
            let mut remove_nodes = vec![false; n + 1];
            for &p in &bad_pickups {
                remove_nodes[p] = true;
                remove_nodes[inst.delivery_of(p)] = true;
            }
            let mut clean_routes = Vec::new();
            for route in &child_routes {
                let clean: Vec<usize> = route
                    .iter()
                    .copied()
                    .filter(|&v| v == 0 || !remove_nodes[v])
                    .collect();
                if clean.len() > 2 {
                    clean_routes.push(clean);
                }
            }
            child_routes = clean_routes;
        }

        // Collect all unassigned pickups (bad pairs + subtour nodes)
        let mut in_routes = vec![false; n + 1];
        for route in &child_routes {
            for &v in &route[1..route.len() - 1] {
                in_routes[v] = true;
            }
        }
        let unassigned: Vec<usize> = inst
            .pickups
            .iter()
            .copied()
            .filter(|&p| !in_routes[p])
            .collect();

        let mut child = Solution {
            routes: child_routes,
            unassigned,
        };

        // Repair with regret insertion
        if !child.unassigned.is_empty() {
            let mut infos: Vec<RouteInfo> = Vec::new();
            lns::regret_insertion(
                inst,
                &mut child,
                CREX_REPAIR_K,
                0.0,
                true,
                rng,
                &mut infos,
                0,
                false,
            );
        }

        // Evaluate
        let cost = child.cost(inst);
        if child.is_feasible(inst) && cost < best_cost {
            best_cost = cost;
            best_child = Some(child);
        }
    }

    best_child
}

#[cfg(test)]
mod tests {
    use super::{eax_crossover, lcs_len, lcs_srex_crossover};
    use crate::instance::tests::make_test_instance;
    use crate::solution::Solution;
    use rand::SeedableRng;
    use rand::rngs::SmallRng;

    #[test]
    fn lcs_length_works_on_simple_sequences() {
        let a = vec![1u64, 2, 3, 4];
        let b = vec![0u64, 2, 3, 4, 5];
        assert_eq!(lcs_len(&a, &b), 3);
    }

    #[test]
    fn lcs_srex_produces_feasible_child_on_tiny_instance() {
        let m = 3usize;
        let n = 2 * m + 1;
        let mut dist = vec![0.0; n * n];
        for i in 0..n {
            for j in 0..n {
                dist[i * n + j] = if i == j { 0.0 } else { 1.0 };
            }
        }
        let demand = vec![0, 1, 1, 1, -1, -1, -1];
        let early = vec![0.0; n];
        let late = vec![10_000.0; n];
        let service = vec![0.0; n];
        let inst = make_test_instance(m, 3, &dist, &demand, &early, &late, &service);

        let parent_a = Solution {
            routes: vec![vec![0, 1, 4, 2, 5, 0], vec![0, 3, 6, 0]],
            unassigned: Vec::new(),
        };
        let parent_b = Solution {
            routes: vec![vec![0, 2, 5, 3, 6, 0], vec![0, 1, 4, 0]],
            unassigned: Vec::new(),
        };

        let mut rng = SmallRng::seed_from_u64(7);
        let child = lcs_srex_crossover(&inst, &parent_a, &parent_b, &mut rng)
            .expect("expected at least one child");
        assert!(child.is_feasible(&inst));
    }

    #[test]
    fn eax_produces_feasible_child_on_tiny_instance() {
        let m = 3usize;
        let n = 2 * m + 1;
        let mut dist = vec![0.0; n * n];
        for i in 0..n {
            for j in 0..n {
                dist[i * n + j] = if i == j { 0.0 } else { 1.0 };
            }
        }
        let demand = vec![0, 1, 1, 1, -1, -1, -1];
        let early = vec![0.0; n];
        let late = vec![10_000.0; n];
        let service = vec![0.0; n];
        let inst = make_test_instance(m, 3, &dist, &demand, &early, &late, &service);

        let parent_a = Solution {
            routes: vec![vec![0, 1, 4, 2, 5, 0], vec![0, 3, 6, 0]],
            unassigned: Vec::new(),
        };
        let parent_b = Solution {
            routes: vec![vec![0, 2, 5, 3, 6, 0], vec![0, 1, 4, 0]],
            unassigned: Vec::new(),
        };

        let mut rng = SmallRng::seed_from_u64(42);
        // EAX should produce at least one feasible child
        if let Some(child) = eax_crossover(&inst, &parent_a, &parent_b, &mut rng) {
            assert!(child.is_feasible(&inst));
        }
        // Also test with different seeds for robustness
        for seed in 0..10 {
            let mut rng2 = SmallRng::seed_from_u64(seed);
            if let Some(child) = eax_crossover(&inst, &parent_a, &parent_b, &mut rng2) {
                assert!(
                    child.is_feasible(&inst),
                    "infeasible child with seed {seed}"
                );
            }
        }
    }
}