rgeometry 0.12.0

use num_traits::*;
use rand::Rng;
use rand::distr::uniform::SampleUniform;
use rand::distr::{Distribution, StandardUniform};
use rand::seq::SliceRandom;
use std::cmp::Ordering;
use std::ops::*;

use crate::data::{Cursor, Point, PointLocation, TriangleView, Vector};
use crate::{Error, Orientation, PolygonScalar, TotalOrd};

use super::Polygon;

#[derive(Debug, Clone, Hash)]
pub struct PolygonConvex<T>(Polygon<T>);

///////////////////////////////////////////////////////////////////////////////
// PolygonConvex

impl<T> PolygonConvex<T>
where
  T: PolygonScalar,
{
  /// Assume that a polygon is convex.
  ///
  /// # Safety
  /// The input polygon has to be strictly convex, ie. no vertices are allowed to
  /// be concave or colinear.
  ///
  /// # Time complexity
  /// $O(1)$
  pub fn new_unchecked(poly: Polygon<T>) -> PolygonConvex<T> {
    PolygonConvex(poly)
  }

  /// Locate a point relative to a convex polygon.
  ///
  /// # Time complexity
  /// $O(\log n)$
  ///
  /// # Examples
  /// <iframe src="https://web.rgeometry.org/wasm/gist/2cb9ff5bd6ce24f395a5ea30280aabee"></iframe>
  ///
  pub fn locate(&self, pt: &Point<T, 2>) -> PointLocation {
    let poly = &self.0;
    let vertices = self.boundary_slice();
    let p0 = poly.point(vertices[0]);
    let mut lower = 1;
    let mut upper = vertices.len() - 1;
    while lower + 1 < upper {
      let middle = (lower + upper) / 2;
      if Point::orient(p0, poly.point(vertices[middle]), pt) == Orientation::CounterClockWise {
        lower = middle;
      } else {
        upper = middle;
      }
    }
    let p1 = poly.point(vertices[lower]);
    let p2 = poly.point(vertices[upper]);
    let triangle = TriangleView::new_unchecked([p0, p1, p2]);
    triangle.locate(pt)
  }

  /// Return the vertex that maximizes the dot product with `direction`.
  ///
  /// The search relies on the monotone sequence of outward-facing edge normals
  /// on a strictly convex polygon to cut the search interval in half until the
  /// supporting edge aligned with `direction` is found.
  ///
  /// If `direction` is the zero vector, the first boundary vertex is returned.
  ///
  /// # Time complexity
  /// $O(\log n)$
  ///
  /// # Examples
  /// ```rust
  /// # use rgeometry::data::{Point, Polygon, PolygonConvex, Vector};
  /// let polygon = PolygonConvex::new_unchecked(Polygon::new_unchecked(vec![
  ///   Point::new([0i32, 0]),
  ///   Point::new([2, 0]),
  ///   Point::new([2, 2]),
  ///   Point::new([0, 2]),
  /// ]));
  /// let cursor = polygon.extreme_in_direction(&Vector([1, 1]));
  /// assert_eq!(cursor.point(), &Point::new([2, 2]));
  /// ```
  #[must_use]
  pub fn extreme_in_direction(&self, direction: &Vector<T, 2>) -> Cursor<'_, T> {
    let polygon: &Polygon<T> = self.into();
    let vertices = polygon.boundary_slice();
    let n = vertices.len();
    assert!(
      n >= 3,
      "PolygonConvex::extreme_in_direction requires at least 3 vertices."
    );

    let direction_coords = &direction.0;
    let zero = T::from_constant(0);
    if direction_coords.iter().all(|coord| coord == &zero) {
      return polygon.cursor(vertices[0]);
    }

    let origin = [T::from_constant(0), T::from_constant(0)];

    let cursor_for_edge = |edge_idx: usize| -> Cursor<'_, T> {
      let vertex_idx = (edge_idx + 1) % n;
      polygon.cursor(vertices[vertex_idx])
    };

    // Get the endpoints of an edge (for computing the outward normal without overflow)
    let edge_endpoints = |edge_idx: usize| -> (&[T; 2], &[T; 2]) {
      let a = &polygon.point(vertices[edge_idx]).array;
      let b = &polygon.point(vertices[(edge_idx + 1) % n]).array;
      (a, b)
    };

    // Reference edge (edge 0)
    let (ref_start, ref_end) = edge_endpoints(0);

    // Compare the angle of edge `edge_idx`'s normal against the direction,
    // using edge 0's normal as the reference angle.
    // This uses overflow-safe primitives that don't compute normals explicitly.
    let compare_normal = |edge_idx: usize| -> Ordering {
      let (edge_start, edge_end) = edge_endpoints(edge_idx);
      // We need to compare:
      // - angle of edge_normal (normal of edge `edge_idx`) around origin
      // - angle of direction_coords around origin
      // using ref_normal (normal of edge 0) as the reference
      //
      // This is: ccw_cmp_around_with(ref_normal, origin, edge_normal, direction)
      //
      // We use the overflow-safe version that takes edge endpoints.
      Orientation::ccw_cmp_around_with_edge_normal_vs_direction(
        ref_start,
        ref_end,
        &origin,
        edge_start,
        edge_end,
        direction_coords,
      )
    };

    match compare_normal(0) {
      Ordering::Equal => return cursor_for_edge(0),
      Ordering::Greater => return cursor_for_edge(n - 1),
      Ordering::Less => {}
    }

    let cmp_last = compare_normal(n - 1);
    if cmp_last != Ordering::Greater {
      return cursor_for_edge(n - 1);
    }

    let mut lo = 0usize;
    let mut hi = n - 1;
    while lo + 1 < hi {
      let mid = (lo + hi) / 2;
      match compare_normal(mid) {
        Ordering::Greater => hi = mid,
        Ordering::Equal => return cursor_for_edge(mid),
        Ordering::Less => lo = mid,
      }
    }

    cursor_for_edge(lo)
  }

  /// Validates the following properties:
  ///  * Each vertex is convex, ie. not concave or colinear.
  ///  * All generate polygon properties hold true (eg. no duplicate points, no self-intersections).
  ///
  /// # Time complexity
  /// $O(n \log n)$
  pub fn validate(&self) -> Result<(), Error> {
    for cursor in self.0.iter_boundary() {
      if cursor.orientation() != Orientation::CounterClockWise {
        return Err(Error::ConvexViolation);
      }
    }
    self.0.validate()
  }

  pub fn cast<U>(self) -> PolygonConvex<U>
  where
    T: Clone + Into<U>,
    U: PolygonScalar,
  {
    PolygonConvex::new_unchecked(self.0.cast())
  }

  pub fn float(self) -> PolygonConvex<f64>
  where
    T: Clone + Into<f64>,
  {
    PolygonConvex::new_unchecked(self.0.float())
  }

  /// $O(1)$
  pub fn polygon(&self) -> &Polygon<T> {
    self.into()
  }

  /// Uniformly sample a random convex polygon.
  ///
  /// The output polygon is rooted in `(0,0)`, grows upwards, and has a height and width of [`T::max_value()`](Bounded::max_value).
  ///
  /// # Time complexity
  /// $O(n \log n)$
  ///
  // # Examples
  // ```no_run
  // # use rgeometry_wasm::playground::*;
  // # use rgeometry::data::*;
  // # clear_screen();
  // # set_viewport(2.0, 2.0);
  // # let convex: PolygonConvex<i8> = {
  // PolygonConvex::random(3, &mut rand::thread_rng())
  // # };
  // # render_polygon(&convex.float().normalize());
  // ```
  // <iframe src="https://web.rgeometry.org/wasm/gist/9abc54a5e2e3d33e3dd1785a71e812d2"></iframe>
  pub fn random<R>(n: usize, rng: &mut R) -> PolygonConvex<T>
  where
    T: Bounded + PolygonScalar + SampleUniform,
    R: Rng + ?Sized,
  {
    let n = n.max(3);
    loop {
      let vs = {
        let mut vs = random_vectors(n, rng);
        Vector::sort_around(&mut vs);
        vs
      };
      let vertices: Vec<Point<T, 2>> = vs
        .into_iter()
        .scan(Point::zero(), |st, vec| {
          *st += vec;
          Some((*st).clone())
        })
        .collect();
      let n_vertices = (*vertices).len();
      debug_assert_eq!(n_vertices, n);
      // FIXME: Use the convex hull algorithm for polygons rather than point sets.
      //        It runs in O(n) rather than O(n log n). Hasn't been implemented, yet, though.
      // If the vertices are all colinear then give up and try again.
      // FIXME: If the RNG always returns zero then we might loop forever.
      //        Maybe limit the number of recursions.
      if let Ok(p) = crate::algorithms::convex_hull(vertices) {
        return p;
      }
    }
  }
}

