dyntri-core 0.10.2

use std::hash::Hash;

use indexmap::IndexSet;

use crate::triangulation::simplices::{canonical_oriented_simplex, canonical_simplex};
use crate::triangulation::{Triangulation2D, Triangulation3D, Triangulation4D};
use narray::NArray;

/// Trait representing a `D`-simplex as an ordered list of its vertices `[usize; D+1]`
///
/// Note: this used to avoid needing to use const generics, which do not yet allow
/// expressions like `N-1`, which is needed for the [`get_face()`] function.
/// One should think of the trait as just being a `[usize; D+1]` array.
pub trait FaceSimplex: NArray<Item = usize> {
    type Boundary: NArray<Item = usize> + Copy;
    type Boundaries: NArray<Item = Self::Boundary>;
    const BOUNDARIES: Self::Boundaries;
}

/// Convenience alias representing the `(D-1)`-subsimplex `Face: [usize; D]` of the
/// `D`-simplex `S: [usize; D+1]`
pub type Face<S> = <<S as FaceSimplex>::Boundary as NArray>::MapOutput<usize>;
/// Convenience alias representing all `(D-1)`-subsimplices `AllFaces: [[usize; D], D+1]`
/// of the `D`-simplex `S: [usize; D+1]`, which has `D+1` faces.
pub type AllFaces<S> = <<S as FaceSimplex>::Boundaries as NArray>::MapOutput<Face<S>>;

/// Returns the `idx`th `face: [usize; DIM]` of the given `simplex: [usize; DIM+1]`
///
/// The `idx`th face means the sub-simplex obtained by deleting the `idx`th vertex
/// from the simplex (which is represented as an ordered set of vertices), i.e.
/// 1st face of [0, 1, 2, 3] -> [0, 2, 3].
pub fn get_face<S: FaceSimplex>(simplex: S, idx: usize) -> Face<S> {
    S::BOUNDARIES[idx].map(move |i| simplex[i])
}

/// Returns all `DIM+1` `faces: [[usize; DIM]; DIM+1]`` of the given
/// simplex: [usize; DIM+1]` in order.
///
/// See also [`get_face()`] to get just a single face, and to clarify what "in order" means.
pub fn get_faces<S: FaceSimplex>(simplex: S) -> AllFaces<S> {
    S::BOUNDARIES.map(|face| face.map(|i| simplex[i]))
}

macro_rules! impl_face_simplex {
    ($dim:literal, $boundary:expr) => {
        impl FaceSimplex for [usize; $dim + 1] {
            type Boundary = [usize; $dim];
            type Boundaries = [[usize; $dim]; $dim + 1];
            const BOUNDARIES: Self::Boundaries = $boundary;
        }
    };
}

impl_face_simplex!(1, [[1], [0]]);
impl_face_simplex!(2, [[1, 2], [0, 2], [0, 1]]);
impl_face_simplex!(3, [[1, 2, 3], [0, 2, 3], [0, 1, 3], [0, 1, 2]]);
impl_face_simplex!(
    4,
    [
        [1, 2, 3, 4],
        [0, 2, 3, 4],
        [0, 1, 3, 4],
        [0, 1, 2, 4],
        [0, 1, 2, 3],
    ]
);

type MappedFaces<S> = <AllFaces<S> as NArray>::MapOutput<Boundary>;

/// Given a collection of N-simplices, returns ( (N-1)-simplex/face list, boundary map),
/// where the [`Boundary`] labels in the boundary map refer to the (N-1)-simplices
/// at the corresponding position in the list.
///
/// Note the face list is `Vec<[usize; N]>`, where each element is the ordered vertex set of the
/// subsimplex. The boundary map is `Vec<[Boundary; N+1]>`, where each element represents the
/// boundaries in order of that simplex given by [`Boundary`] labels. See also [`get_face()`] for
/// details on what "in order" means.
pub fn compute_boundary_map<S, I>(simplices: I) -> (Vec<Face<S>>, Vec<MappedFaces<S>>)
where
    S: FaceSimplex,
    Face<S>: CanonicalSimplex + Hash + Eq,
    I: Iterator<Item = S>,
{
    // Hashset that keeps track of insertion order
    let mut facemap = IndexSet::new();
    let boundary_map = simplices
        .map(|simplex| {
            // Use parity sign in the boundary map for odd and even faces
            let mut sign = true;
            get_faces(simplex).map(|face| {
                // Get unique canonical oriented representation
                let (orientation, face) = face.make_canonical_oriented();
                // Get index (position in facemap) of face (inserting if not yet present)
                let label = facemap.insert_full(face).0;
                // Combining orientation and sign with flip such that true true => true
                let bound = Boundary::new(label, !(orientation ^ sign));
                sign = !sign;
                bound
            })
        })
        .collect();
    let faces = facemap.into_iter().collect();
    (faces, boundary_map)
}

