sphereql_embed/

category.rs

use std::collections::HashMap;

use sphereql_core::{SphericalPoint, angular_distance, cosine_similarity};

use crate::config::PipelineConfig;
use crate::kernel_pca::KernelPcaProjection;
use crate::projection::{PcaProjection, Projection};
use crate::spatial_quality::SpatialQuality;
use crate::types::{Embedding, RadialStrategy};

// ── Category summary ───────────────────────────────────────────────────

/// Aggregate statistics for a single category on the outer sphere.
///
/// Computed from the projected positions of all items in that category.
/// Every category gets a summary regardless of size — this is the
/// foundation of the Category Enrichment Layer.
#[derive(Debug, Clone)]
pub struct CategorySummary {
    /// Category name (as provided by the user).
    pub name: String,
    /// Indices of member items in the pipeline's item list.
    pub member_indices: Vec<usize>,
    /// Mean embedding in high-dimensional space (pre-projection).
    /// Length = embedding dimensionality.
    pub centroid_embedding: Vec<f64>,
    /// The centroid projected onto the outer sphere.
    pub centroid_position: SphericalPoint,
    /// Mean angular distance (radians) of members from the centroid
    /// on the projected sphere. Measures how "spread out" the category is.
    pub angular_spread: f64,
    /// 1.0 / (1.0 + angular_spread). Higher = tighter cluster.
    /// Normalized to (0, 1].
    pub cohesion: f64,
    /// Number of member items.
    pub member_count: usize,
    /// Solid angle of this category's cap on S² (steradians).
    pub cap_area: f64,
    /// Fraction of this category's cap not overlapped by any other. [0, 1].
    pub exclusivity: f64,
    /// Voronoi cell area on S² (steradians).
    pub voronoi_area: f64,
    /// Items per steradian of Voronoi cell.
    pub territorial_efficiency: f64,
    /// Mean territorial-adjusted bridge strength across all outgoing edges.
    /// 0.0 if this category has no neighbors. Populated after the graph
    /// is built in [`CategoryLayer::build`].
    pub bridge_quality: f64,
}

// ── Bridge items ───────────────────────────────────────────────────────

/// Semantic relation type for a bridge edge.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, serde::Serialize, serde::Deserialize)]
pub enum RelationType {
    Studies,
    AppliesTo,
    Enables,
    CausedBy,
    Contains,
    SharedMethod,
    HistoricalLink,
}

impl RelationType {
    pub fn name(self) -> &'static str {
        match self {
            Self::Studies => "studies",
            Self::AppliesTo => "applies_to",
            Self::Enables => "enables",
            Self::CausedBy => "caused_by",
            Self::Contains => "contains",
            Self::SharedMethod => "shared_method",
            Self::HistoricalLink => "historical_link",
        }
    }
}

/// Infer a [`RelationType`] for a bridge from its concept label.
///
/// Source and target category names are accepted for future use but
/// currently unused — keyword cues in the label drive the inference.
fn infer_relation(
    concept_label: &str,
    _source_category: &str,
    _target_category: &str,
) -> Option<RelationType> {
    let label = concept_label.to_lowercase();

    if label.contains("history of")
        || label.contains("historical")
        || label.contains("origins of")
        || label.contains("evolution of")
    {
        return Some(RelationType::HistoricalLink);
    }

    if label.contains("method")
        || label.contains("technique")
        || label.contains("analysis")
        || label.contains("spectroscop")
        || label.contains("microscop")
        || label.contains("imaging")
        || label.contains("modeling")
        || label.contains("simulation")
        || label.contains("measurement")
    {
        return Some(RelationType::SharedMethod);
    }

    if label.contains("applications")
        || label.contains("applied")
        || label.contains(" in ")
        || label.contains(" for ")
    {
        return Some(RelationType::AppliesTo);
    }

    if label.contains("studies")
        || label.contains("research")
        || label.contains("science of")
        || label.contains("theory of")
    {
        return Some(RelationType::Studies);
    }

    if label.starts_with("bio")
        || label.starts_with("geo")
        || label.starts_with("neuro")
        || label.starts_with("astro")
        || label.starts_with("electro")
        || label.starts_with("psycho")
        || label.starts_with("socio")
        || label.starts_with("eco")
    {
        return Some(RelationType::Contains);
    }

    if label.contains("foundation")
        || label.contains("prerequisite")
        || label.contains("basis")
        || label.contains("fundamentals")
        || label.contains("principles")
    {
        return Some(RelationType::Enables);
    }

    None
}

/// Quality classification for a bridge item.
///
/// Assigned after all bridges are collected, comparing each bridge's
/// strength against the corpus-wide median and the pair's territorial
/// separation on S².
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum BridgeClassification {
    /// Strength at or above the median and the source/target caps are
    /// spatially distinct — a real cross-domain connector.
    Genuine,
    /// The two categories overlap heavily on S² (low territorial factor).
    /// This bridge is more shared-territory noise than a genuine connector.
    OverlapArtifact,
    /// Strength below median in a territorially clean pair — a real
    /// connection but not a strong one.
    Weak,
}

/// An item that semantically spans two categories.
///
/// Bridge items are closer to a foreign category's centroid than to the
/// median distance within their own category. They are the conceptual
/// connectors that make cross-domain paths meaningful.
#[derive(Debug, Clone)]
pub struct BridgeItem {
    /// Index of this item in the pipeline's item list.
    pub item_index: usize,
    /// The item's own category index.
    pub source_category: usize,
    /// The foreign category this item bridges toward.
    pub target_category: usize,
    /// Cosine similarity to own category centroid (in high-D space).
    pub affinity_to_source: f64,
    /// Cosine similarity to foreign category centroid (in high-D space).
    pub affinity_to_target: f64,
    /// Bridge strength: harmonic mean of the two affinities.
    /// Higher = equally strong connection to both domains.
    pub bridge_strength: f64,
    /// Quality label assigned after the full bridge set is observed.
    pub classification: BridgeClassification,
    /// Semantic relation type, inferred heuristically. `None` when untyped.
    pub relation: Option<RelationType>,
}

// ── Category graph ─────────────────────────────────────────────────────

/// Edge in the category adjacency graph.
#[derive(Debug, Clone)]
pub struct CategoryEdge {
    /// Index of the neighbor category.
    pub target: usize,
    /// Angular distance between category centroids on the sphere.
    pub centroid_distance: f64,
    /// Number of bridge items connecting these two categories.
    pub bridge_count: usize,
    /// Strongest bridge connecting these two categories (0.0 if no bridges).
    pub max_bridge_strength: f64,
    /// Mean bridge strength across all bridges (0.0 if no bridges).
    pub mean_bridge_strength: f64,
    /// Combined edge weight (lower = more connected).
    /// Computed as `centroid_distance / (1 + bridge_count * mean_bridge_strength)`.
    /// This prefers fewer strong bridges over many weak ones.
    pub weight: f64,
}

/// The full category adjacency graph.
#[derive(Debug, Clone)]
pub struct CategoryGraph {
    /// Adjacency list: `adjacency[i]` contains edges from category i.
    pub adjacency: Vec<Vec<CategoryEdge>>,
    /// Bridge items keyed by (source_category, target_category).
    /// Sorted by descending bridge_strength within each pair.
    pub bridges: HashMap<(usize, usize), Vec<BridgeItem>>,
}

// ── Category-level concept path ────────────────────────────────────────

/// A step in a category-level concept path.
#[derive(Debug, Clone)]
pub struct CategoryPathStep {
    /// Category index.
    pub category_index: usize,
    /// Category name.
    pub category_name: String,
    /// Cumulative distance from the start.
    pub cumulative_distance: f64,
    /// Bridge items connecting this step to the next (empty for the last step).
    pub bridges_to_next: Vec<BridgeItem>,
    /// Spatial confidence of this hop: bridge_strength × territorial_factor.
    /// Higher = more trustworthy transition. 0.0 for the last step.
    pub hop_confidence: f64,
}

/// Result of a category-level concept path query.
#[derive(Debug, Clone)]
pub struct CategoryPath {
    /// Ordered steps from source to target category.
    pub steps: Vec<CategoryPathStep>,
    /// Total path distance.
    pub total_distance: f64,
    /// Product of all hop confidences along the path.
    /// Low values indicate the path routes through shaky connections.
    pub path_confidence: f64,
}

// ── Inner sphere (Phase 2) ─────────────────────────────────────────────

/// The projection type used for a category's inner sphere.
///
/// Wraps either a linear PCA or kernel PCA projection, chosen
/// automatically based on the category's size and measured EVR
/// improvement over the global projection.
#[derive(Clone)]
pub enum InnerProjection {
    /// Standard linear PCA — chosen for categories meeting
    /// [`InnerSphereConfig::min_size`](crate::config::InnerSphereConfig::min_size)
    /// but below
    /// [`kernel_pca_min_size`](crate::config::InnerSphereConfig::kernel_pca_min_size),
    /// or when kernel PCA fails to improve over linear by
    /// [`min_kernel_improvement`](crate::config::InnerSphereConfig::min_kernel_improvement).
    LinearPca(PcaProjection),
    /// Gaussian kernel PCA — chosen for categories meeting
    /// [`kernel_pca_min_size`](crate::config::InnerSphereConfig::kernel_pca_min_size)
    /// where it measurably outperforms linear PCA.
    KernelPca(KernelPcaProjection),
}

