use crate::planar_map::PlanarMap;
use log::debug;
use rand::Rng;

/// Kreweras steps using the same convention as in
/// Bernardi's paper (arxiv:math/0605320)
#[derive(Debug, Clone, Copy)]
pub enum KrewerasStep {
    Left,    // a, in paper
    Down,    // b, in paper
    RightUp, // c, in paper
}

/// Produce a uniform random Kreweras exercusion
///
/// This means a path with [KrewerasStep]s, `Left=(-1, 0)`, `Down=(0, -1)`, `RightUp=(1, 1)`,
/// starting at the origin (0, 0) and ending on the diagonal x+y=0, whilst staying in x+y >= 0
pub fn random_kreweras_execursion<R: Rng>(length: usize, rng: &mut R) -> Vec<KrewerasStep> {
    let mut revwalk = Vec::with_capacity(length);
    let mut h: usize = 0;
    // Effectively performing the walk in reverse, but this gives a proper Kreweras excursion
    // The number of Kreweras walks with steps (1,1), (-1,0), (0,-1) of length
    // n from (x,y) with x+y=h to the diagonal x+y=0 is 2^{(2n+h)/3} (h+1)/(n+1)\binom{n+1}{(n-h)/3},
    // so the probability of a (1,1) step is (h+3)*(n-h)) / ((h+1)*n*3).
    for n in (1..=length).rev() {
        let q1 = (h + 3) as f64;
        let q2 = (n - h) as f64;
        let q3 = (h + 1) as f64;
        let q4 = (3 * n) as f64;
        if rng.random_bool(q1 / q3 * q2 / q4) {
            revwalk.push(KrewerasStep::RightUp);
            h += 2;
        } else {
            let step = if rng.random_bool(0.5) {
                KrewerasStep::Down
            } else {
                KrewerasStep::Left
            };
            revwalk.push(step);
            h -= 1;
        }
    }
    revwalk
}

/// Map a kreweras excursion of length 3n to a loopless triangulation of 2n triangles
///
/// The triangulation is given as a planar map of 6n half-edges.
/// Note: the given `path` must be a valid Kreweras excursion of the PlanarMap will not be
/// a valid triangulation and the function may possibly panic.
fn kreweras_excursion_to_triangulation<I>(path: I) -> PlanarMap
where
    I: ExactSizeIterator + DoubleEndedIterator + Iterator<Item = KrewerasStep>,
{
    // In the worst case scenario the planar map has 4*2n+2 half-edges = (8/3)*3n+2
    let (mut map, root) = PlanarMap::with_capacity(2 + 8 * path.len() / 3);
    let mut current_edge = map.get_next(root);
    let mut free_edges = vec![0; 2 + 8 * path.len() / 3];
    let mut age = 0;
    free_edges[root] = age;
    age += 1;
    for step in path.rev() {
        match step {
            KrewerasStep::Left => {
                let tri_edge = map.add_triangle(current_edge);
                free_edges[tri_edge] = age;
                age += 1;
                current_edge = map.get_next(tri_edge);
            }
            KrewerasStep::Down => {
                let tri_edge = map.add_triangle(current_edge);
                current_edge = tri_edge;
                free_edges[map.get_next(tri_edge)] = age;
                age += 1;
            }
            KrewerasStep::RightUp => {
                let prev_edge = map.get_previous(current_edge);
                let next_edge = map.get_next(current_edge);
                // Check if `prev_edge` is older than `next_edge` (they cannot be equal by construction)
                let prev_older = free_edges[prev_edge] < free_edges[next_edge];
                let (old_edge, new_edge) = if prev_older {
                    (prev_edge, next_edge)
                } else {
                    (next_edge, prev_edge)
                };
                map.glue_edges(current_edge, old_edge);
                current_edge = new_edge;
            }
        }
    }
    // Merge last free edges
    map.glue_edges(current_edge, map.get_next(current_edge));

    map
}