impl PolygonConvex<f64> {
  #[must_use]
  pub fn normalize(&self) -> PolygonConvex<f64> {
    PolygonConvex::new_unchecked(self.0.normalize())
  }
}

impl PolygonConvex<f32> {
  #[must_use]
  pub fn normalize(&self) -> PolygonConvex<f32> {
    PolygonConvex::new_unchecked(self.0.normalize())
  }
}

///////////////////////////////////////////////////////////////////////////////
// Trait Implementations

impl<T> Deref for PolygonConvex<T> {
  type Target = Polygon<T>;
  fn deref(&self) -> &Self::Target {
    &self.0
  }
}

impl<T> From<PolygonConvex<T>> for Polygon<T> {
  fn from(convex: PolygonConvex<T>) -> Polygon<T> {
    convex.0
  }
}

impl<'a, T> From<&'a PolygonConvex<T>> for &'a Polygon<T> {
  fn from(convex: &'a PolygonConvex<T>) -> &'a Polygon<T> {
    &convex.0
  }
}

impl Distribution<PolygonConvex<i64>> for StandardUniform {
  fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> PolygonConvex<i64> {
    PolygonConvex::random(100, rng)
  }
}

///////////////////////////////////////////////////////////////////////////////
// Helper functions

// Property: random_between(n, max, &mut rng).sum::<usize>() == max
fn random_between_iter<T, R>(n: usize, rng: &mut R) -> impl Iterator<Item = T> + use<T, R>
where
  T: PolygonScalar + Bounded + SampleUniform,
  R: Rng + ?Sized,
{
  let zero: T = T::from_constant(0);
  let max: T = Bounded::max_value();
  assert!(n > 0);
  let mut pts = Vec::with_capacity(n);
  while pts.len() < n - 1 {
    pts.push(rng.random_range(zero.clone()..max.clone()));
  }
  pts.sort_unstable_by(TotalOrd::total_cmp);
  pts.push(max);
  pts.into_iter().scan(zero, |from, x| {
    let out = x.clone() - (*from).clone();
    *from = x;
    Some(out)
  })
}