impl std::fmt::Debug for InnerProjection {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::LinearPca(_) => write!(f, "LinearPca"),
            Self::KernelPca(_) => write!(f, "KernelPca"),
        }
    }
}

impl Projection for InnerProjection {
    fn project(&self, embedding: &Embedding) -> SphericalPoint {
        match self {
            Self::LinearPca(p) => p.project(embedding),
            Self::KernelPca(p) => p.project(embedding),
        }
    }
    fn project_rich(&self, embedding: &Embedding) -> crate::types::ProjectedPoint {
        match self {
            Self::LinearPca(p) => p.project_rich(embedding),
            Self::KernelPca(p) => p.project_rich(embedding),
        }
    }
    fn dimensionality(&self) -> usize {
        match self {
            Self::LinearPca(p) => p.dimensionality(),
            Self::KernelPca(p) => p.dimensionality(),
        }
    }
}

/// A category-specific inner sphere with its own optimized projection.
///
/// Only created for categories that meet all of:
/// 1. At least [`InnerSphereConfig::min_size`](crate::config::InnerSphereConfig::min_size) members
/// 2. Inner EVR improves over global subset EVR by ≥
///    [`InnerSphereConfig::min_evr_improvement`](crate::config::InnerSphereConfig::min_evr_improvement)
///
/// The inner sphere gives higher-resolution angular discrimination
/// within the category than the global outer projection can provide.
#[derive(Debug, Clone)]
pub struct InnerSphere {
    /// The category-specific projection (linear PCA or kernel PCA).
    pub projection: InnerProjection,
    /// Positions of member items in the inner sphere's coordinate system.
    /// `inner_positions[i]` corresponds to `member_indices[i]`.
    pub inner_positions: Vec<SphericalPoint>,
    /// Global item indices of the members (same order as inner_positions).
    pub member_indices: Vec<usize>,
    /// Explained variance ratio of the inner projection.
    pub explained_variance_ratio: f64,
    /// Mean certainty of these items under the global (outer) projection.
    /// Baseline for measuring improvement.
    pub global_subset_evr: f64,
    /// `explained_variance_ratio - global_subset_evr`.
    pub evr_improvement: f64,
}

/// A single item from a [`CategoryLayer::drill_down_with_projection`] query.
#[derive(Debug, Clone)]
pub struct DrillDownResult {
    /// Index of the item in the pipeline's global item list.
    pub item_index: usize,
    /// Angular distance to the query in the relevant coordinate system.
    pub distance: f64,
    /// Whether the inner sphere's projection was used (true) or the
    /// outer sphere was used as fallback (false).
    pub used_inner_sphere: bool,
}

/// Stats for a single inner sphere, returned by
/// [`inner_sphere_stats`](CategoryLayer::inner_sphere_stats).
#[derive(Debug, Clone)]
pub struct InnerSphereReport {
    /// Category name.
    pub category_name: String,
    /// Category index.
    pub category_index: usize,
    /// Number of members in the inner sphere.
    pub member_count: usize,
    /// `"LinearPca"` or `"KernelPca"`.
    pub projection_type: &'static str,
    /// Explained variance ratio of the inner projection.
    pub inner_evr: f64,
    /// Mean certainty of members under the global projection.
    pub global_subset_evr: f64,
    /// EVR improvement over global.
    pub evr_improvement: f64,
}

// ── The enrichment layer ───────────────────────────────────────────────

/// Category Enrichment Layer: aggregate statistics, inter-category graph,
/// bridge item detection, and automatic inner spheres over a projected
/// SphereQL corpus.
///
/// This is a read-only structure computed from an existing pipeline's
/// data. It adds category-level reasoning without modifying the
/// underlying projection or spatial index.
#[derive(Debug, Clone)]
pub struct CategoryLayer {
    /// One summary per unique category, in insertion order.
    pub summaries: Vec<CategorySummary>,
    /// Map from category name to index in `summaries`.
    pub name_to_index: HashMap<String, usize>,
    /// The inter-category adjacency graph.
    pub graph: CategoryGraph,
    /// Outer-sphere positions for all items (same indexing as embeddings).
    outer_positions: Vec<SphericalPoint>,
    /// Inner spheres keyed by category index. Only present for categories
    /// that meet the size and EVR-improvement thresholds.
    pub inner_spheres: HashMap<usize, InnerSphere>,
    /// Domain groups discovered from category geometry. Default
    /// `nearest()` routing chooses the closest group via these centroids.
    pub domain_groups: Vec<crate::domain_groups::DomainGroup>,
    /// Inner spheres keyed by domain-group index. Built only when the
    /// group's union of members beats the EVR threshold versus the
    /// outer projection — the realistic-grouping unit per the v2 routing
    /// design (replaces the previous one-per-category default).
    pub group_inner_spheres: HashMap<usize, InnerSphere>,
    /// Pre-computed spatial properties of the category layout on S².
    /// Used for bridge quality scoring, confidence signals, and routing.
    pub spatial_quality: SpatialQuality,
}

impl CategoryLayer {
    /// Build the category enrichment layer from pipeline data.
    ///
    /// - `categories[i]` is the category name for item i.
    /// - `embeddings[i]` is the raw embedding for item i.
    /// - `projected_positions[i]` is the spherical position on the outer sphere.
    /// - `projection` is used to project category centroids and measure
    ///   per-point certainty for inner sphere threshold decisions.
    ///
    /// Inner spheres are automatically constructed for categories that:
    /// 1. Have ≥ 20 members
    /// 2. Show ≥ 0.10 EVR improvement over the global projection
    ///
    /// Categories with ≥ 80 members additionally try kernel PCA and
    /// select it if it improves EVR by ≥ 0.05 over linear PCA.
    ///
    /// O(N·C + C²) for the base layer, plus O(n_c²·d) per inner sphere.
    ///
    /// Uses [`PipelineConfig::default`] for all tunables. Call
    /// [`Self::build_with_config`] to override.
    pub fn build<P: Projection>(
        categories: &[String],
        embeddings: &[Embedding],
        projected_positions: &[SphericalPoint],
        projection: &P,
        evr: f64,
    ) -> Self {
        Self::build_with_config(
            categories,
            embeddings,
            projected_positions,
            projection,
            evr,
            &PipelineConfig::default(),
        )
    }