/// Generate a random triangulation with `size` triangles.
///
/// The triangulation is loopless, i.e. the dual graph contains no loops, and had spherical
/// topology, as it is given by a planar map. However, the triangulation will typically contain
/// "self-energies", i.e. two triangles that are connected along more than edge, giving two-cycles
/// in the dual graph. This means this PlanarMap will typically not be a valid non-degenerate
/// simplicial complex.
///
/// The generation works by using a bijection with random Kreweras excursions.
///
/// # Panics
/// The number of triangles must be even to be able to construct a valid triangulation,
/// as restricted by the Euler characteristic. The function will panic if an odd `size` is given.
///
/// # Examples
/// ```rust
/// # use dyntri_edt2d::generator::random_uniform_triangulation;
/// # use rand::SeedableRng;
/// # use rand_xoshiro::Xoshiro256StarStar;
/// let mut rng = Xoshiro256StarStar::seed_from_u64(42);
/// // Generate an 2D EDT triangulation using direct sampling
/// let triangulation = random_uniform_triangulation(500, &mut rng);
/// ```
pub fn random_uniform_triangulation<R: Rng>(size: usize, rng: &mut R) -> PlanarMap {
    assert!(size % 2 == 0, "The number of triangles is not even.");
    let length = 3 * (size / 2);
    let path = random_kreweras_execursion(length, rng);
    assert!(path.len() == length);
    debug!("Finished Kreweras excursion path");
    kreweras_excursion_to_triangulation(path.into_iter())
}

#[cfg(test)]
mod tests {
    use crate::planar_map::FixedPlanarMap;
    use rand::{SeedableRng, rngs::StdRng};

    use super::*;

    #[derive(Debug, Clone, Copy)]
    struct Point {
        x: i32,
        y: i32,
    }

    impl Point {
        fn new(x: i32, y: i32) -> Self {
            Self { x, y }
        }

        fn origin() -> Self {
            Self::new(0, 0)
        }

        fn step(&mut self, step: KrewerasStep) {
            match step {
                KrewerasStep::Left => self.x -= 1,
                KrewerasStep::Down => self.y -= 1,
                KrewerasStep::RightUp => {
                    self.x += 1;
                    self.y += 1;
                }
            }
        }
    }

    #[test]
    fn test_kreweras_excursion() {
        const LENGTH: usize = 3 * 24;
        let mut rng = StdRng::seed_from_u64(42);
        let path = random_kreweras_execursion(LENGTH, &mut rng);
        println!("{path:?}");

        let mut point = Point::origin();
        for step in path {
            assert!(
                point.x + point.y >= 0,
                "Excusion in forbidden part of the plane"
            );
            point.step(step);
        }
        println!("{point:?}");
        assert!(
            point.x + point.y == 0,
            "Excursion did not end on the correct line"
        );
    }

    #[test]
    fn test_kreweras_to_triangulation() {
        const LENGTH: usize = 3 * 24;
        let mut rng = StdRng::seed_from_u64(42);
        let path = random_kreweras_execursion(LENGTH, &mut rng);
        println!("{path:?}");

        let triangulation = kreweras_excursion_to_triangulation(path.into_iter());
        println!("{}", FixedPlanarMap::from(triangulation));
    }

    #[test]
    fn test_triangulation_graphs() {
        const NTRIANGLES: usize = 100_000;
        let mut rng = StdRng::seed_from_u64(42);
        let triangulation = random_uniform_triangulation(NTRIANGLES, &mut rng);

        let graph_dual = triangulation.construct_dual_graph();
        for i in 0..graph_dual.len() {
            assert_eq!(graph_dual.get_node(i).degree(), 3);
        }

        let graph_vertex = triangulation.construct_vertex_graph();
        const MAX_DEGREE: usize = 15;
        let mut counts = [0; MAX_DEGREE + 1];
        for i in 0..graph_vertex.len() {
            let degree = graph_vertex.get_node(i).degree();
            assert!(degree >= 2);
            counts[degree.min(MAX_DEGREE)] += 1;
        }
        print!("Degree: ");
        for degree in (0..counts.len()).skip(2).take(MAX_DEGREE - 2) {
            print!("{:^5} ", degree);
        }
        println!();
        print!("Counts: ");
        let total: usize = counts.iter().sum();
        for &count in counts.iter().skip(2).take(MAX_DEGREE - 2) {
            print!("{:^5.3} ", count as f64 / total as f64);
        }
        println!();
    }

    #[test]
    fn test_planarmap_to_triangulation() {
        const NTRIANGLES: usize = 10_000;
        let mut rng = StdRng::seed_from_u64(42);
        let planar_map = random_uniform_triangulation(NTRIANGLES, &mut rng);

        let triangulation = planar_map.construct_triangulation();
        assert_eq!(triangulation.num_triangles(), NTRIANGLES);
        assert_eq!(triangulation.num_edges() * 2, planar_map.len()); // Compare number of half edges
        assert_eq!(triangulation.num_vertices(), NTRIANGLES / 2 + 2); // Check Euler char for sphere

        assert!(triangulation.check_biconnectivity());
        assert!(triangulation.check_vertex_consistency());
    }
}
