//! Property-based tests for geometric algorithms
//!
//! These tests verify fundamental geometric invariants and properties using
//! proptest for randomized testing across a wide range of inputs.

use approx::assert_relative_eq;
use proptest::prelude::*;
use scirs2_core::ndarray::{array, Array1, Array2};
use scirs2_spatial::*;

// Strategy for generating valid 2D points
fn point2d_strategy() -> impl Strategy<Value = (f64, f64)> {
    ((-100.0..100.0), (-100.0..100.0))
}

// Strategy for generating arrays of 2D points
fn points2d_array_strategy(n: usize) -> impl Strategy<Value = Array2<f64>> {
    prop::collection::vec(point2d_strategy(), n..=n).prop_map(move |points| {
        let mut arr = Array2::zeros((n, 2));
        for (i, (x, y)) in points.iter().enumerate() {
            arr[[i, 0]] = *x;
            arr[[i, 1]] = *y;
        }
        arr
    })
}

// Strategy for generating non-zero vectors
fn non_zero_vector_strategy(dim: usize) -> impl Strategy<Value = Vec<f64>> {
    prop::collection::vec(-100.0..100.0, dim..=dim).prop_filter("vector must be non-zero", |v| {
        v.iter().any(|&x: &f64| x.abs() > 1e-10)
    })
}

// Test: Triangle inequality for distance metrics
proptest! {
    #[test]
    fn test_triangle_inequality_euclidean(
        p1 in non_zero_vector_strategy(10),
        p2 in non_zero_vector_strategy(10),
        p3 in non_zero_vector_strategy(10)
    ) {
        let d12 = distance::euclidean(&p1, &p2);
        let d23 = distance::euclidean(&p2, &p3);
        let d13 = distance::euclidean(&p1, &p3);

        // Triangle inequality: d(p1, p3) <= d(p1, p2) + d(p2, p3)
        prop_assert!(d13 <= d12 + d23 + 1e-10, "Triangle inequality violated");
    }

    #[test]
    fn test_triangle_inequality_manhattan(
        p1 in non_zero_vector_strategy(10),
        p2 in non_zero_vector_strategy(10),
        p3 in non_zero_vector_strategy(10)
    ) {
        let d12 = distance::manhattan(&p1, &p2);
        let d23 = distance::manhattan(&p2, &p3);
        let d13 = distance::manhattan(&p1, &p3);

        prop_assert!(d13 <= d12 + d23 + 1e-10, "Triangle inequality violated for Manhattan distance");
    }

    // Test: Distance symmetry
    #[test]
    fn test_distance_symmetry_euclidean(
        p1 in non_zero_vector_strategy(10),
        p2 in non_zero_vector_strategy(10)
    ) {
        let d12 = distance::euclidean(&p1, &p2);
        let d21 = distance::euclidean(&p2, &p1);

        assert_relative_eq!(d12, d21, epsilon = 1e-10);
    }

    // Test: Distance non-negativity
    #[test]
    fn test_distance_non_negative(
        p1 in non_zero_vector_strategy(10),
        p2 in non_zero_vector_strategy(10)
    ) {
        let d = distance::euclidean(&p1, &p2);
        prop_assert!(d >= 0.0, "Distance must be non-negative");
    }

    // Test: Identity of indiscernibles (d(x,x) = 0)
    #[test]
    fn test_distance_identity(p in non_zero_vector_strategy(10)) {
        let d = distance::euclidean(&p, &p);
        assert_relative_eq!(d, 0.0, epsilon = 1e-10);
    }
}

// Test: Convex hull properties
proptest! {
    #[test]
    fn test_convex_hull_contains_all_points(points in points2d_array_strategy(10)) {
        if let Ok(hull) = ConvexHull::new(&points.view()) {
            // All original points should be inside or on the hull
            for i in 0..points.nrows() {
                let point = array![points[[i, 0]], points[[i, 1]]];
                // Note: contains() might have tolerance issues, so we check vertices
                let vertex_indices = hull.vertex_indices();
                let is_vertex = vertex_indices.contains(&i);

                // Every point is either a vertex or inside
                // We just check this doesn't panic
                let _ = hull.contains(point.to_vec());
            }
        }
    }

    #[test]
    fn test_convex_hull_vertices_count(points in points2d_array_strategy(20)) {
        if let Ok(hull) = ConvexHull::new(&points.view()) {
            let vertices = hull.vertices();
            // Convex hull of n points in 2D has at most n vertices
            prop_assert!(vertices.len() <= points.nrows(),
                "Convex hull has more vertices than input points");
            // Must have at least 3 vertices for a proper 2D hull (unless degenerate)
            if points.nrows() >= 3 {
                prop_assert!(vertices.len() >= 3 || vertices.len() == points.nrows(),
                    "Expected at least 3 vertices for non-degenerate hull");
            }
        }
    }
}