    /// Configurable variant of [`Self::build`]. Threads a [`PipelineConfig`]
    /// through spatial quality, bridge classification, and inner-sphere
    /// gating.
    pub fn build_with_config<P: Projection>(
        categories: &[String],
        embeddings: &[Embedding],
        projected_positions: &[SphericalPoint],
        projection: &P,
        evr: f64,
        config: &PipelineConfig,
    ) -> Self {
        let n = categories.len();
        // Invariant: the pipeline builder validates length parity before calling
        // build_with_config. A mismatch here is a programmer error, not a
        // recoverable input error.
        assert_eq!(n, embeddings.len());
        assert_eq!(n, projected_positions.len());

        // 1. Discover unique categories in first-seen order, then drop any
        //    whose member count falls below `min_category_size`. Items in
        //    excluded categories remain projected and indexed on the sphere
        //    — they just don't get a CategorySummary, don't produce bridges,
        //    and don't participate in domain groups or spatial quality.
        let min_size = config.min_category_size;

        let mut name_to_index: HashMap<String, usize> = HashMap::new();
        let mut cat_names: Vec<String> = Vec::new();
        let mut cat_members: Vec<Vec<usize>> = Vec::new();

        // Single pass: insert in first-seen order, but defer the size filter.
        let mut raw_index: HashMap<&str, usize> = HashMap::new();
        let mut raw_names: Vec<&str> = Vec::new();
        let mut raw_members: Vec<Vec<usize>> = Vec::new();
        for (i, cat) in categories.iter().enumerate() {
            let idx = if let Some(&idx) = raw_index.get(cat.as_str()) {
                idx
            } else {
                let idx = raw_names.len();
                raw_index.insert(cat.as_str(), idx);
                raw_names.push(cat.as_str());
                raw_members.push(Vec::new());
                idx
            };
            raw_members[idx].push(i);
        }

        for (raw_name, members) in raw_names.iter().zip(raw_members) {
            if members.len() < min_size {
                continue;
            }
            let idx = cat_names.len();
            name_to_index.insert((*raw_name).to_string(), idx);
            cat_names.push((*raw_name).to_string());
            cat_members.push(members);
        }

        let num_cats = cat_names.len();
        let dim = if n > 0 { embeddings[0].dimension() } else { 0 };

        // 2. Compute category summaries
        let mut summaries: Vec<CategorySummary> = Vec::with_capacity(num_cats);

        for (ci, name) in cat_names.iter().enumerate() {
            let members = &cat_members[ci];
            let count = members.len();

            // Centroid in high-D space
            let mut centroid_emb = vec![0.0; dim];
            for &mi in members {
                for (j, &v) in embeddings[mi].values.iter().enumerate() {
                    centroid_emb[j] += v;
                }
            }
            if count > 0 {
                for v in &mut centroid_emb {
                    *v /= count as f64;
                }
            }

            // Project the centroid
            let centroid_embedding_obj = Embedding::new(centroid_emb.clone());
            let centroid_position = projection.project(&centroid_embedding_obj);

            // Angular spread: mean angular distance of members from centroid
            let angular_spread = if count > 1 {
                let total: f64 = members
                    .iter()
                    .map(|&mi| angular_distance(&projected_positions[mi], &centroid_position))
                    .sum();
                total / count as f64
            } else {
                0.0
            };

            let cohesion = 1.0 / (1.0 + angular_spread);

            summaries.push(CategorySummary {
                name: name.clone(),
                member_indices: members.clone(),
                centroid_embedding: centroid_emb,
                centroid_position,
                angular_spread,
                cohesion,
                member_count: count,
                // Spatial fields filled after SpatialQuality is computed
                cap_area: 0.0,
                exclusivity: 0.0,
                voronoi_area: 0.0,
                territorial_efficiency: 0.0,
                // Filled after build_graph() classifies bridges
                bridge_quality: 0.0,
            });
        }

        // 3. Compute spatial quality from category geometry
        let centroids: Vec<SphericalPoint> =
            summaries.iter().map(|s| s.centroid_position).collect();
        let half_angles: Vec<f64> = summaries.iter().map(|s| s.angular_spread).collect();
        let mut spatial_quality =
            SpatialQuality::compute_with_config(&centroids, &half_angles, evr, config);

        // Backfill spatial fields on summaries
        for (i, summary) in summaries.iter_mut().enumerate() {
            summary.cap_area = spatial_quality.cap_areas[i];
            summary.exclusivity = spatial_quality.exclusivities[i];
            summary.voronoi_area = spatial_quality.voronoi_area(i);
            summary.territorial_efficiency =
                spatial_quality.territorial_efficiency(i, summary.member_count);
        }

        // 4. Build category graph + detect bridges (spatially informed)
        let graph = Self::build_graph(&summaries, embeddings, num_cats, &spatial_quality, config);

        // Backfill bridge_quality on summaries using the just-built graph.
        // bridge_quality = mean of (edge.mean_bridge_strength × territorial_factor)
        // over all outgoing edges; 0 for isolated categories.
        for (i, summary) in summaries.iter_mut().enumerate() {
            let edges = &graph.adjacency[i];
            if edges.is_empty() {
                summary.bridge_quality = 0.0;
            } else {
                let total: f64 = edges
                    .iter()
                    .map(|e| {
                        e.mean_bridge_strength * spatial_quality.territorial_factor(i, e.target)
                    })
                    .sum();
                summary.bridge_quality = total / edges.len() as f64;
            }
        }

        // Populate C×C spatial bridge quality matrix on SpatialQuality.
        spatial_quality.set_bridge_quality_matrix(&graph);

        // 5. Build inner spheres for qualifying categories
        let inner_spheres = Self::build_inner_spheres(&summaries, embeddings, projection, config);

        // 6. Detect domain groups + build group-level inner spheres
        // (the v2 default-routing unit). Both are no-ops when fewer than
        // two groups exist.
        let pre_layer = CategoryLayer {
            summaries: summaries.clone(),
            name_to_index: name_to_index.clone(),
            graph: graph.clone(),
            outer_positions: projected_positions.to_vec(),
            inner_spheres: inner_spheres.clone(),
            domain_groups: Vec::new(),
            group_inner_spheres: HashMap::new(),
            spatial_quality: spatial_quality.clone(),
        };
        let domain_groups = crate::domain_groups::detect_domain_groups(
            &pre_layer,
            config.routing.num_domain_groups,
        );
        let group_inner_spheres = Self::build_group_inner_spheres(
            &domain_groups,
            &summaries,
            embeddings,
            projection,
            config,
        );

        CategoryLayer {
            summaries,
            name_to_index,
            graph,
            outer_positions: projected_positions.to_vec(),
            inner_spheres,
            domain_groups,
            group_inner_spheres,
            spatial_quality,
        }
    }

    /// Annotate every bridge with a heuristically inferred [`RelationType`].
    ///
    /// Looks at each bridge item's concept label (indexed into `labels`
    /// by `BridgeItem::item_index`) along with its source and target
    /// category names. Bridges whose label doesn't match any cue remain
    /// `relation: None`.
    pub fn annotate_bridge_relations(&mut self, labels: &[String]) {
        for ((src_cat, tgt_cat), bridges) in self.graph.bridges.iter_mut() {
            let src_name = &self.summaries[*src_cat].name;
            let tgt_name = &self.summaries[*tgt_cat].name;
            for bridge in bridges.iter_mut() {
                if bridge.item_index < labels.len() {
                    bridge.relation =
                        infer_relation(&labels[bridge.item_index], src_name, tgt_name);
                }
            }
        }
    }

    /// Histogram of relation types across every bridge in the graph.
    ///
    /// `None` is included as its own bucket so callers can see how many
    /// bridges remain unlabeled after annotation.
    pub fn relation_census(&self) -> std::collections::HashMap<Option<RelationType>, usize> {
        let mut counts = std::collections::HashMap::new();
        for bridges in self.graph.bridges.values() {
            for b in bridges {
                *counts.entry(b.relation).or_default() += 1;
            }
        }
        counts
    }

    /// Build the inter-category adjacency graph and detect bridge items.
    ///
    /// Bridge detection uses the spatial quality's EVR-adaptive threshold
    /// and exclusivity-weighted strength to prevent bridge inflation at
    /// low EVR and discount bridges between overlapping categories.
    fn build_graph(
        summaries: &[CategorySummary],
        embeddings: &[Embedding],
        num_cats: usize,
        spatial: &SpatialQuality,
        config: &PipelineConfig,
    ) -> CategoryGraph {
        let mut centroid_dists = vec![vec![0.0; num_cats]; num_cats];
        for i in 0..num_cats {
            for j in (i + 1)..num_cats {
                let d = angular_distance(
                    &summaries[i].centroid_position,
                    &summaries[j].centroid_position,
                );
                centroid_dists[i][j] = d;
                centroid_dists[j][i] = d;
            }
        }

        let bridge_threshold = spatial.bridge_threshold;

        let mut bridges: HashMap<(usize, usize), Vec<BridgeItem>> = HashMap::new();

        for (ci, summary) in summaries.iter().enumerate() {
            let centroid_a = &summary.centroid_embedding;

            for &mi in &summary.member_indices {
                let item_emb = &embeddings[mi];
                // Invariant: all embeddings are validated to share the same dimension
                // at pipeline build time, so cosine_similarity cannot fail here.
                let sim_to_own = cosine_similarity(&item_emb.values, centroid_a)
                    .expect("centroid and member share fixed embedding dimensionality");

                for (cj, other_summary) in summaries.iter().enumerate() {
                    if ci == cj {
                        continue;
                    }

                    // Same invariant as sim_to_own above.
                    let sim_to_other =
                        cosine_similarity(&item_emb.values, &other_summary.centroid_embedding)
                            .expect("centroid and member share fixed embedding dimensionality");

                    // EVR-adaptive threshold: stricter when projection is lossy.
                    // Require positive own-affinity too; otherwise the harmonic-mean
                    // strength below collapses or goes negative for a category that
                    // doesn't even hold the item.
                    if sim_to_own > 0.0
                        && sim_to_other > 0.0
                        && sim_to_other > sim_to_own * bridge_threshold
                    {
                        let raw_strength = if sim_to_own + sim_to_other > f64::EPSILON {
                            2.0 * sim_to_own * sim_to_other / (sim_to_own + sim_to_other)
                        } else {
                            0.0
                        };

                        // Discount bridges between heavily overlapping categories.
                        // Two categories that can't be distinguished on S² produce
                        // "bridges" that are really just shared territory.
                        let territorial = spatial.territorial_factor(ci, cj);
                        let bridge_strength = raw_strength * territorial;

                        bridges.entry((ci, cj)).or_default().push(BridgeItem {
                            item_index: mi,
                            source_category: ci,
                            target_category: cj,
                            affinity_to_source: sim_to_own,
                            affinity_to_target: sim_to_other,
                            bridge_strength,
                            // Populated in the classification pass below.
                            classification: BridgeClassification::Weak,
                            relation: None,
                        });
                    }
                }
            }
        }

        // Classification pass: compare each bridge against the corpus-wide
        // median strength and the pair's territorial separation on S².
        //
        // When EVR is below `min_evr_for_classification`, the projection
        // is too lossy for territorial factors to distinguish genuine
        // bridges from overlap artifacts — every factor collapses to
        // near-zero, every bridge gets labeled OverlapArtifact, and the
        // tuner landscape flattens. Skip the territorial check entirely
        // in that regime and leave the default `Weak` label on each
        // bridge (honest uncertainty).
        if spatial.evr >= config.bridges.min_evr_for_classification {
            let mut all_strengths: Vec<f64> = bridges
                .values()
                .flat_map(|list| list.iter().map(|b| b.bridge_strength))
                .collect();
            let median_strength = if all_strengths.is_empty() {
                0.0
            } else {
                all_strengths.sort_by(|a, b| a.total_cmp(b));
                all_strengths[all_strengths.len() / 2]
            };

            let overlap_threshold = config.bridges.overlap_artifact_territorial;
            for list in bridges.values_mut() {
                for b in list.iter_mut() {
                    let tf = spatial.territorial_factor(b.source_category, b.target_category);
                    b.classification = if tf < overlap_threshold {
                        BridgeClassification::OverlapArtifact
                    } else if b.bridge_strength >= median_strength {
                        BridgeClassification::Genuine
                    } else {
                        BridgeClassification::Weak
                    };
                }
            }
        }

        for list in bridges.values_mut() {
            list.sort_by(|a, b| {
                b.bridge_strength
                    .partial_cmp(&a.bridge_strength)
                    .unwrap_or(std::cmp::Ordering::Equal)
            });
        }

        let mut adjacency: Vec<Vec<CategoryEdge>> = vec![Vec::new(); num_cats];
        for i in 0..num_cats {
            for (j, &cd) in centroid_dists[i].iter().enumerate() {
                if i == j {
                    continue;
                }
                let bridge_list = bridges.get(&(i, j));
                let bridge_count = bridge_list.map_or(0, |b| b.len());
                let max_bridge_strength = bridge_list
                    .and_then(|b| b.first().map(|item| item.bridge_strength))
                    .unwrap_or(0.0);
                let mean_bridge_strength = bridge_list
                    .map(|b| {
                        let sum: f64 = b.iter().map(|item| item.bridge_strength).sum();
                        sum / b.len() as f64
                    })
                    .unwrap_or(0.0);

                // Voronoi neighbors get a routing bonus — they're geometrically
                // adjacent on S², so traversal between them is natural even
                // without strong bridges.
                let voronoi_factor = if spatial.are_voronoi_neighbors(i, j) {
                    0.8
                } else {
                    1.0
                };

                let weight =
                    cd * voronoi_factor / (1.0 + bridge_count as f64 * mean_bridge_strength);

                adjacency[i].push(CategoryEdge {
                    target: j,
                    centroid_distance: cd,
                    bridge_count,
                    max_bridge_strength,
                    mean_bridge_strength,
                    weight,
                });
            }
            adjacency[i].sort_by(|a, b| {
                a.weight
                    .partial_cmp(&b.weight)
                    .unwrap_or(std::cmp::Ordering::Equal)
            });
        }

        CategoryGraph { adjacency, bridges }
    }

