moyo 0.16.0

use std::collections::HashMap;

use nalgebra::Vector3;

use crate::base::{
    Linear, MoyoError, Operation, Operations, Rotation, Transformation, lattice_points,
};
use crate::subgroup::finite_group::FiniteGroup;
use crate::subgroup::klassengleiche_subgroup::KlassengleicheSubgroup;

/// Finite affine quotient `G/L` used to enumerate Klassengleiche subgroups.
///
/// `G` is the parent space group generated by the input primitive operations
/// and integer translations, and `L` is the point-group-invariant translation
/// sublattice generated by the columns of `transformation`. The quotient has
/// order `|G/T| * [T:L]`, where `T` is the parent translation subgroup.
///
/// `operations` and `parent_operations` are aligned representatives of the
/// same quotient elements. The former are expressed in the `L` basis with
/// translations reduced modulo integers; the latter are expressed in the
/// parent primitive basis. `point_indices` maps each quotient element to its
/// input point-operation representative.
///
/// A complement to `T/L` in this quotient has an inverse image in `G` whose
/// translation subgroup is exactly `L`, and therefore determines a
/// Klassengleiche subgroup with the prescribed translation lattice.
pub(super) struct AffineQuotient {
    pub(super) group: FiniteGroup,
    operations: Operations,
    pub(super) parent_operations: Operations,
    point_indices: Vec<usize>,
    transformation: Linear,
    sublattice_index: usize,
    point_order: usize,
}

impl AffineQuotient {
    pub(super) fn new(
        prim_operations: &Operations,
        transformation: &Linear,
        epsilon: f64,
    ) -> Result<Self, MoyoError> {
        let determinant = exact_determinant(transformation);
        if determinant <= 0 || determinant > i32::MAX as i128 {
            return Err(MoyoError::InvalidSublatticeTransformationError);
        }
        let sublattice_index = determinant as usize;
        let lattice_points = lattice_points(transformation);
        if lattice_points.len() != sublattice_index {
            return Err(MoyoError::InvalidSublatticeTransformationError);
        }

        let to_sublattice = Transformation::from_linear(*transformation);
        let quotient_order = prim_operations.len() * sublattice_index;
        let mut operations = Vec::with_capacity(quotient_order);
        let mut parent_operations = Vec::with_capacity(quotient_order);
        let mut point_indices = Vec::with_capacity(quotient_order);
        for (point_index, operation) in prim_operations.iter().enumerate() {
            for lattice_point in &lattice_points {
                let parent_operation = Operation::new(
                    operation.rotation,
                    operation.translation + lattice_point.map(|element| element as f64),
                );
                let mut transformed = to_sublattice
                    .transform_operation(&parent_operation)
                    .ok_or(MoyoError::InvalidSublatticeTransformationError)?;
                transformed.translation = transformed
                    .translation
                    .map(|element| reduce_mod_one(element, epsilon));
                operations.push(transformed);
                parent_operations.push(parent_operation);
                point_indices.push(point_index);
            }
        }

        let table = affine_cayley_table(&operations, epsilon)
            .ok_or(MoyoError::InvalidSublatticeTransformationError)?;
        let group = FiniteGroup::from_table(table)
            .ok_or(MoyoError::InvalidSublatticeTransformationError)?;

        Ok(Self {
            group,
            operations,
            parent_operations,
            point_indices,
            transformation: *transformation,
            sublattice_index,
            point_order: prim_operations.len(),
        })
    }

    pub(super) fn is_complement(&self, subgroup: &[usize]) -> bool {
        if subgroup.len() != self.point_order {
            return false;
        }

        let mut covered = vec![false; self.point_order];
        for &index in subgroup {
            let point_index = self.point_indices[index];
            if covered[point_index] {
                return false;
            }
            covered[point_index] = true;
        }
        covered.into_iter().all(|value| value)
    }

    pub(super) fn make_subgroup(&self, indices: &[usize]) -> KlassengleicheSubgroup {
        KlassengleicheSubgroup {
            operations: indices
                .iter()
                .map(|&index| self.operations[index].clone())
                .collect(),
            parent_operations: indices
                .iter()
                .map(|&index| self.parent_operations[index].clone())
                .collect(),
            operation_indices: indices
                .iter()
                .map(|&index| self.point_indices[index])
                .collect(),
            transformation: self.transformation,
            klassengleiche_index: self.sublattice_index,
        }
    }
}

fn exact_determinant(matrix: &Linear) -> i128 {
    let value = |row, column| matrix[(row, column)] as i128;
    value(0, 0) * (value(1, 1) * value(2, 2) - value(1, 2) * value(2, 1))
        - value(0, 1) * (value(1, 0) * value(2, 2) - value(1, 2) * value(2, 0))
        + value(0, 2) * (value(1, 0) * value(2, 1) - value(1, 1) * value(2, 0))
}

fn affine_cayley_table(operations: &Operations, epsilon: f64) -> Option<Vec<Vec<usize>>> {
    let mut rotation_indices = HashMap::<Rotation, Vec<usize>>::new();
    for (index, operation) in operations.iter().enumerate() {
        let candidates = rotation_indices.entry(operation.rotation).or_default();
        if candidates.iter().any(|&candidate| {
            translations_equivalent(
                &operation.translation,
                &operations[candidate].translation,
                epsilon,
            )
        }) {
            return None;
        }
        candidates.push(index);
    }

    let mut table = vec![vec![0; operations.len()]; operations.len()];
    for (lhs_index, lhs) in operations.iter().enumerate() {
        for (rhs_index, rhs) in operations.iter().enumerate() {
            let product = lhs.clone() * rhs.clone();
            let product_index =
                rotation_indices
                    .get(&product.rotation)?
                    .iter()
                    .find(|&&candidate| {
                        translations_equivalent(
                            &product.translation,
                            &operations[candidate].translation,
                            epsilon,
                        )
                    })?;
            table[lhs_index][rhs_index] = *product_index;
        }
    }
    Some(table)
}

pub(super) fn translations_equivalent(
    lhs: &Vector3<f64>,
    rhs: &Vector3<f64>,
    epsilon: f64,
) -> bool {
    (lhs - rhs)
        .iter()
        .all(|&value| (value - value.round()).abs() < epsilon)
}

fn reduce_mod_one(value: f64, epsilon: f64) -> f64 {
    let reduced = value.rem_euclid(1.0);
    if reduced < epsilon || 1.0 - reduced < epsilon {
        0.0
    } else {
        reduced
    }
}