dyntri-core 0.10.2

use serde::{Deserialize, Serialize};
use serde_with::serde_as;

/// Euclidean (N-1)-simplex, built from `N` vertices.
///
/// The intended `neighbours` and `vertices` ordering is such that the vertex corresponding
/// to a neighbouring simplex is the vertex opposite the simplex. Additionally the overall
/// ordering is supposed to be consistent, such that the simplices are *oriented*.
#[serde_as]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
pub struct Simplex<const N: usize> {
    #[serde_as(as = "[_; N]")]
    pub neighbours: [usize; N],
    #[serde_as(as = "[_; N]")]
    pub vertices: [usize; N],
}

impl<const N: usize> Default for Simplex<N> {
    fn default() -> Self {
        Self {
            neighbours: [0; N],
            vertices: [0; N],
        }
    }
}

pub type Triangle = Simplex<3>;
pub type Simplex2 = Simplex<3>;
pub type Tetrahedron = Simplex<4>;
pub type Simplex3 = Simplex<4>;
pub type Simplex4 = Simplex<5>;

/// Methods for working with [`Simplex<N>`] when `neighbours` are encoded by `half-faces`,
/// see [`Triangulation`](super::triangulation::Triangulation)
impl<const N: usize> Simplex<N> {
    /// Returns the label of the neighbouring (N-1)-simplex on the `index` side.
    ///
    /// # Panic
    /// This method will panic `index >= N`
    #[inline]
    pub fn get_neighbour(&self, index: usize) -> usize {
        let nbr = self.neighbours[index];
        nbr / N
    }

    /// Returns the label and backindex of neighbouring (N-1)-simplex on the `index` side.
    ///
    /// # Panic
    /// This method will panic `index >= N`
    #[inline]
    pub fn get_neighbour_backindex(&self, index: usize) -> (usize, usize) {
        let nbr = self.neighbours[index];
        (nbr / N, nbr % N)
    }

    /// Returns the neighbouring (N-1)-simplex labels of the [`Simplex`]
    #[inline]
    pub fn get_neighbours(&self) -> [usize; N] {
        let adjacent_edges = self.neighbours;
        adjacent_edges.map(|nbr| nbr / N)
    }

    /// Returns the neighbouring (N-1)-simplex labels and back indices in tuple pairs
    #[inline]
    pub fn get_neighbours_backindices(&self) -> [(usize, usize); N] {
        let adjacent_edges = self.neighbours;
        adjacent_edges.map(|nbr| (nbr / N, nbr % N))
    }
}

/// Lorentian/causal (N-1)-simplex, built from N vertices.
///
/// The intended `neighbours` and `vertices` ordering is such that the vertex corresponding
/// to a neighbouring simplex is the vertex opposite the simplex. Additionally the overall
/// ordering is supposed to be consistent, such that the simplices are *oriented*.
#[serde_as]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
pub struct CausalSimplex<K, const N: usize> {
    pub t: u16,
    pub kind: K,
    #[serde_as(as = "[_; N]")]
    pub neighbours: [usize; N],
    #[serde_as(as = "[_; N]")]
    pub vertices: [usize; N],
}

impl<K: Default, const N: usize> Default for CausalSimplex<K, N> {
    fn default() -> Self {
        Self {
            t: u16::default(),
            kind: K::default(),
            neighbours: [0; N],
            vertices: [0; N],
        }
    }
}

/// Methods for working with [`CausalSimplex<N>`] when `neighbours` are encoded
/// by `half-faces`, see [`Triangulation`](super::triangulation::Triangulation)
impl<K, const N: usize> CausalSimplex<K, N> {
    /// Returns the label of the neighbouring (N-1)-simplex on the `index` side.
    ///
    /// # Panic
    /// This method will panic `index >= N`
    #[inline]
    pub fn get_neighbour(&self, index: usize) -> usize {
        let nbr = self.neighbours[index];
        nbr / N
    }