    /// Evaluate each category and build inner spheres where they help.
    fn build_inner_spheres<P: Projection>(
        summaries: &[CategorySummary],
        embeddings: &[Embedding],
        projection: &P,
        config: &PipelineConfig,
    ) -> HashMap<usize, InnerSphere> {
        let mut result = HashMap::new();
        let cfg = &config.inner_sphere;

        for (ci, summary) in summaries.iter().enumerate() {
            if summary.member_count < cfg.min_size {
                continue;
            }

            let member_embs: Vec<Embedding> = summary
                .member_indices
                .iter()
                .map(|&i| embeddings[i].clone())
                .collect();

            // Global subset EVR: mean certainty under global projection
            let global_subset_evr: f64 = member_embs
                .iter()
                .map(|e| projection.project_rich(e).certainty)
                .sum::<f64>()
                / member_embs.len() as f64;

            // Fit inner linear PCA. On failure (too few members, dim
            // too low, etc.) skip this category's inner sphere
            // silently — the outer sphere still covers queries.
            let Ok(inner_pca) = PcaProjection::fit(&member_embs, RadialStrategy::Fixed(1.0)) else {
                continue;
            };
            let inner_linear_evr = inner_pca.explained_variance_ratio();

            if inner_linear_evr - global_subset_evr < cfg.min_evr_improvement {
                continue;
            }

            let (inner_proj, inner_evr) = if summary.member_count >= cfg.kernel_pca_min_size {
                // Kernel PCA can fail on degenerate subsets too; fall
                // back to the already-successful linear fit.
                match KernelPcaProjection::fit(&member_embs, RadialStrategy::Fixed(1.0)) {
                    Ok(inner_kpca) => {
                        let kernel_evr = inner_kpca.explained_variance_ratio();
                        if kernel_evr > inner_linear_evr + cfg.min_kernel_improvement {
                            (InnerProjection::KernelPca(inner_kpca), kernel_evr)
                        } else {
                            (InnerProjection::LinearPca(inner_pca), inner_linear_evr)
                        }
                    }
                    Err(_) => (InnerProjection::LinearPca(inner_pca), inner_linear_evr),
                }
            } else {
                (InnerProjection::LinearPca(inner_pca), inner_linear_evr)
            };

            let inner_positions: Vec<SphericalPoint> =
                member_embs.iter().map(|e| inner_proj.project(e)).collect();

            result.insert(
                ci,
                InnerSphere {
                    projection: inner_proj,
                    inner_positions,
                    member_indices: summary.member_indices.clone(),
                    explained_variance_ratio: inner_evr,
                    global_subset_evr,
                    evr_improvement: inner_evr - global_subset_evr,
                },
            );
        }

        result
    }

    /// Fit one inner sphere per domain group. Each group's inner
    /// projection sees the union of its member categories' embeddings.
    ///
    /// Built only when (a) ≥ 2 groups exist, (b) the union has at
    /// least [`InnerSphereConfig::min_size`] members, and (c) the
    /// inner-vs-outer EVR improvement clears the same threshold the
    /// per-category builder uses. Kernel PCA is attempted at the same
    /// `kernel_pca_min_size` cutoff as the per-category path.
    fn build_group_inner_spheres<P: Projection>(
        groups: &[crate::domain_groups::DomainGroup],
        summaries: &[CategorySummary],
        embeddings: &[Embedding],
        projection: &P,
        config: &PipelineConfig,
    ) -> HashMap<usize, InnerSphere> {
        let mut result = HashMap::new();
        if groups.len() < 2 {
            return result;
        }
        let cfg = &config.inner_sphere;

        for (gi, group) in groups.iter().enumerate() {
            let mut member_indices: Vec<usize> = Vec::new();
            for &ci in &group.member_categories {
                member_indices.extend(summaries[ci].member_indices.iter().copied());
            }
            if member_indices.len() < cfg.min_size {
                continue;
            }

            let member_embs: Vec<Embedding> = member_indices
                .iter()
                .map(|&i| embeddings[i].clone())
                .collect();

            let global_subset_evr: f64 = member_embs
                .iter()
                .map(|e| projection.project_rich(e).certainty)
                .sum::<f64>()
                / member_embs.len() as f64;

            let Ok(inner_pca) = PcaProjection::fit(&member_embs, RadialStrategy::Fixed(1.0)) else {
                continue;
            };
            let inner_linear_evr = inner_pca.explained_variance_ratio();
            if inner_linear_evr - global_subset_evr < cfg.min_evr_improvement {
                continue;
            }

            let (inner_proj, inner_evr) = if member_indices.len() >= cfg.kernel_pca_min_size {
                match KernelPcaProjection::fit(&member_embs, RadialStrategy::Fixed(1.0)) {
                    Ok(inner_kpca) => {
                        let kernel_evr = inner_kpca.explained_variance_ratio();
                        if kernel_evr > inner_linear_evr + cfg.min_kernel_improvement {
                            (InnerProjection::KernelPca(inner_kpca), kernel_evr)
                        } else {
                            (InnerProjection::LinearPca(inner_pca), inner_linear_evr)
                        }
                    }
                    Err(_) => (InnerProjection::LinearPca(inner_pca), inner_linear_evr),
                }
            } else {
                (InnerProjection::LinearPca(inner_pca), inner_linear_evr)
            };

            let inner_positions: Vec<SphericalPoint> =
                member_embs.iter().map(|e| inner_proj.project(e)).collect();

            result.insert(
                gi,
                InnerSphere {
                    projection: inner_proj,
                    inner_positions,
                    member_indices,
                    explained_variance_ratio: inner_evr,
                    global_subset_evr,
                    evr_improvement: inner_evr - global_subset_evr,
                },
            );
        }

        result
    }

    // ── Phase 1 query methods (unchanged) ──────────────────────────────

    /// Number of categories.
    pub fn num_categories(&self) -> usize {
        self.summaries.len()
    }

    /// Look up a category by name.
    pub fn get_category(&self, name: &str) -> Option<&CategorySummary> {
        self.name_to_index
            .get(name)
            .map(|&idx| &self.summaries[idx])
    }

    /// Get the k nearest neighbor categories to the given category.
    pub fn category_neighbors(&self, category_name: &str, k: usize) -> Vec<&CategorySummary> {
        let Some(&ci) = self.name_to_index.get(category_name) else {
            return Vec::new();
        };
        self.graph.adjacency[ci]
            .iter()
            .take(k)
            .map(|edge| &self.summaries[edge.target])
            .collect()
    }

    /// Get bridge items between two categories.
    pub fn bridge_items(
        &self,
        source_category: &str,
        target_category: &str,
        max_bridges: usize,
    ) -> Vec<&BridgeItem> {
        let Some(&si) = self.name_to_index.get(source_category) else {
            return Vec::new();
        };
        let Some(&ti) = self.name_to_index.get(target_category) else {
            return Vec::new();
        };
        self.graph
            .bridges
            .get(&(si, ti))
            .map(|list| list.iter().take(max_bridges).collect())
            .unwrap_or_default()
    }