// Property: random_between_zero(10, 100, &mut rng).iter().sum::<isize>() == 0
fn random_between_zero<T, R>(n: usize, rng: &mut R) -> Vec<T>
where
  T: Bounded + PolygonScalar + SampleUniform,
  R: Rng + ?Sized,
{
  assert!(n >= 2);
  let n_positive = rng.random_range(1..n); // [1;n[
  let n_negative = n - n_positive;
  assert!(n_positive + n_negative == n);
  let positive = random_between_iter(n_positive, rng);
  let negative = random_between_iter(n_negative, rng).map(|i: T| -i);
  let mut result: Vec<T> = positive.chain(negative).collect();
  result.shuffle(rng);
  result
}

// Random vectors that sum to zero.
fn random_vectors<T, R>(n: usize, rng: &mut R) -> Vec<Vector<T, 2>>
where
  T: Bounded + PolygonScalar + SampleUniform,
  R: Rng + ?Sized,
{
  random_between_zero(n, rng)
    .into_iter()
    .zip(random_between_zero(n, rng))
    .map(|(a, b)| Vector([a, b]))
    .collect()
}

///////////////////////////////////////////////////////////////////////////////
// Tests

#[cfg(test)]
mod tests {
  use super::*;
  use rand::SeedableRng;
  use rand::rngs::SmallRng;

