//! Tissue domain definition and sampling

use bioheat_shared::{BoundingBox3D, Point3D};
use rand::Rng;
use serde::{Deserialize, Serialize};

/// 3D tissue domain for simulation
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct TissueDomain {
    /// Domain bounds
    pub bounds: BoundingBox3D,
    /// Time range [t_start, t_end] in seconds
    pub time_range: (f32, f32),
}

impl TissueDomain {
    /// Create a new tissue domain
    #[must_use]
    pub fn new(bounds: BoundingBox3D, time_range: (f32, f32)) -> Self {
        Self { bounds, time_range }
    }

    /// Create a centered cubic domain
    #[must_use]
    pub fn centered_cube(half_size: f32, t_end: f32) -> Self {
        Self {
            bounds: BoundingBox3D::centered_cube(Point3D::origin(), half_size),
            time_range: (0.0, t_end),
        }
    }

    /// Sample a random point inside the domain (spatial only)
    #[must_use]
    pub fn sample_interior<R: Rng>(&self, rng: &mut R) -> Point3D {
        Point3D::new(
            rng.gen_range(self.bounds.min.x..self.bounds.max.x),
            rng.gen_range(self.bounds.min.y..self.bounds.max.y),
            rng.gen_range(self.bounds.min.z..self.bounds.max.z),
        )
    }

    /// Sample a random point inside the domain with time
    #[must_use]
    pub fn sample_interior_4d<R: Rng>(&self, rng: &mut R) -> (Point3D, f32) {
        let point = self.sample_interior(rng);
        let t = rng.gen_range(self.time_range.0..self.time_range.1);
        (point, t)
    }

    /// Sample multiple interior points with time
    #[must_use]
    pub fn sample_interior_batch<R: Rng>(&self, rng: &mut R, n: usize) -> Vec<(Point3D, f32)> {
        (0..n).map(|_| self.sample_interior_4d(rng)).collect()
    }

    /// Sample a point on the boundary (one of the 6 faces)
    #[must_use]
    pub fn sample_boundary<R: Rng>(&self, rng: &mut R) -> Point3D {
        // Choose a random face (0-5)
        let face = rng.gen_range(0..6);

        match face {
            0 => Point3D::new(
                self.bounds.min.x, // x = x_min face
                rng.gen_range(self.bounds.min.y..self.bounds.max.y),
                rng.gen_range(self.bounds.min.z..self.bounds.max.z),
            ),
            1 => Point3D::new(
                self.bounds.max.x, // x = x_max face
                rng.gen_range(self.bounds.min.y..self.bounds.max.y),
                rng.gen_range(self.bounds.min.z..self.bounds.max.z),
            ),
            2 => Point3D::new(
                rng.gen_range(self.bounds.min.x..self.bounds.max.x),
                self.bounds.min.y, // y = y_min face
                rng.gen_range(self.bounds.min.z..self.bounds.max.z),
            ),
            3 => Point3D::new(
                rng.gen_range(self.bounds.min.x..self.bounds.max.x),
                self.bounds.max.y, // y = y_max face
                rng.gen_range(self.bounds.min.z..self.bounds.max.z),
            ),
            4 => Point3D::new(
                rng.gen_range(self.bounds.min.x..self.bounds.max.x),
                rng.gen_range(self.bounds.min.y..self.bounds.max.y),
                self.bounds.min.z, // z = z_min face
            ),
            _ => Point3D::new(
                rng.gen_range(self.bounds.min.x..self.bounds.max.x),
                rng.gen_range(self.bounds.min.y..self.bounds.max.y),
                self.bounds.max.z, // z = z_max face
            ),
        }
    }

    /// Sample boundary point with time
    #[must_use]
    pub fn sample_boundary_4d<R: Rng>(&self, rng: &mut R) -> (Point3D, f32) {
        let point = self.sample_boundary(rng);
        let t = rng.gen_range(self.time_range.0..self.time_range.1);
        (point, t)
    }

    /// Sample multiple boundary points with time
    #[must_use]
    pub fn sample_boundary_batch<R: Rng>(&self, rng: &mut R, n: usize) -> Vec<(Point3D, f32)> {
        (0..n).map(|_| self.sample_boundary_4d(rng)).collect()
    }

    /// Sample initial condition points (t = 0)
    #[must_use]
    pub fn sample_initial_batch<R: Rng>(&self, rng: &mut R, n: usize) -> Vec<(Point3D, f32)> {
        (0..n)
            .map(|_| {
                let point = self.sample_interior(rng);
                (point, self.time_range.0)
            })
            .collect()
    }

