gam_solve/

seeding.rs

use ndarray::Array1;
use std::collections::HashSet;

pub use gam_problem::{SeedConfig, SeedRiskProfile};
use crate::estimate::EstimationError;
use gam_problem::OrderedRhoBounds;

fn add_seed_dedup(seeds: &mut Vec<Array1<f64>>, seen: &mut HashSet<Vec<u64>>, seed: Array1<f64>) {
    let key: Vec<u64> = seed.iter().map(|&v| v.to_bits()).collect();
    if seen.insert(key) {
        seeds.push(seed);
    }
}

fn safe_ln_pos(x: f64) -> Option<f64> {
    if x.is_finite() && x > 0.0 {
        Some(x.ln())
    } else {
        None
    }
}

fn spde_rho_triplet_from_log_tau_log_kappa_nu(
    log_tau: f64,
    log_kappa: f64,
    nu: f64,
    bounds: OrderedRhoBounds,
) -> Option<Array1<f64>> {
    if !(nu.is_finite() && nu > 1.0) {
        return None;
    }
    let logc0 = 0.0;
    let logc1 = safe_ln_pos(nu)?;
    let logc2 = safe_ln_pos(0.5 * nu * (nu - 1.0))?;
    let rho0 = bounds.clamp(log_tau + logc0 + 2.0 * nu * log_kappa);
    let rho1 = bounds.clamp(log_tau + logc1 + 2.0 * (nu - 1.0) * log_kappa);
    let rho2 = bounds.clamp(log_tau + logc2 + 2.0 * (nu - 2.0) * log_kappa);
    Some(Array1::from_vec(vec![rho0, rho1, rho2]))
}