// Test: KD-Tree properties
proptest! {
    #[test]
    fn test_kdtree_nearest_neighbor_consistency(points in points2d_array_strategy(50)) {
        if let Ok(tree) = KDTree::new(&points) {
            // Query for nearest neighbor of each point
            for i in 0..points.nrows().min(10) {  // Test first 10 points
                let query = vec![points[[i, 0]], points[[i, 1]]];
                if let Ok((indices, distances)) = tree.query(&query, 1) {
                    // Nearest neighbor should have distance 0 (the point itself)
                    if !distances.is_empty() {
                        prop_assert!(distances[0] < 1e-6,
                            "Point should be its own nearest neighbor with distance ~0");
                    }
                }
            }
        }
    }

    #[test]
    fn test_kdtree_query_returns_k_neighbors(
        points in points2d_array_strategy(100),
        k in 1usize..=10
    ) {
        if let Ok(tree) = KDTree::new(&points) {
            let query = vec![50.0, 50.0];
            if let Ok((indices, distances)) = tree.query(&query, k) {
                // Should return exactly k neighbors (or all points if k > n)
                let expected = k.min(points.nrows());
                prop_assert_eq!(indices.len(), expected,
                    "KD-Tree should return {} neighbors", expected);
                prop_assert_eq!(distances.len(), expected,
                    "Should have {} distances", expected);

                // Distances should be non-decreasing (sorted)
                for i in 1..distances.len() {
                    prop_assert!(distances[i] >= distances[i-1] - 1e-10,
                        "Distances should be sorted");
                }
            }
        }
    }
}

// Test: Distance transform properties
proptest! {
    #[test]
    fn test_distance_transform_feature_pixels_zero(
        rows in 10usize..=50,
        cols in 10usize..=50
    ) {
        let mut binary = Array2::zeros((rows, cols));

        // Add some feature pixels
        for i in (0..rows).step_by(5) {
            for j in (0..cols).step_by(5) {
                binary[[i, j]] = 1;
            }
        }

        if let Ok(distances) = euclidean_distance_transform::<f64>(
            &binary.view(),
            DistanceTransformMetric::Euclidean
        ) {
            // All feature pixels should have distance 0
            for i in 0..rows {
                for j in 0..cols {
                    if binary[[i, j]] != 0 {
                        assert_relative_eq!(distances[[i, j]], 0.0, epsilon = 1e-10);
                    } else {
                        // Background pixels should have positive distance
                        prop_assert!(distances[[i, j]] > 0.0,
                            "Background pixel should have positive distance");
                    }
                }
            }
        }
    }

    #[test]
    fn test_distance_transform_monotonicity(
        rows in 10usize..=30,
        cols in 10usize..=30
    ) {
        // Single feature in center
        let mut binary = Array2::zeros((rows, cols));
        binary[[rows/2, cols/2]] = 1;

        if let Ok(distances) = euclidean_distance_transform::<f64>(
            &binary.view(),
            DistanceTransformMetric::Euclidean
        ) {
            // Distances should increase as we move away from center
            let center_i = rows / 2;
            let center_j = cols / 2;

            // Check that distance increases along a ray
            if center_i > 0 && center_i < rows - 1 {
                let d0 = distances[[center_i, center_j]];
                let d1 = distances[[center_i - 1, center_j]];
                let d2 = distances[[center_i - 2, center_j]];

                prop_assert!(d0 < d1 && d1 <= d2 + 0.1,
                    "Distance should generally increase moving away from feature");
            }
        }
    }
}

// Test: Spatial statistics properties
proptest! {
    #[test]
    fn test_morans_i_bounds(n in 5usize..=50) {
        // Generate random values
        let values = Array1::from_shape_fn(n, |_| rand::random::<f64>() * 10.0);

        // Create simple adjacency matrix
        let mut weights = Array2::zeros((n, n));
        for i in 0..n-1 {
            weights[[i, i+1]] = 1.0;
            weights[[i+1, i]] = 1.0;
        }

        if let Ok(moran) = morans_i(&values.view(), &weights.view()) {
            // Moran's I typically ranges from -1 to 1, but can exceed these bounds slightly
            prop_assert!((-2.0..=2.0).contains(&moran),
                "Moran's I should be within reasonable bounds: {}", moran);
        }
    }

    #[test]
    fn test_gearys_c_positivity(n in 5usize..=50) {
        let values = Array1::from_shape_fn(n, |_| rand::random::<f64>() * 10.0);

        let mut weights = Array2::zeros((n, n));
        for i in 0..n-1 {
            weights[[i, i+1]] = 1.0;
            weights[[i+1, i]] = 1.0;
        }

        if let Ok(geary) = gearys_c(&values.view(), &weights.view()) {
            // Geary's C should be non-negative
            prop_assert!(geary >= 0.0,
                "Geary's C must be non-negative: {}", geary);
        }
    }
}

