glossia 0.2.0

/// M x M grid geometry in the normal plane.
///
/// Maps sequence positions to displacement vectors for encoding
/// repeated palette colors at distinct locations.

/// CIELAB Just-Noticeable Difference — default constellation grid spacing.
pub const EPSILON: f64 = 2.3;

/// M x M grid of sequence positions in the normal plane.
#[derive(Debug, Clone)]
pub struct Constellation {
    pub m: usize,
    pub epsilon: f64,
    pub capacity: usize,
}

impl Constellation {
    pub fn new(m: usize, epsilon: f64) -> Self {
        Constellation { m, epsilon, capacity: m * m }
    }

    /// Create constellation from tube radius and step size.
    ///
    /// Uses inscribed square formula: the M×M grid corners extend to
    /// sqrt(2) × axis_distance, so we use sqrt(2)*r/ε to keep all
    /// grid points within the tube radius.
    pub fn from_radius(radius: f64, epsilon: f64) -> Self {
        let m = ((2.0_f64.sqrt() * radius / epsilon) as usize + 1).max(1);
        Self::new(m, epsilon)
    }

    /// Map sequence position to grid coordinates (a, b).
    pub fn position_to_grid(&self, j: usize) -> (usize, usize) {
        (j / self.m, j % self.m)
    }

    /// Map grid coordinates to sequence position.
    pub fn grid_to_position(&self, a: usize, b: usize) -> usize {
        a * self.m + b
    }

    /// Map grid coordinates to (alpha1, alpha2) displacements.
    pub fn grid_to_displacement(&self, a: usize, b: usize) -> (f64, f64) {
        let center = (self.m - 1) as f64 / 2.0;
        let alpha1 = (a as f64 - center) * self.epsilon;
        let alpha2 = (b as f64 - center) * self.epsilon;
        (alpha1, alpha2)
    }

    /// Snap continuous displacements to nearest grid coordinates.
    pub fn displacement_to_grid(&self, alpha1: f64, alpha2: f64) -> (usize, usize) {
        let center = (self.m - 1) as f64 / 2.0;
        let a = (alpha1 / self.epsilon + center).round() as i64;
        let b = (alpha2 / self.epsilon + center).round() as i64;
        let a = a.clamp(0, self.m as i64 - 1) as usize;
        let b = b.clamp(0, self.m as i64 - 1) as usize;
        (a, b)
    }

    /// Map sequence position to displacement vector components.
    pub fn position_to_displacement(&self, j: usize) -> (f64, f64) {
        let (a, b) = self.position_to_grid(j);
        self.grid_to_displacement(a, b)
    }

    /// Snap displacements and recover sequence position.
    pub fn displacement_to_position(&self, alpha1: f64, alpha2: f64) -> usize {
        let (a, b) = self.displacement_to_grid(alpha1, alpha2);
        self.grid_to_position(a, b)
    }
}

/// Generate center-out ordering for an M×M grid.
///
/// Returns a permutation mapping logical index -> grid position.
/// Index 0 maps to the center position (most noise-robust), and
/// subsequent indices spiral outward.
pub fn center_out_order(m: usize) -> Vec<usize> {
    let n = m * m;
    let center = (m as f64 - 1.0) / 2.0;

    // Build (distance_from_center, row-major position) pairs
    let mut positions: Vec<(f64, usize)> = (0..n).map(|j| {
        let a = (j / m) as f64;
        let b = (j % m) as f64;
        let dist = (a - center).powi(2) + (b - center).powi(2);
        (dist, j)
    }).collect();

    // Sort by distance from center (ties broken by row-major order)
    positions.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap().then(a.1.cmp(&b.1)));

    // The result maps logical index i -> grid position positions[i].1
    positions.iter().map(|(_, pos)| *pos).collect()
}

/// Per-color constellations keyed by palette index.
///
/// Each palette color gets its own Constellation sized to the local
/// tube radius.
#[derive(Debug, Clone)]
pub struct ConstellationMap {
    pub constellations: Vec<Constellation>,
    pub epsilon: f64,
}