fn add_spde_manifold_seeds(
    seeds: &mut Vec<Array1<f64>>,
    seen: &mut HashSet<Vec<u64>>,
    bounds: OrderedRhoBounds,
    heuristic_rhos: Option<&[f64]>,
    primary: &Array1<f64>,
) {
    if primary.len() != 3 {
        return;
    }
    // Broad default manifold grid in (log_tau, log_kappa, nu).
    let tau_anchors = [primary[2], 0.0, -2.0, 2.0];
    let log_kappa_grid = [-2.0, -1.0, 0.0, 1.0, 2.0];
    let nu_grid = [1.25, 1.5, 2.0, 2.5, 3.0, 4.0];
    for &tau in &tau_anchors {
        for &lk in &log_kappa_grid {
            for &nu in &nu_grid {
                if let Some(seed) = spde_rho_triplet_from_log_tau_log_kappa_nu(tau, lk, nu, bounds)
                {
                    add_seed_dedup(seeds, seen, seed);
                }
            }
        }
    }

    // Data-informed anchor: convert the rho seed to lambdas, then invert to
    // (nu, kappa^2, tau) when feasible.
    if let Some(vals) = heuristic_rhos
        && vals.len() == 3
    {
        let l0 = vals[0].exp();
        let l1 = vals[1].exp();
        let l2 = vals[2].exp();
        if l0.is_finite() && l1.is_finite() && l2.is_finite() && l0 > 1e-12 && l2 > 1e-12 {
            let r = (l1 * l1) / (l0 * l2);
            if r > 2.0 {
                let nu = r / (r - 2.0);
                let kappa2 = l1 / ((r - 2.0) * l2);
                if nu.is_finite() && nu > 1.0 && kappa2.is_finite() && kappa2 > 0.0 {
                    let log_kappa = 0.5 * kappa2.ln();
                    let c2 = 0.5 * nu * (nu - 1.0);
                    if c2.is_finite() && c2 > 0.0 {
                        let log_tau = (l2 / (c2 * kappa2.powf(nu - 2.0))).max(1e-12).ln();
                        let local_nu = [nu, (nu - 0.3).max(1.05), nu + 0.3];
                        let local_tau = [log_tau, log_tau - 1.0, log_tau + 1.0];
                        let local_kappa = [log_kappa, log_kappa - 0.5, log_kappa + 0.5];
                        for &t in &local_tau {
                            for &lk in &local_kappa {
                                for &n in &local_nu {
                                    if let Some(seed) =
                                        spde_rho_triplet_from_log_tau_log_kappa_nu(t, lk, n, bounds)
                                    {
                                        add_seed_dedup(seeds, seen, seed);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

fn add_first_order_fallback_seeds(
    seeds: &mut Vec<Array1<f64>>,
    seen: &mut HashSet<Vec<u64>>,
    bounds: OrderedRhoBounds,
    heuristic_rhos: Option<&[f64]>,
) {
    // Degenerate λ2 -> 0 fallback (first-order mass+tension):
    // λ0 = τ κ^2, λ1 = τ, λ2 ≈ 0.
    let rho2_floor = bounds.lower();
    let default_log_kappa = [-2.0, -1.0, 0.0, 1.0];
    let default_log_tau = [0.0, -2.0, 2.0];
    for &t in &default_log_tau {
        for &lk in &default_log_kappa {
            let rho0 = bounds.clamp(t + 2.0 * lk);
            let rho1 = bounds.clamp(t);
            add_seed_dedup(seeds, seen, Array1::from_vec(vec![rho0, rho1, rho2_floor]));
        }
    }
    if let Some(vals) = heuristic_rhos
        && vals.len() == 3
        && vals[0].is_finite()
        && vals[1].is_finite()
    {
        let l0 = vals[0].exp();
        let l1 = vals[1].exp();
        let kappa2 = l0 / l1;
        if kappa2.is_finite() && kappa2 > 0.0 {
            let lk = 0.5 * kappa2.ln();
            let t = vals[1];
            let rho0 = bounds.clamp(t + 2.0 * lk);
            let rho1 = bounds.clamp(t);
            add_seed_dedup(seeds, seen, Array1::from_vec(vec![rho0, rho1, rho2_floor]));
        }
    }
}

fn add_nu2_reverse_manifold_seeds(
    seeds: &mut Vec<Array1<f64>>,
    seen: &mut HashSet<Vec<u64>>,
    bounds: OrderedRhoBounds,
    primary: &Array1<f64>,
) {
    if primary.len() != 3 {
        return;
    }
    let ln_two = 2.0_f64.ln();
    let tau_anchors = [primary[2], 0.0, -2.0, 2.0];
    let log_kappa_grid = [-2.0, -1.0, 0.0, 1.0, 2.0];
    for &tau_rho in &tau_anchors {
        for &log_kappa in &log_kappa_grid {
            // Continuous-order reverse map at nu=2:
            // lambda0 = tau * kappa^4, lambda1 = tau * 2*kappa^2, lambda2 = tau.
            let rho2 = bounds.clamp(tau_rho);
            let rho1 = bounds.clamp(tau_rho + ln_two + 2.0 * log_kappa);
            let rho0 = bounds.clamp(tau_rho + 4.0 * log_kappa);
            add_seed_dedup(seeds, seen, Array1::from_vec(vec![rho0, rho1, rho2]));
        }
    }
}

fn halton(mut index: usize, base: usize) -> f64 {
    let mut f = 1.0_f64;
    let mut r = 0.0_f64;
    while index > 0 {
        f /= base as f64;
        r += f * (index % base) as f64;
        index /= base;
    }
    r
}

fn first_primes(n: usize) -> Vec<usize> {
    let mut primes = Vec::with_capacity(n);
    let mut x = 2usize;
    while primes.len() < n {
        let mut is_prime = true;
        let mut d = 2usize;
        while d * d <= x {
            if x.is_multiple_of(d) {
                is_prime = false;
                break;
            }
            d += 1;
        }
        if is_prime {
            primes.push(x);
        }
        x += 1;
    }
    primes
}

pub fn generate_rho_candidates(
    num_penalties: usize,
    heuristic_rhos: Option<&[f64]>,
    config: &SeedConfig,
) -> Result<Vec<Array1<f64>>, EstimationError> {
    let mut seeds = Vec::new();
    let mut seen: HashSet<Vec<u64>> = HashSet::new();

    // Validate the seed ρ-box ONCE, at the boundary where it enters the
    // candidate lattice, and REFUSE an inverted/non-finite interval rather than
    // silently reordering it (#2379). Every clamp below then operates on an
    // interval that is ordered by construction. This is the same disease and the
    // same cure as the outer-optimizer prepass: an inverted box means two
    // independently-owned constants have drifted apart, and a swap would make the
    // lattice explore a different, silently substituted box.
    let bounds = OrderedRhoBounds::new(config.bounds.0, config.bounds.1)?;
    let max_seeds = config.max_seeds.max(1);
    let risk_shift = config.risk_profile.anchor_rho_shift();

    if num_penalties == 0 {
        add_seed_dedup(&mut seeds, &mut seen, Array1::<f64>::zeros(0));
        return Ok(seeds);
    }

    // Prefer a full heuristic vector (length == k) as the primary anchor.
    // Values are already in the outer optimizer's rho/theta parameter space.
    let num_aux = config.num_auxiliary_trailing.min(num_penalties);
    let num_smoothing = num_penalties - num_aux;
    let aux_initial: Vec<f64> = if num_aux > 0 {
        heuristic_rhos
            .filter(|h| h.len() == num_penalties)
            .map(|h| {
                h[num_smoothing..]
                    .iter()
                    .copied()
                    .map(|v| bounds.clamp(v))
                    .collect()
            })
            .unwrap_or_else(|| vec![0.0; num_aux])
    } else {
        Vec::new()
    };
    let heuristic_rhovec: Option<Array1<f64>> = heuristic_rhos.and_then(|vals| {
        if vals.len() == num_penalties {
            Some(Array1::from_iter(
                vals[..num_smoothing]
                    .iter()
                    .copied()
                    .map(|v| bounds.clamp(v))
                    .chain(
                        vals[num_smoothing..]
                            .iter()
                            .copied()
                            .map(|v| bounds.clamp(v)),
                    ),
            ))
        } else {
            None
        }
    });

    let primary = heuristic_rhovec.clone().unwrap_or_else(|| {
        Array1::<f64>::from_elem(num_penalties, bounds.clamp(risk_shift))
    });
    add_seed_dedup(&mut seeds, &mut seen, primary.clone());
    // Always include neutral baseline independently of heuristic anchor.
    add_seed_dedup(&mut seeds, &mut seen, Array1::zeros(num_penalties));
    // Generalized and survival models can hit PIRLS separation at moderate
    // smoothing levels. Put an aggressively over-smoothed isotropic seed near
    // the front so startup validation can still find a stable basin.
    match config.risk_profile {
        SeedRiskProfile::Gaussian | SeedRiskProfile::GaussianLocationScale => {}
        SeedRiskProfile::GeneralizedLinear | SeedRiskProfile::Survival => {
            add_seed_dedup(
                &mut seeds,
                &mut seen,
                Array1::from_elem(num_penalties, bounds.upper()),
            );
        }
    }
    // For exactly three smoothing penalties (mass/tension/stiffness), inject
    // physically coherent manifold seeds in rho-space:
    // - general SPDE manifold over (log_tau, log_kappa, nu),
    // - nu=2 reverse-map seeds,
    // - first-order fallback seeds (lambda2 near lower bound).
    if num_smoothing == 3 {
        let smoothing_primary =
            Array1::from_vec(primary.iter().take(num_smoothing).copied().collect());
        let smoothing_heuristic_lambdas = heuristic_rhos.and_then(|vals| {
            if vals.len() >= num_smoothing {
                Some(&vals[..num_smoothing])
            } else {
                None
            }
        });
        let mut spde_prefix_seeds = Vec::new();
        let mut spde_prefix_seen: HashSet<Vec<u64>> = HashSet::new();
        // Guarantee a first-order fallback anchor regardless of later truncation.
        add_seed_dedup(
            &mut spde_prefix_seeds,
            &mut spde_prefix_seen,
            Array1::from_vec(vec![primary[0], primary[1], bounds.lower()]),
        );
        // Ensure a nu=2-consistent seed is always present before broader grids.
        add_nu2_reverse_manifold_seeds(
            &mut spde_prefix_seeds,
            &mut spde_prefix_seen,
            bounds,
            &smoothing_primary,
        );
        add_first_order_fallback_seeds(
            &mut spde_prefix_seeds,
            &mut spde_prefix_seen,
            bounds,
            smoothing_heuristic_lambdas,
        );
        add_spde_manifold_seeds(
            &mut spde_prefix_seeds,
            &mut spde_prefix_seen,
            bounds,
            smoothing_heuristic_lambdas,
            &smoothing_primary,
        );
        for prefix_seed in spde_prefix_seeds {
            let mut seed = Array1::<f64>::zeros(num_penalties);
            for i in 0..num_smoothing {
                seed[i] = prefix_seed[i];
            }
            for (i, &v) in aux_initial.iter().enumerate() {
                seed[num_smoothing + i] = v;
            }
            add_seed_dedup(&mut seeds, &mut seen, seed);
        }
    }

    // Broad symmetric baselines around the center to guarantee global coverage.
    for &center in config.risk_profile.baseline_centers() {
        add_seed_dedup(
            &mut seeds,
            &mut seen,
            Array1::from_elem(num_penalties, bounds.clamp(center)),
        );
    }

    let dims_to_touch = num_penalties.min(12);
    let step_base = if num_penalties <= 4 {
        2.0
    } else if num_penalties <= 12 {
        2.5
    } else {
        3.0
    };
    let high_dim_cluster_threshold = 10usize;

    if num_penalties >= high_dim_cluster_threshold {
        // High-dimensional path: probe relative scaling conflicts by clustering
        // penalties into low/high heuristic-magnitude groups.
        let mut sorted_idx: Vec<usize> = (0..num_penalties).collect();
        sorted_idx.sort_by(|&i, &j| primary[i].total_cmp(&primary[j]));

        let cluster_size = (num_penalties / 3).max(1);
        let small_end = cluster_size.min(num_penalties);
        let large_start = num_penalties.saturating_sub(cluster_size);
        let small_cluster = &sorted_idx[..small_end];
        let large_cluster = &sorted_idx[large_start..];

        let small_scale = step_base;
        let large_scale = step_base + 0.75;

        let mut conflict_a = primary.clone();
        for &i in large_cluster {
            conflict_a[i] = bounds.clamp(primary[i] + large_scale);
        }
        for &i in small_cluster {
            conflict_a[i] = bounds.clamp(primary[i] - small_scale);
        }
        add_seed_dedup(&mut seeds, &mut seen, conflict_a);

        let mut conflict_b = primary.clone();
        for &i in large_cluster {
            conflict_b[i] = bounds.clamp(primary[i] - large_scale);
        }
        for &i in small_cluster {
            conflict_b[i] = bounds.clamp(primary[i] + small_scale);
        }
        add_seed_dedup(&mut seeds, &mut seen, conflict_b);

        let mut heavy_up = primary.clone();
        for &i in large_cluster {
            heavy_up[i] = bounds.clamp(primary[i] + large_scale);
        }
        add_seed_dedup(&mut seeds, &mut seen, heavy_up);

        let mut light_down = primary.clone();
        for &i in small_cluster {
            light_down[i] = bounds.clamp(primary[i] - small_scale);
        }
        add_seed_dedup(&mut seeds, &mut seen, light_down);
    } else {
        // Low-dimensional path: coordinate and sparse pair probes are still cheap.
        for i in 0..dims_to_touch {
            let scale = step_base + 0.25 * primary[i].abs().min(8.0);
            for dir in [-1.0, 1.0] {
                let mut s = primary.clone();
                s[i] = bounds.clamp(primary[i] + dir * scale);
                add_seed_dedup(&mut seeds, &mut seen, s);
            }
        }

        let pair_dims = num_penalties.min(6);
        for i in 0..pair_dims {
            for j in (i + 1)..pair_dims {
                let mut s1 = primary.clone();
                s1[i] = bounds.clamp(primary[i] + step_base);
                s1[j] = bounds.clamp(primary[j] - step_base);
                add_seed_dedup(&mut seeds, &mut seen, s1);

                let mut s2 = primary.clone();
                s2[i] = bounds.clamp(primary[i] - step_base);
                s2[j] = bounds.clamp(primary[j] + step_base);
                add_seed_dedup(&mut seeds, &mut seen, s2);
            }
        }
    }

    // Global shrink/expand sweeps from the anchor to probe over/under-smoothing regimes.
    // The flexible (negative-shift) side MUST be probed as densely as the
    // over-smoothing side: the seed-screening proxy is a capped-inner-iteration
    // fit, and an over-smoothed seed converges trivially under that cap (its
    // coefficients collapse into the penalty null space, the LAML is locally
    // flat), so screening systematically ranks over-smoothed seeds first
    // (documented in `rank_seeds_with_screening`). For a GeneralizedLinear /
    // Survival model whose true optimum is flexible (e.g. a smooth Poisson
    // tensor surface that genuinely needs ~10 effective df), a seed grid that
    // only sweeps the over-smoothing side leaves the flexible basin unprobed,
    // so none of the few full-budget solves ever lands in it and the fit
    // over-smooths (#1082/#1373). Symmetric negative shifts give the flexible
    // basin a candidate; the keep-best multi-start then retains it only if it
    // actually scores better, so this can never worsen a fit — it only lets the
    // optimizer SEE the lower-λ basin. Over-smoothed seeds remain present (and
    // earlier in the list) so PIRLS-separation startup stability is unchanged.
    for &shift in config.risk_profile.global_shifts() {
        let swept = primary.mapv(|v| bounds.clamp(v + shift));
        add_seed_dedup(&mut seeds, &mut seen, swept);
    }

    // #1464 over-smoothing probe: an ABSOLUTE high-λ start on every smoothing
    // dimension (auxiliary dims left at the anchor's values; they are re-pinned
    // below). The global shift sweeps above reach only ≈ +4 from the anchor, so
    // a collapsing-kernel smooth whose true REML optimum is a large λ would never
    // be seeded into its over-smoothing basin. This puts a candidate IN it; the
    // keep-best multistart adopts it only when it scores strictly better, so it
    // can never worsen a fit. `None` (the default) skips this entirely.
    if let Some(probe_rho) = config.over_smoothing_probe_rho {
        let mut probe = primary.clone();
        for j in 0..num_smoothing {
            probe[j] = bounds.clamp(probe_rho);
        }
        add_seed_dedup(&mut seeds, &mut seen, probe);
    }

    // Low-discrepancy exploratory seeds around the anchor for basin discovery.
    // These are still deterministic and do not encode any solver-side bias.
    let exploratory = max_seeds.saturating_sub(seeds.len()).min(8);
    if exploratory > 0 {
        let primes = first_primes(num_penalties.max(1));
        let amp = config.risk_profile.exploratory_amplitude();
        for t in 0..exploratory {
            let mut s = primary.clone();
            for i in 0..num_penalties {
                let u = halton(t + 1, primes[i]); // (0,1)
                let centered = 2.0 * u - 1.0; // (-1,1)
                s[i] = bounds.clamp(primary[i] + amp * centered);
            }
            add_seed_dedup(&mut seeds, &mut seen, s);
        }
    }

    // Pin auxiliary trailing dimensions to their initial values in every seed.
    // Auxiliary params (e.g. SAS epsilon, log_delta) live in a different
    // parameter space than log-smoothing rho and must not be swept by the
    // smoothing seeding grid.  After pinning we re-dedup because seeds that
    // differed only in the (now-overwritten) auxiliary dimensions collapse.
    if num_aux > 0 {
        for seed in &mut seeds {
            for (i, &v) in aux_initial.iter().enumerate() {
                seed[num_smoothing + i] = v;
            }
        }
        let mut deduped = Vec::new();
        let mut seen2: HashSet<Vec<u64>> = HashSet::new();
        for seed in seeds {
            let key: Vec<u64> = seed.iter().map(|&v| v.to_bits()).collect();
            if seen2.insert(key) {
                deduped.push(seed);
            }
        }
        seeds = deduped;
    }

    if seeds.len() > max_seeds {
        seeds.truncate(max_seeds);
    }

    if seeds.is_empty() {
        seeds.push(Array1::<f64>::zeros(num_penalties));
    }

    Ok(seeds)
}

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

    /// #2379 (mechanism B): the candidate lattice REFUSES an inverted seed box
    /// with a typed error rather than silently reordering it — the same contract
    /// the outer-optimizer prepass now enforces. A future constant drift that
    /// inverts `SeedConfig.bounds` fails loud here instead of quietly generating
    /// a lattice over a different, substituted box.
    #[test]
    fn generate_rho_candidates_refuses_inverted_bounds() {
        let cfg = SeedConfig {
            bounds: (10.0, -10.0),
            ..SeedConfig::default()
        };
        let err = generate_rho_candidates(3, None, &cfg)
            .expect_err("an inverted seed box must be refused, not swapped");
        assert!(
            matches!(err, EstimationError::InvalidInput(_)),
            "inverted seed box must be a typed InvalidInput, got {err:?}"
        );
    }

    /// The ordered default box still produces a full lattice — the refusal never
    /// fires on any real config (every `SeedConfig.bounds` in the tree is the
    /// ordered constant `(-12.0, 12.0)`).
    #[test]
    fn generate_rho_candidates_accepts_ordered_bounds() {
        let cfg = SeedConfig::default();
        let seeds = generate_rho_candidates(3, None, &cfg).expect("ordered box is accepted");
        assert!(!seeds.is_empty());
    }

    #[test]
    fn uses_full_heuristicvector_as_primary_anchor() {
        let cfg = SeedConfig {
            risk_profile: SeedRiskProfile::Gaussian,
            ..SeedConfig::default()
        };
        let heur = [-2.0, 0.0, 2.0];
        let seeds = generate_rho_candidates(3, Some(&heur), &cfg).expect("ordered seed bounds");
        assert!(!seeds.is_empty());
        let first = &seeds[0];
        assert_eq!(first.len(), 3);
        assert!((first[0] - heur[0]).abs() < 1e-12);
        assert!((first[1] - heur[1]).abs() < 1e-12);
        assert!((first[2] - heur[2]).abs() < 1e-12);
    }

    #[test]
    fn high_dim_uses_cluster_conflict_probeswithout_exploding() {
        let cfg = SeedConfig {
            max_seeds: 18,
            risk_profile: SeedRiskProfile::GeneralizedLinear,
            ..SeedConfig::default()
        };
        let heur = [-6.0, -5.0, -4.0, 0.0, 2.0, 4.0, -3.0, 0.0, 3.0, 5.0];
        let seeds = generate_rho_candidates(10, Some(&heur), &cfg).expect("ordered seed bounds");
        assert!(seeds.len() <= 18);
        // Presence of at least one asymmetric cluster-conflict seed:
        // some coordinates increased while others decreased vs primary.
        let primary = &seeds[0];
        let has_conflict = seeds.iter().skip(1).any(|s| {
            let mut any_up = false;
            let mut any_down = false;
            for i in 0..s.len() {
                if s[i] > primary[i] {
                    any_up = true;
                } else if s[i] < primary[i] {
                    any_down = true;
                }
            }
            any_up && any_down
        });
        assert!(has_conflict);
    }

    #[test]
    fn includes_neutralzero_seed() {
        let cfg = SeedConfig::default();
        let seeds = generate_rho_candidates(5, None, &cfg).expect("ordered seed bounds");
        let haszero = seeds
            .iter()
            .any(|s| s.iter().all(|v| (*v - 0.0).abs() < 1e-12));
        assert!(haszero);
    }

    #[test]
    fn generalized_linear_seeds_include_early_stability_retreat_seed() {
        let cfg = SeedConfig {
            risk_profile: SeedRiskProfile::GeneralizedLinear,
            ..SeedConfig::default()
        };
        let seeds = generate_rho_candidates(3, None, &cfg).expect("ordered seed bounds");
        let retreat = Array1::from_elem(3, cfg.bounds.1);
        let retreat_idx = seeds
            .iter()
            .position(|seed| seed == retreat)
            .expect("generalized-linear seeds should include an upper-bound retreat seed");
        assert!(
            retreat_idx <= 2,
            "retreat seed should be available before broader exploratory seeds: {retreat_idx}"
        );
    }

    #[test]
    fn three_penalty_seeds_include_nu2_reverse_manifold_triplets() {
        let cfg = SeedConfig::default();
        let seeds = generate_rho_candidates(3, None, &cfg).expect("ordered seed bounds");
        let ln4 = 4.0_f64.ln();
        let has_nu2_manifold_seed = seeds
            .iter()
            .any(|s| s.len() == 3 && ((2.0 * s[1] - s[0] - s[2]) - ln4).abs() < 1e-8);
        assert!(has_nu2_manifold_seed);
    }

    #[test]
    fn three_penalty_seeds_include_general_spde_manifold_points() {
        let cfg = SeedConfig::default();
        let heur = [2.0, 10.0, 3.0];
        let seeds = generate_rho_candidates(3, Some(&heur), &cfg).expect("ordered seed bounds");
        let has_non_nu2 = seeds.iter().any(|s| {
            // For nu=2, 2*rho1-rho0-rho2 = ln(4).
            // General nu manifold should include points away from ln(4).
            s.len() == 3 && ((2.0 * s[1] - s[0] - s[2]) - 4.0_f64.ln()).abs() > 1e-3
        });
        assert!(has_non_nu2);
    }

    #[test]
    fn three_penalty_seeds_include_first_order_fallbackwith_rho2_floor() {
        let cfg = SeedConfig {
            bounds: (-12.0, 12.0),
            ..SeedConfig::default()
        };
        let seeds = generate_rho_candidates(3, None, &cfg).expect("ordered seed bounds");
        let has_floor = seeds
            .iter()
            .any(|s| s.len() == 3 && (s[2] - (-12.0)).abs() < 1e-12);
        assert!(has_floor);
    }

    #[test]
    fn auxiliary_trailing_dims_pinned_to_initial_values() {
        // Simulate SAS optimization: 2 smoothing dims + 2 auxiliary dims
        // (epsilon=0, log_delta=0).  The heuristic vector is in rho/theta
        // space for both smoothing and auxiliary dimensions.
        let cfg = SeedConfig {
            num_auxiliary_trailing: 2,
            risk_profile: SeedRiskProfile::GeneralizedLinear,
            ..SeedConfig::default()
        };
        let heur = [0.0, 10.0_f64.ln(), 0.0, 0.0]; // rhos + SAS initials
        let seeds = generate_rho_candidates(4, Some(&heur), &cfg).expect("ordered seed bounds");
        assert!(!seeds.is_empty());
        // EVERY seed must have the auxiliary dims pinned to 0.0.
        for (idx, seed) in seeds.iter().enumerate() {
            assert_eq!(seed.len(), 4);
            assert!(
                (seed[2] - 0.0).abs() < 1e-12 && (seed[3] - 0.0).abs() < 1e-12,
                "seed {} has auxiliary dims [{}, {}], expected [0, 0]",
                idx,
                seed[2],
                seed[3],
            );
        }
        // The smoothing dims should NOT all be zero (some seeds should vary them).
        let has_nonzero_smoothing = seeds
            .iter()
            .any(|s| s[0].abs() > 1e-12 || s[1].abs() > 1e-12);
        assert!(has_nonzero_smoothing);
    }

    #[test]
    fn auxiliary_dims_dedup_collapses_identical_seeds() {
        // With auxiliary pinning, seeds that differed only in aux dims
        // should collapse to a single seed.
        let cfg = SeedConfig {
            num_auxiliary_trailing: 1,
            max_seeds: 32,
            risk_profile: SeedRiskProfile::GeneralizedLinear,
            ..SeedConfig::default()
        };
        let seeds_with_aux = generate_rho_candidates(3, None, &cfg).expect("ordered seed bounds");
        let cfg_no_aux = SeedConfig {
            num_auxiliary_trailing: 0,
            max_seeds: 32,
            risk_profile: SeedRiskProfile::GeneralizedLinear,
            ..SeedConfig::default()
        };
        let seeds_without_aux = generate_rho_candidates(3, None, &cfg_no_aux).expect("ordered seed bounds");
        // Aux pinning causes many seeds to collapse, so fewer unique seeds.
        assert!(seeds_with_aux.len() <= seeds_without_aux.len());
    }
}