#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub enum Parity {
    Negative,
    Positive,
}

/// [Boundary] label containing and `idx` refering to the index of the boundary face (position
/// in its simplex list) and `parity` encoding the orientation of faces as boundary.
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub struct Boundary {
    pub parity: Parity,
    pub idx: usize,
}

impl Boundary {
    pub fn new(idx: usize, parity: bool) -> Self {
        Self {
            parity: if parity {
                Parity::Positive
            } else {
                Parity::Negative
            },
            idx,
        }
    }
}

pub trait CanonicalSimplex {
    /// Orders the vertices of the simplex from small to big to give a canonical ordering for comparison
    fn make_canonical(self) -> Self;

    /// Orders the vertices of the simplex from small to big to give a canonical ordering for comparison,
    /// additionally giving the parity of the permulation to have ordering (true-even, false-odd)
    fn make_canonical_oriented(self) -> (bool, Self);
}

impl<T: Ord, const N: usize> CanonicalSimplex for [T; N] {
    fn make_canonical(self) -> Self {
        canonical_simplex(self)
    }

    fn make_canonical_oriented(self) -> (bool, [T; N]) {
        canonical_oriented_simplex(self)
    }
}

pub type BoundaryMap<const N: usize> = Vec<[Boundary; N]>;

impl Triangulation2D {
    pub fn boundary_maps(&self) -> (BoundaryMap<2>, BoundaryMap<3>) {
        let simplices2_iter = self.simplices.iter().map(|simplex| simplex.vertices);
        let (simplices1, boundmap2) = compute_boundary_map(simplices2_iter);
        let (_, boundmap1) = compute_boundary_map(simplices1.into_iter());

        (boundmap1, boundmap2)
    }
}

impl Triangulation3D {
    pub fn boundary_maps(&self) -> (BoundaryMap<2>, BoundaryMap<3>, BoundaryMap<4>) {
        let simplices3_iter = self.simplices.iter().map(|simplex| simplex.vertices);
        let (simplices2, boundmap3) = compute_boundary_map(simplices3_iter);
        let (simplices1, boundmap2) = compute_boundary_map(simplices2.into_iter());
        let (_, boundmap1) = compute_boundary_map(simplices1.into_iter());

        (boundmap1, boundmap2, boundmap3)
    }
}

impl Triangulation4D {
    pub fn boundary_maps(
        &self,
    ) -> (
        BoundaryMap<2>,
        BoundaryMap<3>,
        BoundaryMap<4>,
        BoundaryMap<5>,
    ) {
        let simplices4_iter = self.simplices.iter().map(|simplex| simplex.vertices);
        let (simplices3, boundmap4) = compute_boundary_map(simplices4_iter);
        let (simplices2, boundmap3) = compute_boundary_map(simplices3.into_iter());
        let (simplices1, boundmap2) = compute_boundary_map(simplices2.into_iter());
        let (_, boundmap1) = compute_boundary_map(simplices1.into_iter());

        (boundmap1, boundmap2, boundmap3, boundmap4)
    }
}

mod narray {
    use std::ops::Index;

    /// Trait used to generalize over all [T; N] implementing only the needed methods,
    /// indexing and mapping.
    pub trait NArray:
        Sized
        + Index<usize, Output = <Self as NArray>::Item>
        + IntoIterator<Item = <Self as NArray>::Item>
    {
        type Item;
        type MapOutput<U>: NArray<Item = U>;

        fn map<U, F: FnMut(<Self as NArray>::Item) -> U>(self, f: F) -> Self::MapOutput<U>;
    }

    impl<T, const N: usize> NArray for [T; N] {
        type Item = T;
        type MapOutput<U> = [U; N];

        fn map<U, F: FnMut(T) -> U>(self, f: F) -> [U; N] {
            self.map(f)
        }
    }
}