//! A fixed-volume implemention of a Markov-chain Monte Carlo simulation
//! for 2D Causal dynamical triangulations.
#![allow(clippy::manual_is_multiple_of)]

use dyntri_core::triangulation::CausalTriangulation2D;
use rand::Rng;
use serde::{Deserialize, Serialize};

use crate::triangulation::fixed_triangulation::FixedTriangulation;

/// The ratio between the flip move and the relocation move
pub const DEFAULT_MOVE_RATIO: f64 = 0.5;

/// The main structure storing the triangulation and being in charge of performing the Monte Carlo
/// moves.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Simulation {
    /// The fixed-volume 2D CDT triangulation
    triangulation: FixedTriangulation,
    move_ratio: f64, // The selection probability between flip and relocation move
}

impl Simulation {
    /// Create a new fixed size 2D CDT simulation with `ntriangles` triangles
    /// and `ntime` timeslices.
    ///
    /// Note: that the triangulations with the flat initial configuration can only have
    /// `ntriangles` given by: `2 * ntime * average_slice_length`,
    /// where the `average_slice_length` is an integer.
    pub fn new(ntriangles: usize, ntime: u16, move_ratio: f64) -> Self {
        let ntime = ntime as usize;
        assert!(
            ntriangles % (2 * ntime) == 0,
            "`ntriangles` must be a multiple of `2*ntime`."
        );
        assert!(
            ntriangles >= 2 * ntime,
            "Cannot create a universe with less than two triangles per timeslice!"
        );
        let start_slice_length = ntriangles / ntime / 2; // In terms of vertices per slice
        Simulation {
            triangulation: FixedTriangulation::new(start_slice_length, ntime as u16),
            move_ratio,
        }
    }

    /// Create a new `Universe` with default parameters
    pub fn default(ntriangles: usize, ntime: u16) -> Self {
        Self::new(ntriangles, ntime, DEFAULT_MOVE_RATIO)
    }
}

/// Implementation on performing the Markov-chain moves
impl Simulation {
    /// Perform a Markov-Chain Monte Carlo move
    pub fn step<R: Rng + ?Sized>(&mut self, rng: &mut R) {
        if rng.random_bool(self.move_ratio) {
            // Flip
            let label = self.triangulation.sample_triangle(rng);
            self.triangulation.flip(label);
        } else {
            // Relocate
            let Some(from) = self.triangulation.sample_order_four(rng) else {
                // If there are no order-four vertices simply reject.
                return;
            };
            let to = self.triangulation.sample_triangle(rng);
            self.triangulation.relocate(from, to);
        }
    }
}

/// Implementation of measurement on the universe
impl Simulation {
    /// Computes the current size of the triangulation
    pub fn ntriangles(&self) -> usize {
        self.triangulation.triangles.len()
    }

    /// Returns the number of timeslices
    pub fn ntime(&self) -> u16 {
        self.triangulation.ntime
    }

    pub fn construct_triangulation(&self) -> CausalTriangulation2D {
        self.triangulation.reconstruct_causal_triangulation2d()
    }
}

#[cfg(test)]
mod tests {
    use rand::SeedableRng;
    use rand_xoshiro::Xoshiro128StarStar;

    use super::*;

    #[test]
    fn test_validity_triangulation() {
        let mut rng = Xoshiro128StarStar::seed_from_u64(42);
        let mut simulation = Simulation::default(200, 10);

        for _ in 0..10_000 {
            simulation.step(&mut rng);
        }

        simulation.triangulation.assert_triangulation();
    }
}