    /// Generate a regular 3D grid of points
    #[must_use]
    pub fn regular_grid(&self, resolution: (usize, usize, usize)) -> Vec<Point3D> {
        let (nx, ny, nz) = resolution;
        let size = self.bounds.size();

        let dx = size.x / (nx - 1).max(1) as f32;
        let dy = size.y / (ny - 1).max(1) as f32;
        let dz = size.z / (nz - 1).max(1) as f32;

        let mut points = Vec::with_capacity(nx * ny * nz);
        for k in 0..nz {
            for j in 0..ny {
                for i in 0..nx {
                    points.push(Point3D::new(
                        self.bounds.min.x + i as f32 * dx,
                        self.bounds.min.y + j as f32 * dy,
                        self.bounds.min.z + k as f32 * dz,
                    ));
                }
            }
        }
        points
    }

    /// Generate a regular grid with time for a specific time value
    #[must_use]
    pub fn regular_grid_at_time(
        &self,
        resolution: (usize, usize, usize),
        t: f32,
    ) -> Vec<(Point3D, f32)> {
        self.regular_grid(resolution)
            .into_iter()
            .map(|p| (p, t))
            .collect()
    }

    /// Check if a point is inside the domain
    #[must_use]
    pub fn contains(&self, point: &Point3D) -> bool {
        self.bounds.contains(point)
    }

    /// Check if a point and time are inside the domain
    #[must_use]
    pub fn contains_4d(&self, point: &Point3D, t: f32) -> bool {
        self.bounds.contains(point) && t >= self.time_range.0 && t <= self.time_range.1
    }

    /// Get the domain volume in m³
    #[must_use]
    pub fn volume(&self) -> f32 {
        self.bounds.volume()
    }

    /// Get the domain volume in mL (cm³)
    #[must_use]
    pub fn volume_ml(&self) -> f32 {
        self.volume() * 1e6
    }
}

impl Default for TissueDomain {
    fn default() -> Self {
        // 10cm x 10cm x 10cm cube, 5 minute simulation
        Self::centered_cube(0.05, 300.0)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_domain_creation() {
        let domain = TissueDomain::centered_cube(0.05, 300.0);
        assert!(domain.volume() > 0.0);
        assert!((domain.time_range.1 - 300.0).abs() < 1e-6);
    }

    #[test]
    fn test_interior_sampling() {
        let domain = TissueDomain::default();
        let mut rng = rand::thread_rng();

        for _ in 0..100 {
            let point = domain.sample_interior(&mut rng);
            assert!(domain.contains(&point));
        }
    }

    #[test]
    fn test_boundary_sampling() {
        let domain = TissueDomain::centered_cube(0.05, 300.0);
        let mut rng = rand::thread_rng();

        for _ in 0..100 {
            let point = domain.sample_boundary(&mut rng);
            // Point should be on one of the faces
            let on_x_face = (point.x - domain.bounds.min.x).abs() < 1e-6
                || (point.x - domain.bounds.max.x).abs() < 1e-6;
            let on_y_face = (point.y - domain.bounds.min.y).abs() < 1e-6
                || (point.y - domain.bounds.max.y).abs() < 1e-6;
            let on_z_face = (point.z - domain.bounds.min.z).abs() < 1e-6
                || (point.z - domain.bounds.max.z).abs() < 1e-6;
            assert!(on_x_face || on_y_face || on_z_face);
        }
    }

    #[test]
    fn test_regular_grid() {
        let domain = TissueDomain::default();
        let grid = domain.regular_grid((10, 10, 10));
        assert_eq!(grid.len(), 1000);

        // All points should be in domain
        for point in &grid {
            assert!(domain.contains(point));
        }
    }

    #[test]
    fn test_initial_sampling() {
        let domain = TissueDomain::default();
        let mut rng = rand::thread_rng();

        let samples = domain.sample_initial_batch(&mut rng, 100);
        for (_, t) in &samples {
            assert!((*t - domain.time_range.0).abs() < 1e-6);
        }
    }

    #[test]
    fn test_volume_ml() {
        let domain = TissueDomain::centered_cube(0.05, 300.0); // 10cm cube = 1000 mL
        assert!((domain.volume_ml() - 1000.0).abs() < 1.0);
    }
}