    /// Find the shortest path between two categories through the category graph.
    pub fn category_path(
        &self,
        source_category: &str,
        target_category: &str,
    ) -> Option<CategoryPath> {
        let &si = self.name_to_index.get(source_category)?;
        let &ti = self.name_to_index.get(target_category)?;
        if si == ti {
            return Some(CategoryPath {
                steps: vec![CategoryPathStep {
                    category_index: si,
                    category_name: self.summaries[si].name.clone(),
                    cumulative_distance: 0.0,
                    bridges_to_next: Vec::new(),
                    hop_confidence: 0.0,
                }],
                total_distance: 0.0,
                path_confidence: 1.0,
            });
        }

        // Dijkstra via binary-heap. Previously a linear scan over
        // `dist` per pop gave O(C²) — fine for tens of categories,
        // sloppy under the larger corpora the tuner now exercises.
        // Match the pattern already used in query.rs::concept_path.
        let n = self.summaries.len();
        let mut dist = vec![f64::INFINITY; n];
        let mut prev: Vec<Option<usize>> = vec![None; n];
        let mut heap = std::collections::BinaryHeap::new();

        dist[si] = 0.0;
        heap.push(HeapEntry {
            dist: 0.0,
            node: si,
        });

        while let Some(HeapEntry { dist: d, node: u }) = heap.pop() {
            if u == ti {
                break;
            }
            // Stale entry — we already relaxed this node to something
            // shorter. Skip without touching neighbors.
            if d > dist[u] {
                continue;
            }
            for edge in &self.graph.adjacency[u] {
                let nd = d + edge.weight;
                if nd < dist[edge.target] {
                    dist[edge.target] = nd;
                    prev[edge.target] = Some(u);
                    heap.push(HeapEntry {
                        dist: nd,
                        node: edge.target,
                    });
                }
            }
        }

        if dist[ti].is_infinite() {
            return None;
        }

        let mut path_indices = Vec::new();
        let mut cur = ti;
        loop {
            path_indices.push(cur);
            match prev[cur] {
                Some(p) => cur = p,
                None => break,
            }
        }
        path_indices.reverse();

        let mut steps = Vec::with_capacity(path_indices.len());
        for (step_idx, &ci) in path_indices.iter().enumerate() {
            let bridges_to_next = if step_idx + 1 < path_indices.len() {
                let next_ci = path_indices[step_idx + 1];
                self.graph
                    .bridges
                    .get(&(ci, next_ci))
                    .map(|list| list.iter().take(3).cloned().collect())
                    .unwrap_or_default()
            } else {
                Vec::new()
            };

            let hop_confidence = if step_idx + 1 < path_indices.len() {
                let next_ci = path_indices[step_idx + 1];
                let edge_strength = self.graph.adjacency[ci]
                    .iter()
                    .find(|e| e.target == next_ci)
                    .map_or(0.0, |e| e.max_bridge_strength);
                let territorial = self.spatial_quality.territorial_factor(ci, next_ci);
                let voronoi_bonus = if self.spatial_quality.are_voronoi_neighbors(ci, next_ci) {
                    1.2
                } else {
                    1.0
                };
                (edge_strength * territorial * voronoi_bonus).min(1.0)
            } else {
                0.0
            };

            steps.push(CategoryPathStep {
                category_index: ci,
                category_name: self.summaries[ci].name.clone(),
                cumulative_distance: dist[ci],
                bridges_to_next,
                hop_confidence,
            });
        }

        let path_confidence = steps
            .iter()
            .take(steps.len().saturating_sub(1))
            .map(|s| s.hop_confidence)
            .fold(1.0, |acc, c| acc * c.max(0.01));

        Some(CategoryPath {
            total_distance: dist[ti],
            steps,
            path_confidence,
        })
    }

    /// Find all categories whose centroid is within `max_angle` radians
    /// of the given embedding's projected position.
    pub fn categories_near_embedding<P: Projection>(
        &self,
        embedding: &Embedding,
        projection: &P,
        max_angle: f64,
    ) -> Vec<(usize, f64)> {
        let pos = projection.project(embedding);
        let mut results: Vec<(usize, f64)> = self
            .summaries
            .iter()
            .enumerate()
            .map(|(i, s)| (i, angular_distance(&pos, &s.centroid_position)))
            .filter(|&(_, d)| d <= max_angle)
            .collect();
        results.sort_by(|a, b| a.1.total_cmp(&b.1));
        results
    }

    /// Certainty-weighted category routing.
    ///
    /// Like [`Self::categories_near_embedding`] but penalizes routes through
    /// low-certainty projection regions. The effective distance is scaled
    /// by `1 / sqrt(certainty)`, so poorly-projected queries don't get
    /// routed to whatever random centroid happens to be angularly close
    /// in the distorted projection.
    ///
    /// Returns `(category_index, raw_angular_distance, effective_distance, certainty)`.
    pub fn categories_near_embedding_weighted<P: Projection>(
        &self,
        embedding: &Embedding,
        projection: &P,
        max_angle: f64,
    ) -> Vec<(usize, f64, f64, f64)> {
        let rich = projection.project_rich(embedding);
        let pos = rich.position;
        let certainty = rich.certainty.max(0.001);

        let mut results: Vec<(usize, f64, f64, f64)> = self
            .summaries
            .iter()
            .enumerate()
            .map(|(i, s)| {
                let raw_dist = angular_distance(&pos, &s.centroid_position);
                // Low certainty inflates distance → avoids noisy routing
                let effective = raw_dist / certainty.sqrt();
                (i, raw_dist, effective, certainty)
            })
            .filter(|&(_, raw, _, _)| raw <= max_angle)
            .collect();
        results.sort_by(|a, b| a.2.total_cmp(&b.2));
        results
    }

    // ── Phase 2 query methods (inner spheres) ──────────────────────────

    /// Whether a given category has an inner sphere.
    pub fn has_inner_sphere(&self, category_name: &str) -> bool {
        self.name_to_index
            .get(category_name)
            .is_some_and(|&ci| self.inner_spheres.contains_key(&ci))
    }

    /// Get the inner sphere for a category, if one exists.
    pub fn get_inner_sphere(&self, category_name: &str) -> Option<&InnerSphere> {
        self.name_to_index
            .get(category_name)
            .and_then(|&ci| self.inner_spheres.get(&ci))
    }

    /// Number of categories that have inner spheres.
    pub fn num_inner_spheres(&self) -> usize {
        self.inner_spheres.len()
    }

    /// Drill down with an explicit outer projection for the fallback case.
    ///
    /// When no inner sphere exists, the query is projected using the
    /// provided projection and compared against stored outer positions.
    pub fn drill_down_with_projection<P: Projection>(
        &self,
        category_name: &str,
        embedding: &Embedding,
        projection: &P,
        k: usize,
    ) -> Vec<DrillDownResult> {
        let Some(&ci) = self.name_to_index.get(category_name) else {
            return Vec::new();
        };
        let summary = &self.summaries[ci];

        if let Some(inner) = self.inner_spheres.get(&ci) {
            let query_pos = inner.projection.project(embedding);
            let mut results: Vec<DrillDownResult> = inner
                .inner_positions
                .iter()
                .enumerate()
                .map(|(local_idx, pos)| DrillDownResult {
                    item_index: inner.member_indices[local_idx],
                    distance: angular_distance(&query_pos, pos),
                    used_inner_sphere: true,
                })
                .collect();
            results.sort_by(|a, b| {
                a.distance
                    .partial_cmp(&b.distance)
                    .unwrap_or(std::cmp::Ordering::Equal)
            });
            results.truncate(k);
            results
        } else {
            let query_pos = projection.project(embedding);
            let mut results: Vec<DrillDownResult> = summary
                .member_indices
                .iter()
                .map(|&mi| DrillDownResult {
                    item_index: mi,
                    distance: angular_distance(&self.outer_positions[mi], &query_pos),
                    used_inner_sphere: false,
                })
                .collect();
            results.sort_by(|a, b| {
                a.distance
                    .partial_cmp(&b.distance)
                    .unwrap_or(std::cmp::Ordering::Equal)
            });
            results.truncate(k);
            results
        }
    }

    /// Drill down into a domain group's inner sphere. Returns an empty
    /// vec if the group index is unknown or has no inner sphere.
    ///
    /// Used by the v2 default `nearest()` path: route to the closest
    /// group, then run k-NN against the group-level inner positions.
    pub fn drill_down_group(
        &self,
        group_index: usize,
        embedding: &Embedding,
        k: usize,
    ) -> Vec<DrillDownResult> {
        let Some(inner) = self.group_inner_spheres.get(&group_index) else {
            return Vec::new();
        };
        let query_pos = inner.projection.project(embedding);
        let mut results: Vec<DrillDownResult> = inner
            .inner_positions
            .iter()
            .enumerate()
            .map(|(local_idx, pos)| DrillDownResult {
                item_index: inner.member_indices[local_idx],
                distance: angular_distance(&query_pos, pos),
                used_inner_sphere: true,
            })
            .collect();
        results.sort_by(|a, b| {
            a.distance
                .partial_cmp(&b.distance)
                .unwrap_or(std::cmp::Ordering::Equal)
        });
        results.truncate(k);
        results
    }