impl ConstellationMap {
    /// Build one Constellation per palette color from local tube radii.
    pub fn new(radii: &[f64], epsilon: f64) -> Self {
        let constellations: Vec<Constellation> = radii.iter()
            .map(|&r| Constellation::from_radius(r, epsilon))
            .collect();
        ConstellationMap { constellations, epsilon }
    }

    pub fn get(&self, palette_index: usize) -> &Constellation {
        &self.constellations[palette_index]
    }

    pub fn len(&self) -> usize {
        self.constellations.len()
    }

    pub fn m_min(&self) -> usize {
        self.constellations.iter().map(|c| c.m).min().unwrap_or(0)
    }

    pub fn m_max(&self) -> usize {
        self.constellations.iter().map(|c| c.m).max().unwrap_or(0)
    }

    pub fn capacity_min(&self) -> usize {
        self.constellations.iter().map(|c| c.capacity).min().unwrap_or(0)
    }

    pub fn capacity_max(&self) -> usize {
        self.constellations.iter().map(|c| c.capacity).max().unwrap_or(0)
    }

    pub fn total_capacity(&self) -> usize {
        self.constellations.iter().map(|c| c.capacity).sum()
    }
}

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

    #[test]
    fn test_constellation_roundtrip() {
        let c = Constellation::new(5, EPSILON);
        for j in 0..c.capacity {
            let (alpha1, alpha2) = c.position_to_displacement(j);
            let recovered = c.displacement_to_position(alpha1, alpha2);
            assert_eq!(j, recovered, "Position {} should round-trip", j);
        }
    }

    #[test]
    fn test_constellation_from_radius() {
        let c = Constellation::from_radius(10.0, EPSILON);
        // M = floor(sqrt(2)*10/2.3) + 1 = floor(6.148) + 1 = 6 + 1 = 7
        assert_eq!(c.m, 7, "M should be 7 for radius=10, epsilon=2.3 (inscribed square)");
        assert_eq!(c.capacity, 49);
    }

    #[test]
    fn test_constellation_center_displacement_is_zero() {
        let c = Constellation::new(5, EPSILON);
        // Center grid position for M=5 is (2,2), which is position 12
        let (alpha1, alpha2) = c.grid_to_displacement(2, 2);
        assert!(alpha1.abs() < 1e-10 && alpha2.abs() < 1e-10,
            "Center should have zero displacement");
    }

    #[test]
    fn test_center_out_order() {
        // For a 3x3 grid, center is position 4 (row 1, col 1)
        let order = center_out_order(3);
        assert_eq!(order.len(), 9);
        // First element should be the center position
        assert_eq!(order[0], 4, "Center-out should start at grid center (1,1) = position 4");
        // All positions should be unique
        let mut sorted = order.clone();
        sorted.sort();
        assert_eq!(sorted, (0..9).collect::<Vec<_>>());
    }

    #[test]
    fn test_center_out_order_m1() {
        let order = center_out_order(1);
        assert_eq!(order, vec![0]);
    }

    #[test]
    fn test_center_out_roundtrip() {
        // Verify that encode(idx) -> grid_pos -> decode -> idx is identity
        for m in [3, 5, 8] {
            let order = center_out_order(m);
            for (idx, &grid_pos) in order.iter().enumerate() {
                let recovered = order.iter().position(|&p| p == grid_pos).unwrap();
                assert_eq!(recovered, idx, "Center-out roundtrip failed for M={}, idx={}", m, idx);
            }
        }
    }

    #[test]
    fn test_constellation_map_basics() {
        let radii = vec![10.0, 15.0, 20.0, 12.0];
        let cmap = ConstellationMap::new(&radii, EPSILON);
        assert_eq!(cmap.len(), 4);
        assert!(cmap.m_min() > 0);
        assert!(cmap.m_max() >= cmap.m_min());
    }
}