  use proptest::prelude::*;
  use proptest::proptest as proptest_block;

  proptest_block! {
    // These traits are usually derived but let's not rely on that.
    #[test]
    fn fuzz_convex_traits(poly: PolygonConvex<i8>) {
      let _ = format!("{:?}", &poly);
      let _ = poly.clone();
    }

    #[test]
    fn fuzz_validate(poly: Polygon<i8>) {
      let convex = PolygonConvex::new_unchecked(poly);
      let _ = convex.validate();
    }

    #[test]
    fn all_random_convex_polygons_are_valid_i8(poly: PolygonConvex<i8>) {
      prop_assert_eq!(poly.validate().err(), None)
    }

    #[test]
    fn random_convex_prop(poly: PolygonConvex<i8>) {
      let (min, max) = poly.bounding_box();
      prop_assert_eq!(min.y_coord(), &0);
      let width = max.x_coord() - min.x_coord();
      let height = max.y_coord() - min.y_coord();
      prop_assert_eq!(width, i8::MAX);
      prop_assert_eq!(height, i8::MAX);
    }

    #[test]
    fn all_random_convex_polygons_are_valid_i64(poly: PolygonConvex<i64>) {
      prop_assert_eq!(poly.validate().err(), None)
    }

    #[test]
    fn sum_to_max(n in 1..1000, seed: u64) {
      let mut rng = SmallRng::seed_from_u64(seed);
      let vecs = random_between_iter::<i8, _>(n as usize, &mut rng);
      prop_assert_eq!(vecs.sum::<i8>(), i8::MAX);

      let vecs = random_between_iter::<i64, _>(n as usize, &mut rng);
      prop_assert_eq!(vecs.sum::<i64>(), i64::MAX);
    }

    #[test]
    fn random_between_zero_properties(n in 2..1000, seed: u64) {
      let mut rng = SmallRng::seed_from_u64(seed);
      let vecs: Vec<i8> = random_between_zero(n as usize, &mut rng);
      prop_assert_eq!(vecs.iter().sum::<i8>(), 0);
      prop_assert_eq!(vecs.len(), n as usize);

      let vecs: Vec<i64> = random_between_zero(n as usize, &mut rng);
      prop_assert_eq!(vecs.iter().sum::<i64>(), 0);
      prop_assert_eq!(vecs.len(), n as usize);
    }

    #[test]
    fn sum_to_zero_vector(n in 2..1000, seed: u64) {
      let mut rng = SmallRng::seed_from_u64(seed);
      let vecs: Vec<Vector<i8, 2>> = random_vectors(n as usize, &mut rng);
      prop_assert_eq!(vecs.into_iter().sum::<Vector<i8, 2>>(), Vector([0, 0]))
    }

    #[test]
    fn extreme_direction_matches_naive(poly: PolygonConvex<i16>, direction: Vector<i16, 2>) {
      let cursor = poly.extreme_in_direction(&direction);
      let expected = poly
        .iter_boundary()
        .max_by(|a, b| direction.cmp_along(a.point(), b.point()))
        .unwrap();
      prop_assert_eq!(cursor.point(), expected.point());
    }

    #[test]
    fn extreme_direction_matches_naive_after_rotation(
      poly: PolygonConvex<i16>,
      direction: Vector<i16, 2>,
      rotate_by: usize,
    ) {
      let rotated = rotate_polygon(poly, rotate_by);
      let cursor = rotated.extreme_in_direction(&direction);
      let expected = rotated
        .iter_boundary()
        .max_by(|a, b| direction.cmp_along(a.point(), b.point()))
        .unwrap();
      prop_assert_eq!(cursor.point(), expected.point());
    }
  }