    /// Index of the domain group whose centroid is nearest to `pos`,
    /// along with the angular distance to the *second*-nearest group's
    /// centroid (or `f64::INFINITY` if there is only one group).
    ///
    /// The default `nearest()` routing in
    /// [`SphereQLPipeline`](crate::pipeline::SphereQLPipeline) drills
    /// into the nearest group's inner sphere only when
    /// `nearest_dist / second_dist < group_routing_alpha`; the second
    /// distance is returned alongside so the caller can apply that
    /// gate without recomputing.
    pub fn nearest_group(&self, pos: &SphericalPoint) -> Option<(usize, f64, f64)> {
        if self.domain_groups.is_empty() {
            return None;
        }
        let mut best = (0usize, f64::INFINITY);
        let mut second = f64::INFINITY;
        for (gi, g) in self.domain_groups.iter().enumerate() {
            let d = angular_distance(pos, &g.centroid);
            if d < best.1 {
                second = best.1;
                best = (gi, d);
            } else if d < second {
                second = d;
            }
        }
        Some((best.0, best.1, second))
    }

    /// Report which categories have inner spheres, their projection type,
    /// and EVR metrics.
    pub fn inner_sphere_stats(&self) -> Vec<InnerSphereReport> {
        let mut reports: Vec<InnerSphereReport> = self
            .inner_spheres
            .iter()
            .map(|(&ci, inner)| {
                let proj_type = match &inner.projection {
                    InnerProjection::LinearPca(_) => "LinearPca",
                    InnerProjection::KernelPca(_) => "KernelPca",
                };
                InnerSphereReport {
                    category_name: self.summaries[ci].name.clone(),
                    category_index: ci,
                    member_count: inner.member_indices.len(),
                    projection_type: proj_type,
                    inner_evr: inner.explained_variance_ratio,
                    global_subset_evr: inner.global_subset_evr,
                    evr_improvement: inner.evr_improvement,
                }
            })
            .collect();
        reports.sort_by_key(|r| r.category_index);
        reports
    }
}

// ── Helpers ────────────────────────────────────────────────────────────

/// Min-heap entry keyed on `dist` for [`CategoryLayer::category_path`]'s
/// Dijkstra. `BinaryHeap` is a max-heap, so `Ord` is reversed.
/// `total_cmp` is NaN-safe (NaN sorts to the end).
#[derive(PartialEq)]
struct HeapEntry {
    dist: f64,
    node: usize,
}

impl Eq for HeapEntry {}

impl PartialOrd for HeapEntry {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for HeapEntry {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        other.dist.total_cmp(&self.dist)
    }
}