    /// Returns the label and backindex of neighbouring (N-1)-simplex on the `index` side.
    ///
    /// # Panic
    /// This method will panic `index >= N`
    #[inline]
    pub fn get_neighbour_backindex(&self, index: usize) -> (usize, usize) {
        let nbr = self.neighbours[index];
        (nbr / N, nbr % N)
    }

    /// Returns the neighbouring (N-1)-simplex labels of the [`Simplex`]
    #[inline]
    pub fn get_neighbours(&self) -> [usize; N] {
        let adjacent_edges = self.neighbours;
        adjacent_edges.map(|nbr| nbr / N)
    }

    /// Returns the neighbouring (N-1)-simplex labels and back indices in tuple pairs
    #[inline]
    pub fn get_neighbours_backindices(&self) -> [(usize, usize); N] {
        let adjacent_edges = self.neighbours;
        adjacent_edges.map(|nbr| (nbr / N, nbr % N))
    }
}

/// Kind/type/variant specification of a triangle/2-simplex.
///
/// Variant `S21` refers to a triangle with 2 vertices at the previous time-slice
/// to the other one, where variant `S12` is the opposite.
#[derive(
    Debug, Clone, Copy, PartialEq, Eq, Hash, Default, PartialOrd, Ord, Serialize, Deserialize,
)]
pub enum CausalSimplex2Kind {
    #[default]
    S21 = 0,
    S12 = 1,
}
pub type CausalTriangle = CausalSimplex<CausalSimplex2Kind, 3>;
pub type CausalSimplex2 = CausalSimplex<CausalSimplex2Kind, 3>;

/// Kind/type/variant specification of a tetrahedron/3-simplex.
///
/// Variant `S31` refers to a tetrahedron with 3 vertices at the previous time-slice
/// to 1 vertex in the next, and variant `S13` is the opposite.
/// Variant `S22` refers to a tetrahedron with 2 vertices at both time-slices.
#[derive(
    Debug, Clone, Copy, PartialEq, Eq, Hash, Default, PartialOrd, Ord, Serialize, Deserialize,
)]
pub enum CausalSimplex3Kind {
    #[default]
    S31 = 0,
    S22 = 1,
    S13 = 2,
}
pub type CausalTetrahedron = CausalSimplex<CausalSimplex3Kind, 4>;
pub type CausalSimplex3 = CausalSimplex<CausalSimplex3Kind, 4>;

/// Kind/type/variant specification of a 4-simplex.
///
/// Variants `S41` and `S14` refer to a 4-simplex with 4 vertices at the previous time-slice
/// to 1 vertex in the next and the reverse respectivally, i.e tertahedron to point.
/// Variants `S32` and `S23` refer to a 4-simplex with 3 vertices at the previous time-slice
/// to 2 vertices in the next and the reverse respectivally, i.e. triangle to edge.
#[derive(
    Debug, Clone, Copy, PartialEq, Eq, Hash, Default, PartialOrd, Ord, Serialize, Deserialize,
)]
pub enum CausalSimplex4Kind {
    #[default]
    S41 = 0,
    S32 = 1,
    S23 = 2,
    S14 = 3,
}
pub type CausalSimplex4 = CausalSimplex<CausalSimplex4Kind, 5>;

impl<K, const N: usize> From<CausalSimplex<K, N>> for Simplex<N> {
    /// Convert [`CausalSimplex`] to [`Simplex`] by simply dropping `time` and `kind` information.
    fn from(value: CausalSimplex<K, N>) -> Self {
        Self {
            neighbours: value.neighbours,
            vertices: value.vertices,
        }
    }
}

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default, PartialOrd, Ord)]
pub enum CausalOrientation {
    Past = 0,
    #[default]
    Space = 1,
    Future = 2,
}