// Test: Variogram properties
proptest! {
    #[test]
    fn test_variogram_non_negative(n in 10usize..=50) {
        let coords = Array2::from_shape_fn((n, 2), |_| rand::random::<f64>() * 100.0);
        let values = Array1::from_shape_fn(n, |_| rand::random::<f64>() * 10.0);

        if let Ok((lags, gamma)) = experimental_variogram(
            &coords.view(),
            &values.view(),
            10,
            None
        ) {
            // All variogram values should be non-negative
            for &g in gamma.iter() {
                prop_assert!(g >= 0.0,
                    "Variogram values must be non-negative: {}", g);
            }

            // Lags should be positive
            for &lag in lags.iter() {
                prop_assert!(lag > 0.0,
                    "Lag distances must be positive: {}", lag);
            }
        }
    }

    #[test]
    fn test_variogram_monotonicity(n in 20usize..=50) {
        // Create spatially autocorrelated data: values depend on position
        // so nearby points have similar values and the variogram increases.
        let coords = Array2::from_shape_fn((n, 2), |(i, _j)| {
            // Use deterministic-ish spread based on index to get varied coords
            let base = (i as f64) * 100.0 / (n as f64);
            base + rand::random::<f64>() * 5.0
        });
        // Values are a smooth function of position (+ small noise) to ensure
        // spatial autocorrelation: nearby coords => similar values.
        let values = Array1::from_shape_fn(n, |i| {
            let x = coords[[i, 0]];
            let y = coords[[i, 1]];
            (x + y) * 0.1 + rand::random::<f64>() * 0.5
        });

        if let Ok((_lags, gamma)) = experimental_variogram(
            &coords.view(),
            &values.view(),
            10,
            None
        ) {
            // For spatially autocorrelated data, the variogram should generally
            // increase with distance. We use a very relaxed check:
            // second half mean should be non-trivially positive.
            if gamma.len() >= 3 {
                let first_half_mean = gamma.slice(s![..gamma.len()/2]).mean();
                let second_half_mean = gamma.slice(s![gamma.len()/2..]).mean();

                if let (Some(first), Some(second)) = (first_half_mean, second_half_mean) {
                    // With spatially autocorrelated data, the second half
                    // (larger lags) should have at least ~30% of the first-half
                    // variance. This is very generous to avoid flakiness.
                    prop_assert!(second >= first * 0.3,
                        "Variogram should show general increasing trend: \
                         first_half={}, second_half={}", first, second);
                }
            }
        }
    }
}

// Test: Coordinate transformation properties
proptest! {
    #[test]
    fn test_utm_roundtrip_accuracy(
        lat in -80.0..84.0,
        lon in -180.0..180.0
    ) {
        if let Ok((easting, northing, zone_num, _zone_letter)) = geographic_to_utm(lat, lon) {
            // Create UTMZone from zone number and hemisphere
            // Note: zone_letter is a UTM latitude band letter (C-X), not N/S hemisphere indicator.
            // Hemisphere is determined by the sign of the latitude.
            let zone = scirs2_spatial::projections::UTMZone {
                number: zone_num as u8,
                north: lat >= 0.0,
            };
            if let Ok((lat2, lon2)) = scirs2_spatial::projections::utm_to_geographic(easting, northing, zone) {
                // Roundtrip should preserve coordinates within tolerance
                assert_relative_eq!(lat, lat2, epsilon = 1e-6);
                assert_relative_eq!(lon, lon2, epsilon = 1e-6);
            }
        }
    }

    #[test]
    fn test_web_mercator_roundtrip(
        lat in -85.0..85.0,
        lon in -180.0..180.0
    ) {
        if let Ok((x, y)) = geographic_to_web_mercator(lat, lon) {
            let (lat2, lon2) = web_mercator_to_geographic(x, y);
            assert_relative_eq!(lat, lat2, epsilon = 1e-6);
            assert_relative_eq!(lon, lon2, epsilon = 1e-6);
        }
    }

    #[test]
    fn test_projection_preserves_order(
        lat1 in 0.0..80.0,
        lat2 in 0.0..80.0,
        lon in -180.0..180.0
    ) {
        prop_assume!(lat2 > lat1 + 0.1); // Ensure meaningful difference

        if let Ok((_, n1, _, _)) = geographic_to_utm(lat1, lon) {
            if let Ok((_, n2, _, _)) = geographic_to_utm(lat2, lon) {
                // Higher latitude should give higher northing
                prop_assert!(n2 > n1,
                    "UTM northing should increase with latitude");
            }
        }
    }
}

// Helper for slice operations
use scirs2_core::ndarray::s;

#[test]
fn test_property_based_framework() {
    // This test just ensures the property test framework is working
    assert_eq!(2 + 2, 4);
}