// ── Tests ──────────────────────────────────────────────────────────────

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

    fn emb(vals: &[f64]) -> Embedding {
        Embedding::new(vals.to_vec())
    }

    // --- Phase 1 test helpers ---

    fn test_corpus() -> (Vec<String>, Vec<Embedding>) {
        let categories = vec![
            "science".into(),
            "science".into(),
            "science".into(),
            "science".into(),
            "cooking".into(),
            "cooking".into(),
            "cooking".into(),
            "cooking".into(),
            "music".into(),
            "music".into(),
            "music".into(),
            "music".into(),
        ];
        let embeddings = vec![
            emb(&[1.0, 0.1, 0.0, 0.05, 0.02]),
            emb(&[0.9, 0.15, 0.05, 0.03, 0.01]),
            emb(&[0.95, 0.05, 0.1, 0.04, 0.03]),
            emb(&[0.85, 0.2, 0.0, 0.06, 0.01]),
            emb(&[0.1, 1.0, 0.0, 0.02, 0.05]),
            emb(&[0.15, 0.9, 0.05, 0.03, 0.04]),
            emb(&[0.05, 0.95, 0.1, 0.01, 0.06]),
            emb(&[0.2, 0.85, 0.0, 0.04, 0.03]),
            emb(&[0.0, 0.1, 1.0, 0.05, 0.02]),
            emb(&[0.05, 0.15, 0.9, 0.03, 0.01]),
            emb(&[0.1, 0.05, 0.95, 0.04, 0.03]),
            emb(&[0.0, 0.2, 0.85, 0.06, 0.01]),
        ];
        (categories, embeddings)
    }

    fn build_test_layer() -> (CategoryLayer, Vec<Embedding>, PcaProjection) {
        let (categories, embeddings) = test_corpus();
        let pca = PcaProjection::fit(&embeddings, RadialStrategy::Fixed(1.0)).unwrap();
        let projected: Vec<SphericalPoint> = embeddings.iter().map(|e| pca.project(e)).collect();
        let evr = pca.explained_variance_ratio();
        let layer = CategoryLayer::build(&categories, &embeddings, &projected, &pca, evr);
        (layer, embeddings, pca)
    }

    // --- Phase 2 test helpers ---

    fn large_category_corpus() -> (Vec<String>, Vec<Embedding>) {
        let mut categories = Vec::new();
        let mut embeddings = Vec::new();

        for i in 0..25 {
            categories.push("big".into());
            let t = i as f64 / 25.0;
            let mut v = vec![0.0; 10];
            v[0] = 1.0 + 0.3 * (t * std::f64::consts::TAU).sin();
            v[1] = 0.5 + 0.3 * (t * std::f64::consts::TAU).cos();
            v[2] = 0.2 * t;
            for (d, slot) in v.iter_mut().enumerate().take(10).skip(3) {
                *slot = 0.01 * ((i * 7 + d) as f64 % 1.0);
            }
            embeddings.push(emb(&v));
        }

        for i in 0..4 {
            categories.push("small_a".into());
            let mut v = vec![0.0; 10];
            v[5] = 1.0 + 0.1 * i as f64;
            v[6] = 0.05;
            embeddings.push(emb(&v));
        }

        for i in 0..4 {
            categories.push("small_b".into());
            let mut v = vec![0.0; 10];
            v[8] = 1.0 + 0.1 * i as f64;
            v[9] = 0.05;
            embeddings.push(emb(&v));
        }

        (categories, embeddings)
    }

    fn build_large_test_layer() -> (CategoryLayer, Vec<Embedding>, PcaProjection) {
        let (categories, embeddings) = large_category_corpus();
        let pca = PcaProjection::fit(&embeddings, RadialStrategy::Fixed(1.0)).unwrap();
        let projected: Vec<SphericalPoint> = embeddings.iter().map(|e| pca.project(e)).collect();
        let evr = pca.explained_variance_ratio();
        let layer = CategoryLayer::build(&categories, &embeddings, &projected, &pca, evr);
        (layer, embeddings, pca)
    }

    // ======== Phase 1 tests (unchanged) ========

    #[test]
    fn builds_correct_number_of_categories() {
        let (layer, _, _) = build_test_layer();
        assert_eq!(layer.num_categories(), 3);
    }

    #[test]
    fn category_names_correct() {
        let (layer, _, _) = build_test_layer();
        let names: Vec<&str> = layer.summaries.iter().map(|s| s.name.as_str()).collect();
        assert!(names.contains(&"science"));
        assert!(names.contains(&"cooking"));
        assert!(names.contains(&"music"));
    }

    #[test]
    fn member_counts_correct() {
        let (layer, _, _) = build_test_layer();
        for summary in &layer.summaries {
            assert_eq!(summary.member_count, 4);
            assert_eq!(summary.member_indices.len(), 4);
        }
    }

    #[test]
    fn centroid_embedding_is_mean() {
        let (layer, embeddings, _) = build_test_layer();
        let science = layer.get_category("science").unwrap();
        let mut expected = vec![0.0; 5];
        for emb in embeddings.iter().take(4) {
            for (j, &v) in emb.values.iter().enumerate() {
                expected[j] += v;
            }
        }
        for v in &mut expected {
            *v /= 4.0;
        }
        for (j, (&actual, &exp)) in science
            .centroid_embedding
            .iter()
            .zip(expected.iter())
            .enumerate()
        {
            assert!(
                (actual - exp).abs() < 1e-10,
                "centroid dim {j}: {actual} != {exp}"
            );
        }
    }

    #[test]
    fn angular_spread_is_nonnegative() {
        let (layer, _, _) = build_test_layer();
        for s in &layer.summaries {
            assert!(s.angular_spread >= 0.0);
        }
    }

    #[test]
    fn cohesion_in_range() {
        let (layer, _, _) = build_test_layer();
        for s in &layer.summaries {
            assert!(s.cohesion > 0.0 && s.cohesion <= 1.0);
        }
    }

    #[test]
    fn graph_has_edges_for_all_pairs() {
        let (layer, _, _) = build_test_layer();
        for (i, edges) in layer.graph.adjacency.iter().enumerate() {
            assert_eq!(edges.len(), layer.num_categories() - 1, "cat {i}");
        }
    }

    #[test]
    fn edge_weights_positive() {
        let (layer, _, _) = build_test_layer();
        for edges in &layer.graph.adjacency {
            for e in edges {
                assert!(e.weight > 0.0);
                assert!(e.centroid_distance > 0.0);
            }
        }
    }

    #[test]
    fn edges_sorted_by_weight() {
        let (layer, _, _) = build_test_layer();
        for edges in &layer.graph.adjacency {
            for w in edges.windows(2) {
                assert!(w[0].weight <= w[1].weight);
            }
        }
    }

    #[test]
    fn edge_bridge_strength_fields_populated() {
        let (layer, _, _) = build_test_layer();
        for edges in &layer.graph.adjacency {
            for e in edges {
                assert!(e.max_bridge_strength >= 0.0 && e.max_bridge_strength <= 1.0);
                assert!(e.mean_bridge_strength >= 0.0 && e.mean_bridge_strength <= 1.0);
                assert!(e.mean_bridge_strength <= e.max_bridge_strength + 1e-10);
                if e.bridge_count > 0 {
                    assert!(e.max_bridge_strength > 0.0);
                    assert!(e.mean_bridge_strength > 0.0);
                } else {
                    assert!(e.max_bridge_strength == 0.0);
                    assert!(e.mean_bridge_strength == 0.0);
                }
            }
        }
    }

    #[test]
    fn edge_weight_incorporates_bridge_strength() {
        let (layer, _, _) = build_test_layer();
        for edges in &layer.graph.adjacency {
            for e in edges {
                // The Voronoi bonus can only reduce weight, so weight ≤ base formula.
                // The territorial factor can also reduce bridge strength.
                let base_no_bonus =
                    e.centroid_distance / (1.0 + e.bridge_count as f64 * e.mean_bridge_strength);
                assert!(
                    e.weight <= base_no_bonus + 1e-10,
                    "weight {:.6} should be ≤ base {:.6} (Voronoi bonus reduces it)",
                    e.weight,
                    base_no_bonus,
                );
                assert!(e.weight > 0.0, "weight must be positive");
            }
        }
    }

    #[test]
    fn get_category_by_name() {
        let (layer, _, _) = build_test_layer();
        assert!(layer.get_category("science").is_some());
        assert!(layer.get_category("astrology").is_none());
    }

    #[test]
    fn category_neighbors_returns_sorted() {
        let (layer, _, _) = build_test_layer();
        assert_eq!(layer.category_neighbors("science", 2).len(), 2);
    }

    #[test]
    fn category_neighbors_k_larger_than_available() {
        let (layer, _, _) = build_test_layer();
        assert_eq!(layer.category_neighbors("science", 100).len(), 2);
    }

    #[test]
    fn category_neighbors_unknown_returns_empty() {
        let (layer, _, _) = build_test_layer();
        assert!(layer.category_neighbors("nonexistent", 5).is_empty());
    }

    #[test]
    fn bridge_items_detected() {
        let (layer, _, _) = build_test_layer();
        let _ = layer.bridge_items("science", "cooking", 10);
    }

    #[test]
    fn bridge_items_unknown_category_returns_empty() {
        let (layer, _, _) = build_test_layer();
        assert!(layer.bridge_items("science", "nonexistent", 10).is_empty());
    }

    #[test]
    fn bridge_classification_populated() {
        let (layer, _, _) = build_test_layer();
        for bridges in layer.graph.bridges.values() {
            for b in bridges {
                // Every bridge should have one of the three classifications.
                assert!(
                    b.classification == BridgeClassification::Genuine
                        || b.classification == BridgeClassification::OverlapArtifact
                        || b.classification == BridgeClassification::Weak
                );
            }
        }
    }

    #[test]
    fn low_evr_skips_territorial_classification() {
        // When `spatial.evr` falls below `min_evr_for_classification`,
        // every bridge should fall back to `Weak`. We force the gate
        // by raising the threshold above the measured EVR rather than
        // synthesizing a low-EVR projection.
        let (categories, embeddings) = test_corpus();
        let pca = PcaProjection::fit(&embeddings, RadialStrategy::Fixed(1.0)).unwrap();
        let projected: Vec<SphericalPoint> = embeddings.iter().map(|e| pca.project(e)).collect();

        let mut config = PipelineConfig::default();
        // Set the gate above the natural ceiling so the early-skip
        // branch is exercised regardless of corpus EVR.
        config.bridges.min_evr_for_classification = 1.5;

        let layer = CategoryLayer::build_with_config(
            &categories,
            &embeddings,
            &projected,
            &pca,
            0.10,
            &config,
        );

        for bridges in layer.graph.bridges.values() {
            for b in bridges {
                assert_eq!(
                    b.classification,
                    BridgeClassification::Weak,
                    "low EVR should label every bridge Weak, got {:?}",
                    b.classification
                );
            }
        }
    }

    #[test]
    fn bridge_quality_nonnegative() {
        let (layer, _, _) = build_test_layer();
        for s in &layer.summaries {
            assert!(
                s.bridge_quality >= 0.0,
                "{} has negative bridge_quality",
                s.name
            );
        }
    }

    #[test]
    fn bridge_quality_matrix_symmetric_ish() {
        let (layer, _, _) = build_test_layer();
        let m = &layer.spatial_quality.bridge_quality_matrix;
        let n = m.len();
        assert_eq!(n, layer.num_categories());
        for (i, row) in m.iter().enumerate() {
            assert_eq!(row.len(), n);
            assert_eq!(row[i], 0.0, "diagonal should be zero");
        }
    }

    #[test]
    fn bridge_strength_in_valid_range() {
        let (layer, _, _) = build_test_layer();
        for list in layer.graph.bridges.values() {
            for b in list {
                assert!(b.bridge_strength >= 0.0 && b.bridge_strength <= 1.0);
            }
        }
    }

    #[test]
    fn bridges_sorted_by_strength() {
        let (layer, _, _) = build_test_layer();
        for list in layer.graph.bridges.values() {
            for w in list.windows(2) {
                assert!(w[0].bridge_strength >= w[1].bridge_strength);
            }
        }
    }

    #[test]
    fn category_path_same_category() {
        let (layer, _, _) = build_test_layer();
        let path = layer.category_path("science", "science").unwrap();
        assert_eq!(path.steps.len(), 1);
        assert!(path.total_distance.abs() < 1e-12);
    }

    #[test]
    fn category_path_adjacent() {
        let (layer, _, _) = build_test_layer();
        let path = layer.category_path("science", "cooking").unwrap();
        assert!(path.steps.len() >= 2);
        assert_eq!(path.steps.first().unwrap().category_name, "science");
        assert_eq!(path.steps.last().unwrap().category_name, "cooking");
        assert!(path.total_distance > 0.0);
    }

    #[test]
    fn category_path_unknown_returns_none() {
        let (layer, _, _) = build_test_layer();
        assert!(layer.category_path("science", "nonexistent").is_none());
    }

    #[test]
    fn category_path_distances_monotonic() {
        let (layer, _, _) = build_test_layer();
        let path = layer.category_path("science", "music").unwrap();
        for w in path.steps.windows(2) {
            assert!(w[1].cumulative_distance >= w[0].cumulative_distance);
        }
    }

    #[test]
    fn categories_near_embedding_finds_correct() {
        let (layer, _, pca) = build_test_layer();
        let near = layer.categories_near_embedding(
            &emb(&[1.0, 0.0, 0.0, 0.0, 0.0]),
            &pca,
            std::f64::consts::PI,
        );
        assert!(!near.is_empty());
        assert_eq!(layer.summaries[near[0].0].name, "science");
    }

    #[test]
    fn categories_near_embedding_sorted_by_distance() {
        let (layer, _, pca) = build_test_layer();
        let near = layer.categories_near_embedding(
            &emb(&[0.5, 0.5, 0.5, 0.0, 0.0]),
            &pca,
            std::f64::consts::PI,
        );
        for w in near.windows(2) {
            assert!(w[0].1 <= w[1].1);
        }
    }

    #[test]
    fn categories_near_embedding_respects_threshold() {
        let (layer, _, pca) = build_test_layer();
        let near = layer.categories_near_embedding(&emb(&[1.0, 0.0, 0.0, 0.0, 0.0]), &pca, 0.01);
        for &(_, d) in &near {
            assert!(d <= 0.01);
        }
    }

    #[test]
    fn cosine_similarity_identical() {
        assert!(
            (cosine_similarity(&[1.0, 0.0, 0.0], &[1.0, 0.0, 0.0]).unwrap() - 1.0).abs() < 1e-12
        );
    }

    #[test]
    fn cosine_similarity_orthogonal() {
        assert!(
            cosine_similarity(&[1.0, 0.0, 0.0], &[0.0, 1.0, 0.0])
                .unwrap()
                .abs()
                < 1e-12
        );
    }

    #[test]
    fn cosine_similarity_opposite() {
        assert!(
            (cosine_similarity(&[1.0, 0.0, 0.0], &[-1.0, 0.0, 0.0]).unwrap() + 1.0).abs() < 1e-12
        );
    }

    #[test]
    fn cosine_similarity_zero_vector() {
        assert!(
            cosine_similarity(&[0.0, 0.0, 0.0], &[1.0, 0.0, 0.0])
                .unwrap()
                .abs()
                < 1e-12
        );
    }

    // ======== Phase 2 tests (inner spheres) ========

    #[test]
    fn small_categories_get_no_inner_sphere() {
        let (layer, _, _) = build_test_layer();
        assert_eq!(layer.num_inner_spheres(), 0);
        assert!(!layer.has_inner_sphere("science"));
    }

    #[test]
    fn large_category_may_get_inner_sphere() {
        let (layer, _, _) = build_large_test_layer();
        assert!(!layer.has_inner_sphere("small_a"));
        assert!(!layer.has_inner_sphere("small_b"));
        let _ = layer.has_inner_sphere("big");
    }

    #[test]
    fn inner_sphere_stats_count_matches() {
        let (layer, _, _) = build_large_test_layer();
        assert_eq!(layer.inner_sphere_stats().len(), layer.num_inner_spheres());
    }

    #[test]
    fn inner_sphere_stats_sorted_by_index() {
        let (layer, _, _) = build_large_test_layer();
        let stats = layer.inner_sphere_stats();
        for w in stats.windows(2) {
            assert!(w[0].category_index <= w[1].category_index);
        }
    }

    #[test]
    fn inner_sphere_evr_improvement_positive() {
        let (layer, _, _) = build_large_test_layer();
        let min_improvement = PipelineConfig::default().inner_sphere.min_evr_improvement;
        for inner in layer.inner_spheres.values() {
            assert!(inner.evr_improvement >= min_improvement);
        }
    }

    #[test]
    fn inner_sphere_positions_match_member_count() {
        let (layer, _, _) = build_large_test_layer();
        for (&ci, inner) in &layer.inner_spheres {
            assert_eq!(inner.inner_positions.len(), inner.member_indices.len());
            assert_eq!(inner.member_indices.len(), layer.summaries[ci].member_count);
        }
    }

    #[test]
    fn inner_sphere_member_indices_valid() {
        let (layer, _, _) = build_large_test_layer();
        let total = layer.outer_positions.len();
        for inner in layer.inner_spheres.values() {
            for &mi in &inner.member_indices {
                assert!(mi < total);
            }
        }
    }

    #[test]
    fn inner_sphere_report_projection_type_valid() {
        let (layer, _, _) = build_large_test_layer();
        for r in layer.inner_sphere_stats() {
            assert!(r.projection_type == "LinearPca" || r.projection_type == "KernelPca");
        }
    }

    #[test]
    fn inner_sphere_evr_in_range() {
        let (layer, _, _) = build_large_test_layer();
        for inner in layer.inner_spheres.values() {
            assert!(inner.explained_variance_ratio >= 0.0 && inner.explained_variance_ratio <= 1.0);
            assert!(inner.global_subset_evr >= 0.0 && inner.global_subset_evr <= 1.0);
        }
    }

    #[test]
    fn has_inner_sphere_unknown_category() {
        let (layer, _, _) = build_test_layer();
        assert!(!layer.has_inner_sphere("nonexistent"));
    }

    #[test]
    fn get_inner_sphere_returns_none_for_small() {
        let (layer, _, _) = build_test_layer();
        assert!(layer.get_inner_sphere("science").is_none());
    }

    #[test]
    fn drill_down_returns_results() {
        let (layer, _, pca) = build_large_test_layer();
        let q = emb(&[1.0, 0.5, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
        let results = layer.drill_down_with_projection("big", &q, &pca, 5);
        assert!(!results.is_empty());
        assert!(results.len() <= 5);
    }

    #[test]
    fn drill_down_sorted_by_distance() {
        let (layer, _, pca) = build_large_test_layer();
        let q = emb(&[1.0, 0.5, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
        let results = layer.drill_down_with_projection("big", &q, &pca, 10);
        for w in results.windows(2) {
            assert!(w[0].distance <= w[1].distance);
        }
    }

    #[test]
    fn drill_down_unknown_category_empty() {
        let (layer, _, pca) = build_large_test_layer();
        assert!(
            layer
                .drill_down_with_projection("nonexistent", &emb(&[1.0; 10]), &pca, 5)
                .is_empty()
        );
    }

    #[test]
    fn drill_down_item_indices_valid() {
        let (layer, _, pca) = build_large_test_layer();
        let q = emb(&[1.0, 0.5, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
        let total = layer.outer_positions.len();
        for r in layer.drill_down_with_projection("big", &q, &pca, 25) {
            assert!(r.item_index < total);
        }
    }

    #[test]
    fn drill_down_small_category_uses_outer() {
        let (layer, _, pca) = build_large_test_layer();
        let q = emb(&[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0]);
        for r in layer.drill_down_with_projection("small_a", &q, &pca, 4) {
            assert!(!r.used_inner_sphere);
        }
    }

    #[test]
    fn drill_down_distances_nonnegative() {
        let (layer, _, pca) = build_large_test_layer();
        for r in layer.drill_down_with_projection("big", &emb(&[1.0; 10]), &pca, 10) {
            assert!(r.distance >= 0.0);
        }
    }

    #[test]
    fn inner_projection_enum_debug() {
        let corpus: Vec<Embedding> = (0..5)
            .map(|i| emb(&[i as f64, 0.0, 0.0, 0.0, 0.0]))
            .collect();
        let pca = PcaProjection::fit(&corpus, RadialStrategy::Fixed(1.0)).unwrap();
        assert_eq!(
            format!("{:?}", InnerProjection::LinearPca(pca)),
            "LinearPca"
        );
    }

    #[test]
    fn inner_projection_projects_correctly() {
        let corpus: Vec<Embedding> = (0..5)
            .map(|i| emb(&[i as f64, 0.0, 0.0, 0.0, 0.0]))
            .collect();
        let pca = PcaProjection::fit(&corpus, RadialStrategy::Fixed(1.0)).unwrap();
        let proj = InnerProjection::LinearPca(pca.clone());
        let e = emb(&[1.0, 0.0, 0.0, 0.0, 0.0]);
        let sp_enum = proj.project(&e);
        let sp_direct = pca.project(&e);
        assert!((sp_enum.theta - sp_direct.theta).abs() < 1e-12);
        assert!((sp_enum.phi - sp_direct.phi).abs() < 1e-12);
    }

    #[test]
    fn inner_projection_dimensionality() {
        let corpus: Vec<Embedding> = (0..5)
            .map(|i| emb(&[i as f64, 0.0, 0.0, 0.0, 0.0]))
            .collect();
        let pca = PcaProjection::fit(&corpus, RadialStrategy::Fixed(1.0)).unwrap();
        assert_eq!(InnerProjection::LinearPca(pca).dimensionality(), 5);
    }

    // ======== min_category_size filtering ========

    #[test]
    fn min_category_size_excludes_small_categories() {
        let (categories, embeddings) = test_corpus(); // 3 categories × 4 items each
        let pca = PcaProjection::fit(&embeddings, RadialStrategy::Fixed(1.0)).unwrap();
        let projected: Vec<SphericalPoint> = embeddings.iter().map(|e| pca.project(e)).collect();
        let evr = pca.explained_variance_ratio();

        let config = PipelineConfig {
            min_category_size: 5,
            ..Default::default()
        };
        let layer = CategoryLayer::build_with_config(
            &categories,
            &embeddings,
            &projected,
            &pca,
            evr,
            &config,
        );

        assert_eq!(
            layer.num_categories(),
            0,
            "all categories below threshold should be excluded"
        );
        assert!(layer.graph.adjacency.is_empty());
        assert!(layer.graph.bridges.is_empty());
    }

    #[test]
    fn min_category_size_one_includes_everything() {
        let (categories, embeddings) = test_corpus();
        let pca = PcaProjection::fit(&embeddings, RadialStrategy::Fixed(1.0)).unwrap();
        let projected: Vec<SphericalPoint> = embeddings.iter().map(|e| pca.project(e)).collect();
        let evr = pca.explained_variance_ratio();

        let config = PipelineConfig {
            min_category_size: 1,
            ..Default::default()
        };
        let layer = CategoryLayer::build_with_config(
            &categories,
            &embeddings,
            &projected,
            &pca,
            evr,
            &config,
        );

        assert_eq!(
            layer.num_categories(),
            3,
            "all categories should be included at min_size=1"
        );
    }

    #[test]
    fn min_category_size_partial_filter() {
        // One large category (10 items) and two small ones (2 each).
        let mut categories = Vec::new();
        let mut embeddings = Vec::new();
        for i in 0..10 {
            categories.push("big".to_string());
            let mut v = vec![0.0; 5];
            v[0] = 1.0 + i as f64 * 0.01;
            embeddings.push(emb(&v));
        }
        for i in 0..2 {
            categories.push("small_a".to_string());
            let mut v = vec![0.0; 5];
            v[1] = 1.0 + i as f64 * 0.01;
            embeddings.push(emb(&v));
        }
        for i in 0..2 {
            categories.push("small_b".to_string());
            let mut v = vec![0.0; 5];
            v[2] = 1.0 + i as f64 * 0.01;
            embeddings.push(emb(&v));
        }

        let pca = PcaProjection::fit(&embeddings, RadialStrategy::Fixed(1.0)).unwrap();
        let projected: Vec<SphericalPoint> = embeddings.iter().map(|e| pca.project(e)).collect();
        let evr = pca.explained_variance_ratio();

        let config = PipelineConfig {
            min_category_size: 5,
            ..Default::default()
        };
        let layer = CategoryLayer::build_with_config(
            &categories,
            &embeddings,
            &projected,
            &pca,
            evr,
            &config,
        );

        assert_eq!(layer.num_categories(), 1);
        assert_eq!(layer.summaries[0].name, "big");
        assert_eq!(layer.summaries[0].member_count, 10);
    }

    #[test]
    fn infer_relation_historical() {
        assert_eq!(
            infer_relation("History of technology", "history", "engineering"),
            Some(RelationType::HistoricalLink)
        );
    }

    #[test]
    fn infer_relation_method() {
        assert_eq!(
            infer_relation("Spectroscopy", "chemistry", "physics"),
            Some(RelationType::SharedMethod)
        );
    }

    #[test]
    fn infer_relation_none() {
        assert_eq!(infer_relation("Quantum computing", "physics", "cs"), None);
    }

    #[test]
    fn relation_type_name_stable() {
        assert_eq!(RelationType::Studies.name(), "studies");
        assert_eq!(RelationType::HistoricalLink.name(), "historical_link");
    }
}