  fn rotate_polygon(poly: PolygonConvex<i16>, rotate_by: usize) -> PolygonConvex<i16> {
    let n = poly.boundary_slice().len();
    let shift = rotate_by % n;
    if shift == 0 {
      return poly;
    }
    let points: Vec<Point<i16, 2>> = poly.iter_boundary().map(|c| *c.point()).collect();
    let rotated: Vec<Point<i16, 2>> = points.iter().cycle().skip(shift).take(n).cloned().collect();
    PolygonConvex::new_unchecked(Polygon::new_unchecked(rotated))
  }

  #[test]
  fn extreme_direction_zero_vector_returns_first_vertex() {
    let polygon = sample_square();
    let cursor = polygon.extreme_in_direction(&Vector([0, 0]));
    let first = polygon.boundary_slice()[0];
    assert_eq!(cursor.point(), polygon.point(first));
  }

  #[test]
  fn extreme_direction_prefers_vertex_after_supporting_edge() {
    let polygon = sample_square();
    let cursor = polygon.extreme_in_direction(&Vector([1, 0]));
    let expected = polygon.boundary_slice()[2];
    assert_eq!(cursor.point(), polygon.point(expected));
  }

  fn sample_square() -> PolygonConvex<i32> {
    PolygonConvex::new_unchecked(Polygon::new_unchecked(vec![
      Point::new([0, 0]),
      Point::new([2, 0]),
      Point::new([2, 2]),
      Point::new([0, 2]),
    ]))
  }

  #[test]
  fn extreme_direction_diagonal_on_square() {
    // This is the doctest case - direction [1, 1] on a square
    let polygon = sample_square();
    let direction = Vector([1, 1]);
    let cursor = polygon.extreme_in_direction(&direction);
    // The point that maximizes dot product with [1,1] is [2,2]
    // because 2*1 + 2*1 = 4 is the maximum dot product
    assert_eq!(cursor.point(), &Point::new([2, 2]));
  }

  // Test case demonstrating potential overflow in extreme_in_direction for i8.
  // The issue: when computing edge vectors like (127 - (-127)) = 254, this overflows i8.
  // Similarly, computing the normal [-edge_y, edge_x] can overflow.
  #[test]
  fn extreme_direction_i8_overflow_case() {
    // A wide triangle that spans nearly the full i8 range.
    // Edge from (-100, 0) to (100, 0) has vector (200, 0) which overflows i8.
    let polygon: PolygonConvex<i8> = PolygonConvex::new_unchecked(Polygon::new_unchecked(vec![
      Point::new([-100i8, 0]),
      Point::new([100i8, 0]),
      Point::new([0i8, 100]),
    ]));
    let direction = Vector([1i8, 0]);
    let cursor = polygon.extreme_in_direction(&direction);
    // The rightmost point should be (100, 0)
    assert_eq!(cursor.point(), &Point::new([100i8, 0]));
  }

  #[test]
  fn extreme_direction_i8_matches_naive() {
    // Test that extreme_in_direction gives the same result as the naive O(n) approach
    // for a polygon that might cause overflow.
    let polygon: PolygonConvex<i8> = PolygonConvex::new_unchecked(Polygon::new_unchecked(vec![
      Point::new([-100i8, 0]),
      Point::new([100i8, 0]),
      Point::new([0i8, 100]),
    ]));
    let direction = Vector([1i8, 1]);
    let cursor = polygon.extreme_in_direction(&direction);
    let expected = polygon
      .iter_boundary()
      .max_by(|a, b| direction.cmp_along(a.point(), b.point()))
      .unwrap();
    assert_eq!(cursor.point(), expected.point());
  }
}