pub trait CausalSimplexKind {
    /// Return spacetime orientation representing the orientation of the simplex
    ///
    /// Consider a path starting inside the simplex moving through the simplex faces:
    /// The orientation is `Past` if there is a face moving through which goes towards the past;
    /// it is `Future` if there is a face moving through which goes towards the future;
    /// and it is `Space` if moving through all faces moves in spacelike direction.
    /// Note: it is indeed only possible for CDT simplices to be only one of these 3.
    fn orientation(&self) -> CausalOrientation;
    fn signature(&self) -> (u8, u8);
}

impl CausalSimplexKind for CausalSimplex2Kind {
    fn orientation(&self) -> CausalOrientation {
        match self {
            CausalSimplex2Kind::S21 => CausalOrientation::Past,
            CausalSimplex2Kind::S12 => CausalOrientation::Future,
        }
    }

    fn signature(&self) -> (u8, u8) {
        match self {
            CausalSimplex2Kind::S21 => (2, 1),
            CausalSimplex2Kind::S12 => (1, 2),
        }
    }
}

impl CausalSimplexKind for CausalSimplex3Kind {
    fn orientation(&self) -> CausalOrientation {
        match self {
            CausalSimplex3Kind::S31 => CausalOrientation::Past,
            CausalSimplex3Kind::S22 => CausalOrientation::Space,
            CausalSimplex3Kind::S13 => CausalOrientation::Future,
        }
    }

    fn signature(&self) -> (u8, u8) {
        match self {
            CausalSimplex3Kind::S31 => (3, 1),
            CausalSimplex3Kind::S22 => (2, 2),
            CausalSimplex3Kind::S13 => (1, 3),
        }
    }
}

impl CausalSimplexKind for CausalSimplex4Kind {
    fn orientation(&self) -> CausalOrientation {
        match self {
            CausalSimplex4Kind::S41 => CausalOrientation::Past,
            CausalSimplex4Kind::S32 => CausalOrientation::Space,
            CausalSimplex4Kind::S23 => CausalOrientation::Space,
            CausalSimplex4Kind::S14 => CausalOrientation::Future,
        }
    }
    fn signature(&self) -> (u8, u8) {
        match self {
            CausalSimplex4Kind::S41 => (4, 1),
            CausalSimplex4Kind::S32 => (3, 2),
            CausalSimplex4Kind::S23 => (2, 3),
            CausalSimplex4Kind::S14 => (1, 4),
        }
    }
}

/// Orders the vertices of the simplex from small to big to give a canonical ordering for comparison
pub(crate) fn canonical_simplex<T: Ord, const N: usize>(mut simplex: [T; N]) -> [T; N] {
    simplex.sort();
    simplex
}

/// 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)
pub(crate) fn canonical_oriented_simplex<T: Ord, const N: usize>(
    mut simplex: [T; N],
) -> (bool, [T; N]) {
    // Use bubble sort
    let mut parity = true;
    loop {
        let mut unchanged = true;
        for i in 0..(N - 1) {
            if simplex[i] > simplex[i + 1] {
                parity = !parity;
                simplex.swap(i, i + 1);
                unchanged = false;
            }
        }
        if unchanged {
            break;
        }
    }
    (parity, simplex)
}

#[cfg(test)]
mod tests {
    use crate::triangulation::simplices::canonical_oriented_simplex;

    #[test]
    fn test_canonical_oriented_simplex() {
        assert_eq!(
            canonical_oriented_simplex([0, 1, 2, 3, 4]),
            (true, [0, 1, 2, 3, 4])
        );
        assert_eq!(
            canonical_oriented_simplex([1, 0, 2, 3, 4]),
            (false, [0, 1, 2, 3, 4])
        );
        assert_eq!(
            canonical_oriented_simplex([3, 1, 2, 0, 4]),
            (false, [0, 1, 2, 3, 4])
        );
        assert_eq!(
            canonical_oriented_simplex([4, 1, 2, 3, 0]),
            (false, [0, 1, 2, 3, 4])
        );
        assert_eq!(
            canonical_oriented_simplex([1, 3, 0, 4, 2]),
            (true, [0, 1, 2, 3, 4])
        );
    }
}