curve_sampling/

lib.rs

//! This module provide a collection of routines to perform adaptive
//! sampling of curves as well as manipulating those samplings.
//!
//! As of now, the following is available:
//! - uniform sampling of the graph of functions ℝ → ℝ and the image
//!   of functions ℝ → ℝ² (see [`Sampling::uniform`]);
//! - adaptive sampling of the graph functions ℝ → ℝ and the image of
//!   functions ℝ → ℝ² (see [`Sampling::fun`]).
//!
//! Samplings can be saved as a list of coordinates x-y, one point per
//! line with blank lines to indicate that the path is interrupted,
//! with [`Sampling::write`].  They can also be saved in a format
//! usable from the [LaTeX][] package [TikZ][] using
//! [`Sampling::latex`].  This format allows to add arrows to the
//! curves to indicate their orientation.
//!
//! [LaTeX]: https://www.latex-project.org/
//! [TikZ]: https://tikz.dev/

use std::{
    cell::RefCell,
    cmp::Ordering,
    fmt::{self, Display, Formatter},
    io::{self, Write},
    marker::PhantomData,
    ops::ControlFlow,
};
use rgb::*;

// mod fibonacci_heap;
// use fibonacci_heap as pq;
// mod priority_queue;
// use priority_queue as pq;
mod approx_priority_queue;
use approx_priority_queue as pq;

mod linked_list;
use linked_list as list;
use list::List;

////////////////////////////////////////////////////////////////////////
//
// Bounding box

/// A two dimensional rectangle \[`xmin`, `xmax`\] × \[`ymin`, `ymax`\].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct BoundingBox {
    pub xmin: f64,
    pub xmax: f64,
    pub ymin: f64,
    pub ymax: f64,
}

impl BoundingBox {
    #[inline]
    pub(crate) fn empty() -> Self {
        Self { xmin: f64::INFINITY,  xmax: f64::NEG_INFINITY,
               ymin: f64::INFINITY,  ymax: f64::NEG_INFINITY }
    }

    /// Return `true` if the bounding box has a non-empty interior.
    ///
    /// # Example
    ///
    /// ```
    /// use curve_sampling::Sampling;
    /// let s = Sampling::from_iter([0.]);
    /// assert!(s.bounding_box().is_empty());
    /// ```
    #[inline]
    pub fn is_empty(&self) -> bool {
        !(self.xmin < self.xmax && self.ymin < self.ymax) // NAN ⟹ empty
    }

    /// Return `true` if the point `p` belongs to the box `self`
    /// (possibly on the boundary).
    #[inline]
    fn contains<D>(&self, p: &Point<D>) -> bool {
        let [x, y] = p.xy;
        self.xmin <= x && x <= self.xmax && self.ymin <= y && y <= self.ymax
    }

    /// Return the smaller bounding-box containing both `self` and
    /// `other`.
    #[inline]
    pub fn hull(&self, other: &Self) -> Self {
        BoundingBox { xmin: self.xmin.min(other.xmin),
                      xmax: self.xmax.max(other.xmax),
                      ymin: self.ymin.min(other.ymin),
                      ymax: self.ymax.max(other.ymax) }
    }
}


////////////////////////////////////////////////////////////////////////
//
// Sampling datastructure

/// A 2D point with coordinates (`x`, `y`) supposed to be given by a
/// function evaluated at "time" `t`.  A path is a sequence of points.
/// Points that are invalid (see [`Point::is_valid`]) are considered
/// cuts in the path (the line is interrupted).  As most 2 cuts may
/// follow each other (to track ranges of t outside the domain of the
/// function) except at the boundary of the path where at most 1 cut
/// is allowed.
#[derive(Debug)]
struct Point<D> {
    t: f64, // "time" parameter, ALWAYS finite.
    xy: [f64; 2],
    data: D,
    cost: f64, // Cache the cost of the point (not the segment).  If
               // the point is not valid, the cost has no meaning.

    // Pointer to the priority queue node if the point is the initial
    // point of a segment.  `None` otherwise.
    witness: RefCell<Option<pq::Witness<list::Witness<Point<D>>>>>,
}

/// Essential elements of a point for the algorithm.  This is to pass
/// to functions to work around the fact that `Point` is not copy
/// (when we borrow a point, we implicitly borrow the whole path).
//
// It is actually slower to use `Coord` in `Point`.
#[derive(Clone, Copy, Debug)]
struct Coord {
    t: f64,
    xy: [f64; 2],
    cost: f64,
}

impl<D: Clone> Clone for Point<D> {
    // As a cloned point may likely go to another path, reset the
    // priority queue pointer.
    fn clone(&self) -> Self {
        Self {
            t: self.t,  xy: self.xy,
            data: self.data.clone(),
            cost: self.cost,
            witness: None.into(),
        }
    }
}

impl<D> Point<D> {
    /// Return a new point.  `t` is assumed to be finite.
    #[inline]
    fn new_unchecked(t: f64, xy: [f64; 2], data: D) -> Self {
        Point { t, xy, data,  cost: 0.,
            witness: None.into() }
    }

    /// Return a point considered to be a cut in the path.
    #[inline]
    fn cut(t: f64, data: D) -> Self {
        Point { t, xy: [f64::NAN; 2], data,
                cost: 0.,
            witness: None.into() }
    }

    /// Return `true` if the point is valid — otherwise it is a cut.
    #[inline]
    fn is_valid(&self) -> bool {
        self.xy.iter().all(|z| z.is_finite())
    }

    fn cmp_t_incr(p1: &Self, p2: &Self) -> Ordering {
        p1.t.total_cmp(&p2.t)
    }

    fn cmp_t_decr(p1: &Self, p2: &Self) -> Ordering {
        p1.t.total_cmp(&p2.t).reverse()
    }

    /// Transform the data as optional.
    #[inline]
    fn opt(self) -> Point<Option<D>> {
        Point {
            t: self.t,  xy: self.xy,
            data: Some(self.data),
            cost: self.cost,
            witness: None.into(),
        }
    }

    /// Project on the essential information for cost computation.
    ///
    /// This is used to work around the fact that Point<_> can not be
    /// copy because <T> and the `witness`.
    #[inline]
    fn coord(&self) -> Coord {
        Coord { t: self.t, xy: self.xy, cost: self.cost }
    }
}

/// A sampling of a 2D curve.  This can be thought of as a 2D path,
/// with possible "cuts" due to discontinuities or leaving the domain
/// of the (parametric) function describing the path.
///
/// Additional data of type `D` may be attached to each point.
pub struct Sampling<D> {
    // Priority queue of segments.  The value points to the first
    // `Point` of the segment (thus its .next pointer must not be
    // null).  At least one of the endpoint of all segments must be
    // valid (see `Point::is_valid`) If `pq` is empty but `begin` is
    // not null, it means that the costs need to to be updated.
    pq: PQ<D>,
    path: List<Point<D>>,
    vp: Option<BoundingBox>, // viewport (zone of interest)
}

type PQ<D> = pq::PQ<list::Witness<Point<D>>>;

impl<D: Clone> Clone for Sampling<D> {
    fn clone(&self) -> Self {
        // It is guessed that one clones the sampling to transform it.
        // Thus start with an empty queue.
        Self { pq: PQ::new(),
               path: self.path.clone(),
               vp: self.vp }
    }
}

/// `t` is the length of total range of time, `x` and `y` are the
/// dimensions of the bounding box.
#[derive(Debug, Clone, Copy)]
struct Lengths {
    t: f64,
    x: f64,
    y: f64,
}

impl<D> Sampling<D> {
    /// Return `true` if the sampling contains no point.
    #[inline]
    pub fn is_empty(&self) -> bool { self.path.is_empty() }

    /// Create an empty sampling.
    #[inline]
    pub(crate) fn empty() -> Self {
        Self { pq: PQ::new(),
               path: List::new(),
               vp: None }
    }

    #[inline]
    pub(crate) fn singleton(p: Point<D>) -> Self {
        debug_assert!(p.t.is_finite() && p.xy.iter().all(|z| z.is_finite()));
        let mut path = List::new();
        path.push_back(p);
        Self { pq: PQ::new(), path, vp: None }
    }

    fn from_list(path: List<Point<D>>) -> Self {
        Self {
            pq: PQ::new(),
            path,
            vp: None,
        }
    }

    #[inline]
    pub(crate) fn set_vp(&mut self, bb: BoundingBox) {
        self.vp = Some(bb);
    }

    /// Return the length of the "time interval" as well as the
    /// lengths of the viewport.
    pub(crate) fn lengths(&self) -> Option<Lengths> {
        if self.is_empty() { return None }
        let p0 = self.path.first().unwrap();
        let p1 = self.path.last().unwrap();
        let (len_x, len_y) = match self.vp {
            Some(vp) => (vp.xmax - vp.xmin, vp.ymax - vp.ymin),
            None => (1., 1.),

        };
        Some(Lengths { t: (p1.t - p0.t).abs(), x: len_x, y: len_y })
    }

    /// Return the number of points and "cuts" in the sampling.
    pub fn len(&self) -> usize {
        self.path.len()
    }

    /// Returns an iterator on the points (and cuts) of the path.
    /// More precisely, a path is made of continuous segments whose
    /// points are given by contiguous values `[x,y]` with both `x`
    /// and `y` not NaN, interspaced by "cuts" `[f64::NAN; 2]`.  Two
    /// cuts never follow each other.  Isolated points `p` are given
    /// by ... `[f64::NAN; 2]`, `p`, `None`,...
    pub fn iter(&self) -> Iter<'_, D> {
        Iter { iter: self.iter_data() }
    }

    /// Same as [`Self::iter`] but also provides a reference to the
    /// data associated to each point.
    pub fn iter_data(&self) -> IterData<'_, D> {
        IterData {
            path: self.path.iter(),
            prev_is_cut: true,
        }
    }

    /// Returns an iterator that allows to modify the points and cuts
    /// of the path.  Unlike [`Self::iter`], this iterates on all the
    /// nodes even if several cuts (i.e., node with a non finite
    /// coordinate) follow each other.  The data associated to each
    /// point is unmodified.
    pub fn iter_mut(&mut self) -> IterMut<'_, D> {
        IterMut { path: self.path.iter_mut() }
    }

    /// Consumes the sampling and return an iterator on the curve
    /// points and (owned) associated data.
    pub fn into_iter_data(self) -> IntoIterData<D> {
        IntoIterData { path: self.path.into_iter() }
    }

    /// Return the X-coordinates of the sampling.
    #[inline]
    pub fn x(&self) -> Vec<f64> {
        let mut v = Vec::with_capacity(self.path.len());
        for [x, _] in self.iter() { v.push(x) }
        v
    }

    /// Return the Y-coordinates of the sampling.
    #[inline]
    pub fn y(&self) -> Vec<f64> {
        let mut v = Vec::with_capacity(self.path.len());
        for [_, y] in self.iter() { v.push(y) }
        v
    }
}


/// Iterator on the curve points (and cuts) together with the attached data.
///
/// Created by [`Sampling::iter_data`].
pub struct IterData<'a, D> {
    path: list::Iter<'a, Point<D>>,
    prev_is_cut: bool,
}

impl<'a, D> Iterator for IterData<'a, D> {
    type Item = ([f64; 2], &'a D);

    fn next(&mut self) -> Option<Self::Item> {
        match self.path.next() {
            None => None,
            Some(p) => {
                if p.is_valid() {
                    self.prev_is_cut = false;
                    Some((p.xy, &p.data))
                } else if self.prev_is_cut {
                    // Find the next valid point.
                    // FIXME: use `try_find` once stabilized.
                    let r = self.path.try_fold((), |_, p| {
                        if p.is_valid() {
                            ControlFlow::Break(((), p))
                        } else {
                            ControlFlow::Continue(())
                        }
                    });
                    match r {
                        ControlFlow::Continue(_) => {
                            // Iterator exhausted
                            None
                        }
                        ControlFlow::Break((_, p)) => {
                            self.prev_is_cut = false;
                            Some((p.xy, &p.data))
                        }
                    }
                } else {
                    self.prev_is_cut = true;
                    Some(([f64::NAN; 2], &p.data))
                }
            }
        }
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        (0, Some(self.path.len()))
    }
}

/// Iterator on the curve points (and cuts).
///
/// Created by [`Sampling::iter`].
pub struct Iter<'a, D> {
    iter: IterData<'a, D>,
}

impl<'a, D> Iterator for Iter<'a, D> {
    type Item = [f64; 2];

    #[inline]
    fn next(&mut self) -> Option<Self::Item> {
        self.iter.next().map(|(xy, _)| xy)
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        self.iter.size_hint()
    }
}


/// Mutable iterator on the curve points (and cuts).
///
/// Created by [`Sampling::iter_mut`].
pub struct IterMut<'a, D> {
    path: list::IterMut<'a, Point<D>>,
}

impl<'a, D> Iterator for IterMut<'a, D> {
    type Item = &'a mut [f64; 2];

    fn next(&mut self) -> Option<Self::Item> {
        self.path.next().map(|p| &mut p.xy)
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        (self.path.len(), Some(self.path.len()))
    }
}

/// Iterator returning the curve points and owned data.
///
/// Created by [`Sampling::into_iter_data`].
pub struct IntoIterData<D> {
    path: list::IntoIter<Point<D>>,
}

impl<D> Iterator for IntoIterData<D> {
    type Item = ([f64; 2], D);

    fn next(&mut self) -> Option<Self::Item> {
         self.path.next().map(|p| (p.xy, p.data))
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        (self.path.len(), Some(self.path.len()))
    }
}


/// Intersection of a segment with the bounding box.
#[derive(Debug)]
enum Intersection<D> {
    Empty,
    Pt(Point<D>),
    Seg(Point<D>, Point<D>),
}

impl<D> Sampling<D> {
    /// Return the smallest rectangle enclosing all the points of the
    /// sampling `self`.  If the path is empty, the "min" fields of
    /// the bounding box are set to +∞ and "max" fields to -∞.
    ///
    /// # Examples
    ///
    /// ```
    /// use curve_sampling::{Sampling, BoundingBox};
    /// let s = Sampling::from_iter([3., 4.]);
    /// let bb = BoundingBox { xmin: 0., xmax: 1., ymin: 3., ymax: 4. };
    /// assert_eq!(s.bounding_box(), bb);
    /// ```
    ///
    /// ```
    /// use curve_sampling::Sampling;
    /// let s = Sampling::from_iter([0.]);
    /// assert!(s.bounding_box().is_empty());
    /// ```
    pub fn bounding_box(&self) -> BoundingBox {
        let mut points = self.iter().skip_while(|[x,_]| x.is_nan());
        let mut bb = match &points.next() {
            Some([x, y]) => BoundingBox { xmin: *x,  xmax: *x,
                                         ymin: *y,  ymax: *y },
            None => return BoundingBox::empty()
        };
        for [x, y] in points {
            if x < bb.xmin { bb.xmin = x }
            else if bb.xmax < x { bb.xmax = x };
            if y < bb.ymin { bb.ymin = y }
            else if bb.ymax < y { bb.ymax = y };
        }
        bb
    }

    /// Transpose in place the x and y coordinates of the sampling.
    pub fn transpose(&mut self) -> &mut Self {
        for p in self.path.iter_mut() {
            let [x, y] = p.xy;
            p.xy = [y, x];
        }
        self
    }

    /// Assume `p0` ∈ `bb` and `p1` ∉ `bb`.  Return the point
    /// intersecting the boundary of `bb`.  If the intersection point
    /// is the same as `p0`, return `None`
    #[inline]
    #[must_use]
    fn intersect(
        p0: Coord, p1: Coord, bb: BoundingBox
    ) -> Option<Point<Option<D>>> {
        let [x0, y0] = p0.xy;
        let [x1, y1] = p1.xy;
        let mut t = 1.; // t ∈ [0, 1]
        let dx = x1 - x0; // May be 0.
        let r = (if dx >= 0. {bb.xmax} else {bb.xmin} - x0) / dx;
        if r < t { t = r } // ⟹ r finite (as r ≥ 0 or NaN)
        let dy = y1 - y0; // May be 0.
        let r = (if dy >= 0. {bb.ymax} else {bb.ymin} - y0) / dy;
        if r < t { t = r };
        if t <= 1e-14 {
            None
        } else {
            Some(Point::new_unchecked(
                p0.t + t * (p1.t - p0.t),
                [x0 + t * dx, y0 + t * dy],
            None))
        }
    }

    /// Assume `p0` ∉ `bb` and `p1` ∉ `bb` (thus, not even on the
    /// boundary of `bb`) and `p0` ≠ `p1`.  Return the endpoints of
    /// the segment intersecting the boundary of `bb` if any.  The
    /// "parameter direction" of `p0`, `p1` is preserved.
    #[inline]
    #[must_use]
    fn intersect_seg(
        p0: Coord, p1: Coord, bb: BoundingBox,
    ) -> Intersection<Option<D>> {
        let [x0, y0] = p0.xy;
        let [x1, y1] = p1.xy;
        let mut t0 = 0.; // t0 ∈ [0, 1]
        let mut t1 = 1.; // t1 ∈ [0, 1]
        let dx = x1 - x0; // may be 0.
        // r0 ≤ r1 or NAN if on x-boundary lines.
        let (r0, r1) = if dx >= 0. {
            ((bb.xmin - x0) / dx, (bb.xmax - x0) / dx)
        } else {
            ((bb.xmax - x0) / dx, (bb.xmin - x0) / dx)
        };
        if r0 > 1. || r1 < 0. { return Intersection::Empty }
        if r0 > 0. { t0 = r0 } // if r0 is NAN, keep the whole segment
        if r1 < 1. { t1 = r1 }
        let dy = y1 - y0; // may be 0.
        let (r0, r1) = if dy >= 0. {
            ((bb.ymin - y0) / dy, (bb.ymax - y0) / dy)
        } else {
            ((bb.ymax - y0) / dy, (bb.ymin - y0) / dy)
        };
        if r0 > t1 || r1 < t0 { return Intersection::Empty }
        if r0 > t0 { t0 = r0 }
        if r1 < t1 { t1 = r1 }
        if t0 < t1 { // segment not reduced to a point
            let dt = p1.t - p0.t;
            let q0 = Point::new_unchecked(p0.t + t0 * dt,
                [x0 + t0 * dx, y0 + t0 * dy], None);
            let q1 = Point::new_unchecked(p0.t + t1 * dt,
                [x0 + t1 * dx, y0 + t1 * dy], None);
            Intersection::Seg(q0, q1)
        } else if t0 == t1 {
            let q0 = Point::new_unchecked(
                p0.t + t0 * (p1.t - p0.t),
                [x0 + t0 * dx, y0 + t0 * dy], None);
            Intersection::Pt(q0)
        } else {
            Intersection::Empty
        }
    }

    /// Returns the sampling `self` but clipped to the 2D box `bb`.  A
    /// path that crosses the boundary will get additional nodes at
    /// the points of crossing and the part outside the bounding box
    /// will be dropped.  (Thus a path entirely out of the bounding
    /// box will be removed.)
    #[must_use]
    pub fn clip(self, bb: BoundingBox) -> Sampling<Option<D>> {
        let mut path = Vec::new();
        // First point of the current segment, if any.
        let mut p0_opt: Option<Coord> = None;
        let mut p0_inside = false;
        let mut prev_cut = true; // New path ("virtual" cut)
        for p1 in self.path.into_iter() {
            if prev_cut {
                // A cut was pushed at an earlier step (and the path
                // has not subsequently intersected `bb`).  This may
                // be because the original path was just cut (`p0_opt
                // = None`) or because the previous segment was cut
                // (`p0_opt = Some(p0)` with `p0` ∉ `bb`) or because
                // there was an earlier cut and the next point `p0` is
                // the possible new start.
                match (p0_opt, p1.is_valid()) {
                    (Some(p0), true) => {
                        let p1_inside = bb.contains(&p1);
                        p0_opt = Some(p1.coord());
                        if p0_inside { // p0 ∈ bb with cut before
                            // p0 is already in `s`.
                            if p1_inside {
                                path.push(p1.opt());
                                prev_cut = false;
                            } else if
                                let Some(p) = Self::intersect(p0, p1.coord(), bb) {
                                    let t = p.t;
                                    path.push(p);
                                    path.push(Point::cut(t, None));
                                } else {
                                    let c = Point::cut(p0.t, None);
                                    path.push(c);
                                }
                        } else if p1_inside { // p0 ∉ bb, p1 ∈ bb
                            if let Some(p) = Self::intersect(p1.coord(), p0, bb) {
                                path.push(p); // p ≠ p1
                            }
                            path.push(p1.opt());
                            prev_cut = false;
                        } else { // p0, p1 ∉ bb but maybe intersection
                            match Self::intersect_seg(p0, p1.coord(), bb) {
                                Intersection::Seg(q0, q1) => {
                                    let t1 = q1.t;
                                    path.push(q0);
                                    path.push(q1);
                                    path.push(Point::cut(t1, None));
                                }
                                Intersection::Pt(p) => {
                                    let t = p.t;
                                    path.push(p);
                                    path.push(Point::cut(t, None));
                                }
                                Intersection::Empty => (),
                            }
                        }
                        p0_inside = p1_inside;
                    }
                    (None, true) => {
                        p0_inside = bb.contains(&p1);
                        p0_opt = Some(p1.coord());
                        if p0_inside {
                            // We keep `p0`.
                            path.push(p1.opt());
                        }
                    }
                    (Some(p0), false) => {
                        if p0_inside {
                            // `p0` is already in the `s.path`.
                            path.push(Point::cut(p0.t, None));
                        }
                        p0_opt = None;
                    }
                    (None, false) => (), // p0_opt == None
                }
            } else {
                // Previous step was not a cut which also implies that
                // `p0_opt = Some(p0)` with `p0` ∈ `bb`, and `p0` is
                // already present in the final sampling.
                let p0 = p0_opt.expect("No cut ⟹ previous point");
                debug_assert!(p0_inside);
                if p1.is_valid() {
                    p0_opt = Some(p1.coord());
                    p0_inside = bb.contains(&p1); // update for next step
                    if p0_inside { // p0, p1 ∈ bb
                        path.push(p1.opt());
                    } else { // p0 ∈ bb, p1 ∉ bb
                        if let Some(p) = Self::intersect(p0, p1.coord(), bb) {
                            let t = p.t;
                            path.push(p);
                            path.push(Point::cut(t, None));
                        } else {
                            path.push(Point::cut(p0.t, None));
                        }
                        prev_cut = true;
                    }
                } else { // p1 is invalid (i.e., represent a cut)
                    p0_opt = None;
                    path.push(p1.opt());
                    prev_cut = true
                }
            }
        }
        if prev_cut { path.pop(); }
        let mut s = Sampling::from_list(path.into());
        s.set_vp(bb);
        s
    }

    /// Create a sampling from an iterator of points.
    ///
    /// # Preconditions
    /// - It is assumed that the `Point.t` field is sorted in increasing
    ///   or decreasing order.
    /// - The invariant "`p.y` is finite ⇒ `p.x` is finite" is assumed.
    fn from_point_iterator<P>(points: P) -> Self
    where P: IntoIterator<Item = Point<D>> {
        let mut s = Sampling::empty();
        let mut points = points.into_iter();
        macro_rules! skip_until_last_cut { () => {
            let mut cut = None;
            let mut first_pt = None;
            while let Some(p) = points.next() {
                if p.is_valid() { first_pt = Some(p); break; }
                cut = Some(p);
            }
            match (cut, first_pt) {
                (_, None) => {
                    return s
                }
                (None, Some(p)) => {
                    s.path.push_back(p);
                }
                (Some(c), Some(p)) => {
                    s.path.push_back(Point::cut(c.t, c.data));
                    s.path.push_back(p);
                }
            }
        }}
        skip_until_last_cut!();
        while let Some(p) = points.next() {
            if p.is_valid() {
                s.path.push_back(p);
            } else {
                s.path.push_back(Point::cut(p.t, p.data));
                skip_until_last_cut!();
            }
        }
        s
    }
}

impl<D, Y> FromIterator<Y> for Sampling<D>
where Y: Img<D> {
    /// Return an sampling from the Y-values (possibly with data).
    /// The X-values are the indices of the elements (see
    /// [`Iterator::enumerate`]).  Points with non-finite coordinates
    /// are interpreted as cuts.
    fn from_iter<P>(points: P) -> Self
    where P: IntoIterator<Item = Y> {
        Sampling::from_point_iterator(
            points.into_iter().enumerate().map(|(i, y)| {
                let t = i as f64;
                let p = y.into_point(t);
                let [x, y] = p.xy;
                if x.is_finite() && y.is_finite() {
                    p
                } else {
                    Point::cut(t, p.data)
                } }))
    }
}

struct MergeSortedIterator<I1, I2, CMP>
where I1: Iterator,
      I2: Iterator,
{
    i1: std::iter::Peekable<I1>,
    i2: std::iter::Peekable<I2>,
    cmp: CMP,
}

impl<I1, I2, CMP> MergeSortedIterator<I1, I2, CMP>
where I1: Iterator,
      I2: Iterator,
      CMP: Fn(&I1::Item, &I2::Item) -> Ordering,
{
    /// Create a new iterator returning the values of `i1` and `i2` in
    /// the order given by `cmp`.
    ///
    /// # Precondition
    /// The iterators `i1` and `i2` are assumed to be sorted, both
    /// with the order `cmp`.
    fn new(i1: I1, i2: I2, cmp: CMP) -> Self {
        Self {
            i1: i1.peekable(),
            i2: i2.peekable(),
            cmp,
        }
    }
}

impl<I1, I2, CMP> Iterator for MergeSortedIterator<I1, I2, CMP>
where I1: Iterator,
      I2: Iterator<Item = I1::Item>,
      CMP: Fn(&I1::Item, &I1::Item) -> Ordering,
{
    type Item = I1::Item;

    fn next(&mut self) -> Option<Self::Item> {
        if let Some(p1) = self.i1.peek() {
            if let Some(p2) = self.i2.peek() {
                match (self.cmp)(p1, p2) {
                    Ordering::Less => self.i1.next(),
                    Ordering::Greater => self.i2.next(),
                    Ordering::Equal => {
                        // For a given time, the data should be equal.
                        drop(self.i2.next());
                        self.i1.next()
                    },
                }
            } else {
                self.i1.next()
            }
        } else {
            self.i2.next()
        }
    }
}

////////////////////////////////////////////////////////////////////////
//
// Acceptable return values to generate samplings, possibly with
// attached data.

/// Suitable images of functions that can be used to generate samplings.
///
/// See e.g. [`Sampling::uniform`] and [`Sampling::fun`].
//
// This trait is not implementable from the outside so it can uphold
// `Point` invariants for example `Point.t` must be the `t` provided
// and be finite,...
pub trait Img<D> {
    /// Transform a map `f64` → `self` into the internal value `Point`.
    /// It will only be used with finite `t`.
    #[allow(private_interfaces)]
    fn into_point(self, t: f64) -> Point<D>;
}

/// Type indicating that no data is associated with the sampling.
#[derive(Clone, Debug)]
pub struct NoData {
    // Non-constructible outside this library (i.e. can only be used
    // as a type annotation).
    marker: PhantomData<()>,
}

impl NoData {
    fn new() -> Self { Self { marker: PhantomData } }
}

impl Img<NoData> for f64 {
    #[allow(private_interfaces)]
    #[inline]
    fn into_point(self, t: f64) -> Point<NoData> {
        Point::new_unchecked(t, [t, self], NoData::new())
    }
}

impl Img<NoData> for [f64; 2] {
    #[allow(private_interfaces)]
    #[inline]
    fn into_point(self, t: f64) -> Point<NoData> {
        Point::new_unchecked(t, self, NoData::new())
    }
}

impl Img<NoData> for (f64, f64) {
    #[allow(private_interfaces)]
    #[inline]
    fn into_point(self, t: f64) -> Point<NoData> {
        Point::new_unchecked(t, [self.0, self.1], NoData::new())
    }
}

#[cfg(feature = "num-complex")]
impl Img<num_complex::Complex<f64>> for num_complex::Complex<f64> {
    #[allow(private_interfaces)]
    #[inline]
    fn into_point(self, t: f64) -> Point<num_complex::Complex<f64>> {
        Point::new_unchecked(t, [self.re, self.im], self)
    }
}

/// Wrapper to indicate that a value is a data associated with a point.
///
/// See [`Sampling::fun`] for an example.
pub struct Data<D>(pub D);

/// Values `(y, Data(d))` are acceptable images (see [`Sampling::fun`]).
impl<D> Img<D> for (f64, Data::<D>) {
    #[allow(private_interfaces)]
    #[inline]
    fn into_point(self, t: f64) -> Point<D> {
        Point::new_unchecked(t, [t, self.0], self.1.0)
    }
}

/// Values `([x,y], Data(d))` are acceptable images (see [`Sampling::fun`]).
impl<D> Img<D> for ([f64; 2], Data::<D>) {
    #[allow(private_interfaces)]
    #[inline]
    fn into_point(self, t: f64) -> Point<D> {
        Point::new_unchecked(t, self.0, self.1.0)
    }
}

#[cfg(feature = "num-complex")]
impl<D> Img<D> for (num_complex::Complex<f64>, Data::<D>) {
    #[allow(private_interfaces)]
    #[inline]
    fn into_point(self, t: f64) -> Point<D> {
        Point::new_unchecked(t, [self.0.re, self.0.im], self.1.0)
    }
}


////////////////////////////////////////////////////////////////////////
//
// Defining a sampling with standard options & checks

/// Define a structure with standard fields, standard options, and a
/// function to generate it.
macro_rules! new_sampling_fn {
    // Function to init the struct.
    ($(#[$docfn: meta])*, $(#[$docfn_extra: meta])* $fun: ident -> $ft: ty,
     // The structure to hold the options (and other fields).
     $(#[$doc: meta])* $struct: ident,
     $wrap_f: ident
    ) => {
        impl<D> Sampling<D> {
            $(#[$docfn])*
            ///
            /// Panics if `a` or `b` is not finite.
            ///
            $(#[$docfn_extra])*
            #[must_use]
            #[allow(deprecated)]
            pub fn $fun<F, Y>(f: F, a: f64, b: f64) -> $struct<F, D>
            where F: FnMut(f64) -> Y,
                  Y: Img<D> {
                if !a.is_finite() {
                    panic!("curve_sampling::{}: a = {} must be finite",
                           stringify!($fun), a);
                }
                if !b.is_finite() {
                    panic!("curve_sampling::{}: b = {} must be finite",
                           stringify!($fun), b);
                }
                $struct { f,
                          a, b,  // Order of `a`, `b` reflect orientation
                          n: 100,
                          viewport: None,
                          init: vec![],
                          init_pt: vec![],
                }
            }
        }

        $(#[$doc])*
        pub struct $struct<F, D> {
            f: F,  a: f64,  b: f64,
            n: usize,
            viewport: Option<BoundingBox>,
            init: Vec<f64>,
            init_pt: Vec<Point<D>>,
        }

        #[allow(deprecated)]
        impl<F, Y, D> $struct<F, D>
        where F: FnMut(f64) -> Y,
              Y: Img<D>,
        {
            /// Set the maximum number of evaluations of the function
            /// to build the sampling.  Panic if `n < 2`.  An initial
            /// exploration of the function always uses at least 10
            /// evaluations so, for `n < 10`, this is the number of
            /// evaluations performed.
            pub fn n(mut self, n: usize) -> Self {
                if n < 2 {
                    panic!("curve_sampling: n = {} must at least be 2", n)
                }
                self.n = n;
                self
            }

            /// Set the zone of interest for the sampling.  Segments
            /// that end up outside this box will not be refined.
            pub fn viewport(mut self, vp: BoundingBox) -> Self {
                self.viewport = Some(vp);
                self
            }

            /// Add initial values of `t` such that `f(t)` (see [`
            #[doc = stringify!(Sampling::$fun)]
            /// `]) must be included into the sampling in addition to
            /// the `n` evaluations.  Only the values between `a` and
            /// `b` are taken into account (other values are ignored).
            pub fn init<'a, I>(mut self, ts: I) -> Self
            where I: IntoIterator<Item = &'a f64> {
                for &t in ts {
                    if self.a <= t && t <= self.b { // ⟹ t is finite
                        self.init.push(t);
                    }
                }
                self
            }

            /// Add initial points `(t, f(t))` to include into the
            /// sampling in addition to the `n` evaluations.  This
            /// allows you to use previous evaluations of `f`.  Only
            /// the couples with first coordinate `t` between `a` and
            /// `b` (see [`
            #[doc = stringify!(Sampling::$fun)]
            /// `]) are considered (other values are ignored).
            pub fn init_pt<I>(mut self, pts: I) -> Self
            where I: IntoIterator<Item = (f64, Y)>,
            {
                for (t, y) in pts {
                    if self.a <= t && t <= self.b { // ⟹ p.0 = t is finite
                        self.init_pt.push(y.into_point(t));
                    }
                }
                self
            }

            /// Evaluate the function at all initial values that where
            /// provided through [`Self::init`] and add the results to
            /// [`self.init_pt`].  Moreover sort all points by
            /// increasing or decreasing order of `Point`s `t`
            /// according to the order of `self.a` and `self.b`.
            fn eval_init_and_sort(&mut self) {
                self.init_pt.reserve(self.init.len());
                // `t` ∈ \[`a`, `b`\] already checked by [`init`] and
                // [`init_pt`].
                for &t in &self.init {
                    self.init_pt.push((self.f)(t).into_point(t))
                }
                if self.a <= self.b {
                    self.init_pt.sort_unstable_by(Point::cmp_t_incr);
                } else {
                    self.init_pt.sort_unstable_by(Point::cmp_t_decr);
                }
            }
        }
    }
}

////////////////////////////////////////////////////////////////////////
//
// Uniform sampling

new_sampling_fn!(
    /// Create a sampling for the graph of `f` on the interval
    /// \[`a`, `b`\] with evenly spaced values of the argument.
    ///
    /// See [`Sampling::fun`] for more documentation on the possible
    /// return values of `f`.
    ,
    /// # Example
    ///
    /// ```
    /// use std::{fs::File, io::BufWriter};
    /// use curve_sampling::Sampling;
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// let s = Sampling::uniform(|x| x.sin(), 0., 4.).build();
    /// s.write(&mut BufWriter::new(File::create("target/uniform.dat")?))?;
    /// # Ok(()) }
    /// ```
    uniform -> f64,
    /// Options for uniform sampling graphs of function ℝ → ℝ.
    /// See [`Sampling::uniform`].
    Uniform,
    FnPoint);

impl<F, Y, D> Uniform<F, D>
where F: FnMut(f64) -> Y,
      Y: Img<D> {
    /// Return a uniform sampling of the function.
    pub fn build(mut self) -> Sampling<D> {
        if self.a == self.b {
            // `a` is finite by previous tests
            let p = (self.f)(self.a).into_point(self.a);
            if p.is_valid() {
                return Sampling::singleton(p);
            } else {
                return Sampling::empty()
            }
        }
        self.eval_init_and_sort();
        let dt = (self.b - self.a) / (self.n - 1) as f64;
        let uniform = (0 .. self.n).into_iter().map(|i| {
            let t = self.a + i as f64 * dt;
            (self.f)(t).into_point(t)
        });
        let cmp = if self.a <= self.b {
            Point::cmp_t_incr
        } else {
            Point::cmp_t_decr
        };
        let points = MergeSortedIterator::new(
            self.init_pt.into_iter(),
            uniform,
            cmp);
        Sampling::from_point_iterator(points)
    }
}

////////////////////////////////////////////////////////////////////////
//
// Cost

mod cost {
    use super::{Point, Coord, Lengths};

    // The cost of a point is a measure of the curvature at this
    // point.  This requires segments before and after the point.  In
    // case the point is a cut, or first, or last, it has a cost of 0.
    // If it is an endpoint of a segment with the other point a cut,
    // the cost is set to [`HANGING_NODE`] because the segment with
    // the invalid point needs to be cut of too long to better
    // determine the boundary.  Only the absolute value of the cost
    // matters for the priority queue (see `segment`), its sign must
    // reflect the orientation of the angle.
    //
    // The cost of a point is apportioned to the segments of which it is
    // an endpoint according to their relative lengths.  More precisely,
    // the cost c of a point p is distributed on the segments s1 and s2
    // (of respective lengths l1 and l2) it is an endpoint of as
    //
    //   c * l1/(l1+l2) for s1 and c * l2/(l1+l2) for s2.
    //
    // In order to be able to update the cost of s1 without accessing
    // s2, p.cost holds c/(l1+l2).

    /// Cost for new "hanging" nodes — nodes created splitting a
    /// segment with an invalid endpoint.  Note that this cost will be
    /// multiplied by a function of `dt` in `segment` so it must be
    /// set high enough to ensure proper resolution of the endpoints
    /// of the domain.
    pub const HANGING_NODE: f64 = 5e5;

    /// Set the cost of the middle point `pm`.  Assumes `pm` is a
    /// valid point.
    #[inline]
    pub(crate) fn set_middle<D>(
        p0: Coord, pm: &mut Point<D>, p1: Coord, len: Lengths
    ) {
        let [x0, y0] = p0.xy;
        let [xm, ym] = pm.xy;
        let [x1, y1] = p1.xy;
        let dx0m = (x0 - xm) / len.x;
        let dy0m = (y0 - ym) / len.y;
        let dx1m = (x1 - xm) / len.x;
        let dy1m = (y1 - ym) / len.y;
        let len0m = dx0m.hypot(dy0m);
        let len1m = dx1m.hypot(dy1m);
        if len0m == 0. || len1m == 0. {
            pm.cost = 0.; // Do not subdivide
        } else {
            let dx = - dx0m * dx1m - dy0m * dy1m;
            let dy = dy0m * dx1m - dx0m * dy1m;
            pm.cost = dy.atan2(dx); // ∈ [-π, π]
            debug_assert!(!pm.cost.is_nan());
        }
    }

    /// Compute the cost of a segment according to the costs of its
    /// endpoints unless `in_vp` is false in which case -∞ is returned.
    #[inline]
    pub(crate) fn segment_vp(
        p0: Coord,
        p1: Coord,
        len: Lengths,
        in_vp: bool
    ) -> f64 {
        if ! in_vp { return f64::NEG_INFINITY }
        segment(p0, p1, len)
    }

    /// Compute the cost of a segment according to the costs of its
    /// endpoints.
    #[inline]
    pub(crate) fn segment(p0: Coord, p1: Coord, len: Lengths) -> f64 {
        let dt = (p1.t - p0.t).abs() / len.t; // ∈ [0,1]
        debug_assert!((0. ..=1.).contains(&dt), "dt = {dt} ∉ [0,1]");
        // Put less efforts when `dt` is small.  For functions, the
        // Y-variation may be large but, if it happens for a small range
        // of `t`, there is no point in adding indistinguishable details.
        let [x0, y0] = p0.xy;
        let [x1, y1] = p1.xy;
        let dx = ((x1 - x0) / len.x).abs();
        let dy = ((y1 - y0) / len.y).abs();
        let mut cost = p0.cost.abs() + p1.cost.abs(); // ≥ 0
        if p0.cost * p1.cost < 0. {
            // zigzag are bad on a large scale but less important on a
            // small scale.
            if dx <= 0.01 && dy <= 0.01 { cost *= 0.5 }
            else if !(dx <= 0.05 && dy <= 0.05) { cost *= 8. }
        }
        if dt >= 0.8 { cost }
        else {
            let dt = dt / 0.8;
            dt * dt * (6. + (-8. + 3. * dt) * dt) * cost
        }
    }

}

/// Compute the cost of the segment \[`p0`, `p1`\] (taking into
/// account `in_vp`) and push it to the queue `pq`.  `p0` is pointed
/// to by a list-witness so one can update it from the PQ elements.
/// `p0` is updated with the PQ-witness.
fn push_segment<D>(
    pq: &mut PQ<D>,
    p0: &Point<D>, w0: &list::Witness<Point<D>>, p1: Coord,
    len: Lengths, in_vp: bool
) {
    // FIXME: do we really want to push the segment when `!in_vp`?
    // In not all points are in the queue, one must beware of the
    // validity of witnesses though.
    let cost_segment = cost::segment_vp(p0.coord(), p1, len, in_vp);
    // The segment is referred to by its first point.
    let w = pq.push(cost_segment, w0.clone());
    *p0.witness.borrow_mut() = Some(w);
}

/// With the (new) costs in points `p0` and `p1`, update the position
/// of the segment \[`p0`, `p1`\] in the priority queue of `self`.
///
/// # Safety
/// `p0` must be in `pq`, otherwise it is UB.
unsafe fn update_segment<D>(
    pq: &mut PQ<D>,
    p0: &Point<D>,
    p1: Coord,
    len: Lengths
) {
    match &*p0.witness.borrow() {
        Some(w) => {
            let priority = cost::segment(p0.coord(), p1, len);
            pq.update_priority(w, priority)
        }
        None => panic!("Sampling::update_segment: unset witness"),
    }
}

////////////////////////////////////////////////////////////////////////
//
// Function sampling

/// Update the cost of all points in the sampling and add segments
/// to the priority queue.
fn compute<D>(s: &mut Sampling<D>, in_vp: impl Fn(&Point<D>) -> bool) {
    if let Some(len) = s.lengths() {
        s.path.first_mut().unwrap().cost = 0.;
        s.path.last_mut().unwrap().cost = 0.;
        let segments = unsafe { s.path.iter_segments_mut() };
        for (p0, w0, pm, p1) in segments {
            let p0_is_valid = p0.is_valid();
            let p0_in_vp = p0_is_valid && in_vp(p0);
            let pm_is_valid = pm.is_valid();
            let pm_in_vp;
            if let Some(p1) = p1 {
                if pm_is_valid {
                    pm_in_vp = in_vp(pm);
                    if p0_is_valid && p1.is_valid() {
                        let c0 = p0.coord();
                        let c1 = p1.coord();
                        cost::set_middle(c0, pm, c1, len);
                    } else {
                        pm.cost = cost::HANGING_NODE;
                    }
                } else { // pm is the location of a cut
                    pm_in_vp = false;
                    pm.cost = 0.;
                }
                if p0_is_valid || pm_is_valid {
                    // Add segment [p0, pm] to the PQ and set `p0` witness.
                    push_segment(&mut s.pq, p0, &w0, pm.coord(), len,
                        p0_in_vp || pm_in_vp);
                }
            } else { // Last segment (no `p1`)
                if p0_is_valid || pm_is_valid {
                    let mut vp = p0_in_vp;
                    if pm.is_valid() { vp = vp || in_vp(pm) };
                    push_segment(&mut s.pq, p0, &w0, pm.coord(), len, vp);
                }
            }
        }
    }
}

fn refine_gen<D>(
    s: &mut Sampling<D>, n: usize,
    mut f: impl FnMut(f64) -> Point<D>,
    in_vp: impl Fn(&Point<D>) -> bool
) {
    let len = match s.lengths() {
        Some(lengths) => lengths,
        None => return };
    macro_rules! r { ($x: ident) => { unsafe { s.path.get(&$x) } } }
    macro_rules! m { ($x: ident) => { unsafe { s.path.get_mut(&$x) } } }
    macro_rules! prev { ($x: ident) => { unsafe { s.path.prev(& $x) } } }
    macro_rules! next { ($x: ident) => { unsafe { s.path.next(& $x) } } }
    for _ in 0 .. n {
        let mut p0: list::Witness<Point<D>> = match s.pq.pop() {
            None => break,
            Some(p) => p };
        *r!(p0).witness.borrow_mut() = None; // PQ element it points to just popped.
        let p1 = next!(p0).unwrap();
        // Refine the segment [p0, p1] inserting a middle point `pm`.
        let t = (r!(p0).t + r!(p1).t) * 0.5;
        let mut pm = f(t);
        if r!(p0).is_valid() {
            if r!(p1).is_valid() {
                let pm = unsafe { s.path.insert_after(&mut p0, pm) };
                let mut pm_in_vp = false;
                if r!(pm).is_valid() {
                    pm_in_vp = in_vp(r!(pm));
                    let c0 = r!(p0).coord();
                    let c1 = r!(p1).coord();
                    cost::set_middle(c0, m!(pm), c1, len);
                    let cm = r!(pm).coord();
                    if let Some(p_1) = prev!(p0) && r!(p_1).is_valid() {
                        let c_1 = r!(p_1).coord();
                        cost::set_middle(c_1, m!(p0), cm, len);
                        let p_1 = r!(p_1);
                        let p0 = r!(p0).coord();
                        unsafe { update_segment(&mut s.pq, p_1, p0, len) }
                    }
                    if let Some(p2) = next!(p1) && r!(p2).is_valid() {
                        let c2 = r!(p2).coord();
                        cost::set_middle(cm, m!(p1), c2, len);
                        let p1 = r!(p1);
                        unsafe { update_segment(&mut s.pq, p1, c2, len) }
                    }
                } else { // `pm` invalid ⟹ cut between `p0` and `p1`
                    m!(p0).cost = cost::HANGING_NODE;
                    m!(pm).cost = 0.;
                    m!(p1).cost = cost::HANGING_NODE;
                    if let Some(p_1) = prev!(p0) {
                        let p_1 = r!(p_1);
                        let p0 = r!(p0).coord();
                        unsafe { update_segment(&mut s.pq, p_1, p0, len) }
                    }
                    if let Some(p2) = next!(p1) {
                        let p1 = r!(p1);
                        let p2 = r!(p2).coord();
                        unsafe { update_segment(&mut s.pq, p1, p2, len) }
                    }
                }
                let vp = pm_in_vp || in_vp(r!(p0));
                let cm = r!(pm).coord();
                push_segment(&mut s.pq, r!(p0), &p0, cm, len, vp);
                let vp = pm_in_vp || in_vp(r!(p1));
                let c1 = r!(p1).coord();
                push_segment(&mut s.pq, r!(pm), &pm, c1, len, vp);
            } else { // `p0` valid, `p1` invalid (i.e. is a cut)
                // Thus `p0` is a hanging node.
                if pm.is_valid() {
                    pm.cost = cost::HANGING_NODE;
                    let pm = unsafe { s.path.insert_after(&mut p0, pm) };
                    if let Some(p_1) = prev!(p0) && r!(p_1).is_valid() {
                        let c_1 = r!(p_1).coord();
                        let pm = r!(pm).coord();
                        cost::set_middle(c_1, m!(p0), pm, len);
                        let p_1 = r!(p_1);
                        let p0 = r!(p0).coord();
                        unsafe{ update_segment(&mut s.pq, p_1, p0, len) }
                    }
                    let pm_in_vp = in_vp(r!(pm));
                    let vp = pm_in_vp || in_vp(r!(p0));
                    let cm = r!(pm).coord();
                    let c1 = r!(p1).coord();
                    push_segment(&mut s.pq, r!(p0), &p0, cm, len, vp);
                    push_segment(&mut s.pq, r!(pm), &pm, c1, len, pm_in_vp)
                } else { // `pm` invalid
                    // Insert only \[`p0`, `pm`\] and forget
                    // \[`pm`, `p1`\].  The cost of `p0` stays
                    // `cost::HANGING_NODE`.  We can see this as
                    // reducing the uncertainty of the boundary in the
                    // segment \[`p0`, `p1`\].
                    pm.cost = 0.;
                    let pm = {
                        if r!(p1).witness.borrow().is_none() {
                            // `p1` is not part of a segment.  One can
                            // replace it by `pm`.
                            unsafe { *s.path.get_mut(&p1) = pm } ;
                            p1 // witness for `pm` now.
                        } else {
                            unsafe { s.path.insert_after(&mut p0, pm) }
                        }
                    };
                    let vp = in_vp(r!(p0));
                    let cm = r!(pm).coord();
                    push_segment(&mut s.pq, r!(p0), &p0, cm, len, vp)
                }
            }
        } else { // `p0` invalid (i.e., cut) ⟹ `p1` valid
            debug_assert!(r!(p1).is_valid());
            if pm.is_valid() {
                pm.cost = cost::HANGING_NODE;
                let pm = unsafe { s.path.insert_after(&mut p0, pm) };
                if let Some(p2) = next!(p1) && r!(p2).is_valid() {
                    let pm = r!(pm).coord();
                    let p2 = r!(p2).coord();
                    cost::set_middle(pm, m!(p1), p2, len);
                    let p1 = r!(p1);
                    unsafe{ update_segment(&mut s.pq, p1, p2, len) }
                }
                let pm_in_vp = in_vp(r!(pm));
                let cm = r!(pm).coord();
                let c1 = r!(p1).coord();
                push_segment(&mut s.pq, r!(p0), &p0, cm, len, pm_in_vp);
                let in_vp = pm_in_vp || in_vp(r!(p1));
                push_segment(&mut s.pq, r!(pm), &pm, c1, len, in_vp);
            } else { // `pm` invalid ⟹ drop segment \[`p0`, `pm`\].
                // Cost of `p1` stays `cost::HANGING_NODE`.
                pm.cost = 0.;
                let pm = {
                    if let Some(p_1) = prev!(p0) {
                        if r!(p_1).is_valid() {
                            unsafe { s.path.insert_after(&mut p0, pm) }
                        } else {
                            // `p_1` is the cut ending the previous segment.
                            unsafe { *s.path.get_mut(&p0) = pm };
                            p0
                        }
                    } else {
                        unsafe { s.path.insert_after(&mut p0, pm) }
                    }
                };
                let vp = in_vp(r!(p1));
                let c1 = r!(p1).coord();
                push_segment(&mut s.pq, r!(pm), &pm, c1, len, vp)
            }
        }
    }
}

struct AlmostUniformSampling<D, F>
where F: FnMut(f64) -> Point<D> {
    f: F,
    a: f64,
    b: f64,
    n: usize,
    i: usize,
    dt: f64,
    random: u32,
}

/// Return an iterator that yields evaluations of `f` at `n` almost
/// equispaced values in the interval \[`a`, `b`\].
///
/// Remark: The randomization seed only depends on the parameter `n`
/// for simple reproducibility.
fn almost_uniform_sampling<D, F: FnMut(f64) -> Point<D>>(
    f: F,
    a: f64, b: f64, n: usize
) -> AlmostUniformSampling<D, F> {
    let n = n.max(4);
    let dt = (b - a) / (n - 3) as f64;
    let random = (n as u32).wrapping_mul(1_000_000);
    AlmostUniformSampling { f, a, b, n, i: 0, dt, random }
}

impl<D, F: FnMut(f64) -> Point<D>> Iterator for AlmostUniformSampling<D, F> {
    type Item = Point<D>;

    fn next(&mut self) -> Option<Point<D>> {
        const NORMALIZE_01: f64 = 1. / u32::MAX as f64;

        if self.i >= self.n { return None }
        let f = &mut self.f;
        let v =
            if self.i == 0 { f(self.a) }
            else if self.i == 1 { f(self.a + 0.0625 * self.dt) }
            else if self.i == self.n - 2 { f(self.b - 0.0625 * self.dt) }
            else if self.i == self.n - 1 { f(self.b) }
            else { // i ∈ {2,...,n-3}
                // Pseudorandom number generator from the "Xorshift
                // RNGs" paper by George Marsaglia.  See
                // https://matklad.github.io/2023/01/04/on-random-numbers.html
                // and
                // https://github.com/rust-lang/rust/blob/1.55.0/library/core/src/slice/sort.rs#L559-L573
                self.random ^= self.random << 13;
                self.random ^= self.random >> 17;
                self.random ^= self.random << 5;
                let rand = self.random as f64 * NORMALIZE_01;
                let j = (self.i - 1) as f64 + rand * 0.125 - 0.0625;
                f(self.a + j * self.dt)
            };
        self.i += 1;
        Some(v)
    }
}

impl<D> Sampling<D> {
    /// Return a sampling from the initial list of `points`.
    ///
    /// # Preconditions
    /// - All `points` are in \[`a`, `b`\].
    /// - `points` is sorted in increasing order (if `a` ≤ `b`) or
    ///   decreasing otherwise.
    fn build(
        points: Vec<Point<D>>,
        mut f: impl FnMut(f64) -> Point<D>,
        a: f64, b: f64, n: usize,
        viewport: Option<BoundingBox>
    ) -> Self {
        if a == b {
            let p = f(a);
            if p.is_valid() {
                return Sampling::singleton(p);
            } else {
                return Sampling::empty()
            }
        }
        // Uniform sampling requires ≥ 4 points but actually makes no
        // sense with less than 10 points.  Thus, if `n < 10`, more
        // than `n` evaluations are performed and no refinement is done
        // (`n - n0` must not underflow).
        let n0 = (n / 10).max(10);
        let uniform = almost_uniform_sampling(&mut f, a, b, n0);
        let points = MergeSortedIterator::new(
            points.into_iter(),
            uniform,
            if a <= b { Point::cmp_t_incr } else { Point::cmp_t_decr });
        let mut s = Sampling::from_point_iterator(points);
        match viewport {
            Some(vp) => {
                let in_vp = |p: &Point<_>| vp.contains(p);
                compute(&mut s, in_vp);
                refine_gen(&mut s, n.saturating_sub(n0), &mut f, in_vp)
            }
            None => {
                compute(&mut s, |_| true);
                refine_gen(&mut s, n.saturating_sub(n0), &mut f, |_| true)
            }
        }
        s
    }
}


new_sampling_fn!(
    /// When `f` returns `f64` (resp. `[f64; 2]`), construct a
    /// sampling of the *graph* (resp. the *image*) of `f` on the
    /// interval \[`a`, `b`\] by evaluating `f` at `n` points (see
    /// [`Fun::n`]).  The function `f` may also provide additional
    /// data by returning `(y, Data(d))` of type `(f64, Data::<D>)` or
    /// `([x, y], Data(d))` of type `([f64; 2], Data::<D>)` (the type
    /// parameter `D` will be inferred from the type of the data `d`).
    ///
    /// If the feature `num-complex` is enabled, complex numbers may
    /// also be used.
    ,
    /// # Examples
    ///
    /// ```
    /// use std::{fs::File, io::BufWriter};
    /// use curve_sampling::Sampling;
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// let s = Sampling::fun(|x| x.sin(), 0., 4.).build();
    /// s.write(&mut BufWriter::new(File::create("target/fun1.dat")?))?;
    /// let s = Sampling::fun(|x| [x.sin(), x.cos()], 0., 4.).build();
    /// s.write(&mut BufWriter::new(File::create("target/fun2.dat")?))?;
    /// let s = Sampling::fun(|x| x.sin_cos(), 0., 4.).build();
    /// s.write(&mut BufWriter::new(File::create("target/fun3.dat")?))?;
    /// # Ok(()) }
    /// ```
    ///
    /// Sampling with associated data:
    ///
    /// ```
    /// use curve_sampling::{Sampling, Data};
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// fn f(x: f64) -> [f64; 3] {
    ///     [x, x*x, x*x*x]
    /// }
    /// let s = Sampling::fun(|x| {
    ///     let fx = f(x);
    ///     ([fx[0], fx[1]], Data(fx))
    /// }, 0., 1.).build();
    /// // ...
    /// # Ok(()) }
    /// ```
    ///
    #[cfg_attr(feature = "num-complex", doc = "
With the feature `num-complex`, complex-valued functions may be used directly.

```
use num_complex::Complex64 as C64;
use curve_sampling::{Sampling, Data};
# fn main() -> Result<(), Box<dyn std::error::Error>> {
fn f(t: f64) -> C64 {
    C64::i().powf(t)
}
let s = Sampling::fun(f, 0., 1.).build();
// ...
# Ok(()) }
```
")]
    fun -> f64,
    /// Options for sampling graphs of functions ℝ → ℝ and ℝ → ℝ².
    /// See [`Sampling::fun`].
    Fun,
    FnPoint);

impl<F, Y, D> Fun<F, D>
where F: FnMut(f64) -> Y,
      Y: Img<D> {
    /// Return the sampling.
    pub fn build(mut self) -> Sampling<D> {
        self.eval_init_and_sort();
        Sampling::build(
            self.init_pt,
            |x| (self.f)(x).into_point(x),
            self.a, self.b, self.n, self.viewport)
    }
}


new_sampling_fn!(
    #[deprecated(note = "Use Sampling::fun instead")]
    /// Create a sampling of the *image* of `f` on the interval
    /// \[`a`, `b`\] by evaluating `f` at `n` points.
    ,
    /// # Example
    ///
    /// ```
    /// use std::{fs::File, io::BufWriter};
    /// use curve_sampling::Sampling;
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// let s = Sampling::param(|t| [t.sin(), t.cos()], 0., 4.).build();
    /// s.write(&mut BufWriter::new(File::create("target/param.dat")?))?;
    /// # Ok(()) }
    /// ```
    param -> [f64; 2],
    #[deprecated(note = "See curve_sampling::Fun instead")]
    /// Options for sampling the image of functions ℝ → ℝ².
    /// See [`Sampling::param`].
    Param,
    ParamPoint);

#[allow(deprecated)]
impl<F, Y, D> Param<F, D>
where F: FnMut(f64) -> Y,
      Y: Img<D> {
    /// Return the sampling.
    pub fn build(mut self) -> Sampling<D> {
        self.eval_init_and_sort();
        Sampling::build(
            self.init_pt,
            |t| (self.f)(t).into_point(t),
            self.a, self.b, self.n, self.viewport)
    }
}


////////////////////////////////////////////////////////////////////////
//
// Output

/// LaTeX output, created by [`Sampling::latex`].
///
/// # Example
///
/// ```
/// use std::fs::File;
/// use curve_sampling::Sampling;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let s = Sampling::from_iter([[0., 0.], [1., 1.]]);
/// s.latex().write(&mut File::create("target/sampling.tex")?)?;
/// # Ok(()) }
/// ```
pub struct LaTeX<'a, D> {
    sampling: &'a Sampling<D>,
    n: usize,
    color: Option<RGB8>,
    arrow: Option<&'a str>,
    arrow_pos: Option<f64>,
}

impl<'a, D> LaTeX<'a, D> {
    #[inline]
    fn new(s: &'a Sampling<D>) -> Self {
        Self { sampling: s,  n: 20_000,  color: None,
               arrow: None,  arrow_pos: None }
    }

    /// Set the maximum number of points of a PGF path to `n`.  If it
    /// contains more than `n` points, the sampling curve is drawn as
    /// several PGF paths.  Default: 20_000.
    pub fn n(&mut self, n: usize) -> &mut Self {
        self.n = n;
        self
    }

    /// Set the color of the curve to `color`.  If not specified the
    /// active LaTeX color will be used.
    pub fn color(&mut self, color: RGB8) -> &mut Self {
        self.color = Some(color);
        self
    }

    /// Set the type of arrow to draw to `arrow`.
    /// See the documentation of `\pgfsetarrowsend` in the
    /// [TikZ manual](https://tikz.dev/base-actions#sec-103.3).
    pub fn arrow(&mut self, arrow: &'a str) -> &mut Self {
        self.arrow = Some(arrow);
        self
    }

    /// The position of the arrow as a percent of the curve length (in
    /// the interval \[0.,1.\]).  If [`LaTeX::arrow`] is specified but
    /// not this, it defaults to `0.5`.
    pub fn arrow_pos(&mut self, arrow_pos: f64) -> &mut Self {
        if ! arrow_pos.is_finite() {
            panic!("curve_sampling::LaTeX::arrow_pos: \
                    position must be finite");
        }
        self.arrow_pos = Some(arrow_pos.clamp(0., 1.));
        self
    }

    /// Write the sampling with lines segments.
    fn write_with_lines(&self, f: &mut impl Write) -> Result<(), io::Error> {
        let mut n = 0;
        let mut new_path = true;
        for [x,y] in self.sampling.iter() {
            if x.is_nan() {
                writeln!(f, "\\pgfusepath{{stroke}}")?;
                n = 0;
                new_path = true;
            } else {
                n += 1;
                if new_path {
                    writeln!(f, "\\pgfpathmoveto{{\\pgfpointxy\
                                 {{{x}}}{{{y}}}}}")?
                } else if n >= self.n {
                    write!(f, "\\pgfpathlineto{{\\pgfpointxy\
                               {{{x}}}{{{y}}}}}\n\
                               \\pgfusepath{{stroke}}\n\
                               \\pgfpathmoveto{{\\pgfpointxy\
                               {{{x}}}{{{y}}}}}\n")?;
                    n = 0;
                } else {
                    writeln!(f, "\\pgfpathlineto{{\\pgfpointxy\
                                 {{{x}}}{{{y}}}}}")?
                }
                new_path = false;
            }
        }
        Ok(())
    }

    /// Write the path, each continuous sub-path with an arrow.
    fn write_with_arrows(&self, f: &mut impl Write, arrow: &str,
                         arrow_pos: f64) -> Result<(), io::Error> {
        // Compute the length of all sub-paths.
        let mut lens = vec![];
        let mut prev_pt: Option<[f64; 2]> = None;
        let mut cur_len = 0.;
        for p @ [x, y] in self.sampling.iter() {
            if x.is_nan() {
                lens.push(arrow_pos * cur_len);
                prev_pt = None;
                cur_len = 0.;
            } else {
                if let Some([x0, y0]) = prev_pt {
                    cur_len += (x - x0).hypot(y - y0);
                }
                prev_pt = Some(p);
            }
        }
        lens.push(arrow_pos * cur_len);
        if lens.is_empty() { return Ok(()) }
        let mut lens = lens.iter();
        let mut rem_len = *lens.next().unwrap(); // remaining before arrow
        prev_pt = None;
        let mut n = 0;
        for p @ [x, y] in self.sampling.iter() {
            if x.is_nan() {
                writeln!(f, "\\pgfusepath{{stroke}}")?;
                rem_len = *lens.next().unwrap();
                prev_pt = None;
            } else {
                n += 1;
                if let Some([x0, y0]) = prev_pt {
                    let dx = x - x0;
                    let dy = y - y0;
                    let l = dx.hypot(dy);
                    if rem_len <= l {
                        writeln!(f, "\\pgfusepath{{stroke}}")?;
                        // Drawing a long path with an arrow specified is
                        // extremely expensive.  Just draw the current segment.
                        let pct = rem_len / l;
                        if pct <= 1e-14 {
                            write!(f, "\\pgfsetarrowsstart{{{arrow}}}\n\
                                       \\pgfpathmoveto{{\\pgfpointxy
                                       {{{x0}}}{{{y0}}}}}\n\
                                       \\pgfpathlineto{{\\pgfpointxy\
                                       {{{x}}}{{{y}}}}}\n\
                                       \\pgfusepath{{stroke}}\n\
                                       \\pgfsetarrowsend{{}}\n\
                                       \\pgfpathmoveto{{\\pgfpointxy\
                                       {{{x}}}{{{y}}}}}\n")?;
                        } else {
                            let xm = x0 + pct * dx;
                            let ym = y0 + pct * dy;
                            write!(f, "\\pgfsetarrowsend{{{arrow}}}\n\
                                       \\pgfpathmoveto{{\\pgfpointxy\
                                       {{{x0}}}{{{y0}}}}}\n\
                                       \\pgfpathlineto{{\\pgfpointxy\
                                       {{{xm}}}{{{ym}}}}}\n\
                                       \\pgfusepath{{stroke}}\n\
                                       \\pgfsetarrowsend{{}}\n\
                                       \\pgfpathmoveto{{\\pgfpointxy\
                                       {{{xm}}}{{{ym}}}}}\n\
                                       \\pgfpathlineto{{\\pgfpointxy\
                                       {{{x}}}{{{y}}}}}\n")?;
                            n = 2;
                        }
                        rem_len = f64::INFINITY; // No more arrow for this sub-path
                    } else if n >= self.n {
                        write!(f, "\\pgfpathlineto{{\\pgfpointxy\
                                   {{{x}}}{{{y}}}}}\n\
                                   \\pgfusepath{{stroke}}\n\
                                   \\pgfpathmoveto{{\\pgfpointxy\
                                   {{{x}}}{{{y}}}}}\n")?;
                        n = 0;
                    } else {
                        writeln!(f, "\\pgfpathlineto{{\\pgfpointxy\
                                     {{{x}}}{{{y}}}}}")?
                    }
                    rem_len -= l;
                } else {
                    // No previous point.  New sub-path.
                    writeln!(f, "\\pgfpathmoveto{{\\pgfpointxy\
                                 {{{x}}}{{{y}}}}}")?
                }
                prev_pt = Some(p);
            }
        }
        Ok(())
    }

    /// Write the sampling to the formatter as PGF/TikZ commands.
    pub fn write(&self, f: &mut impl Write) -> Result<(), io::Error> {
        writeln!(f, "% Written by the Rust curve-sampling crate.")?;
        writeln!(f, "\\begin{{pgfscope}}")?;
        if let Some(RGB8 {r, g, b}) = self.color {
            write!(f, "\\definecolor{{RustCurveSamplingColor}}{{RGB}}\
                       {{{r},{g},{b}}}\n\
                       \\pgfsetstrokecolor{{RustCurveSamplingColor}}\n")?
        }
        match (self.arrow, self.arrow_pos) {
            (None, None) => self.write_with_lines(f)?,
            (Some(arrow), None) =>
                self.write_with_arrows(f, arrow, 0.5)?,
            (None, Some(arrow_pos)) =>
                self.write_with_arrows(f, ">",arrow_pos)?,
            (Some(arrow), Some(arrow_pos)) =>
                self.write_with_arrows(f, arrow, arrow_pos)?,
        }
        write!(f, "\\pgfusepath{{stroke}}\n\\end{{pgfscope}}\n")
    }
}

/// # Output
impl<D> Sampling<D> {
    /// Write the sampling `self` using PGF/TikZ commands.
    pub fn latex(&self) -> LaTeX<'_, D> { LaTeX::new(self) }

    /// Write the sampling to `f` in a tabular form: each point is
    /// written as "x y" on a single line (in scientific notation).
    /// If the path is interrupted, a blank line is printed.  This
    /// format is compatible with Gnuplot.
    pub fn write(&self, f: &mut impl Write) -> Result<(), io::Error> {
        for [x, y] in self.iter() {
            if x.is_nan() {
                writeln!(f)?
            } else {
                writeln!(f, "{x:e} {y:e}")?
            }
        }
        Ok(())
    }
}

impl<D> Display for Sampling<D> {
    /// Display the sampling in a tabular form: each point is written
    /// as "x y" on a single line (in scientific notation).  If the
    /// path is interrupted, a blank line is printed.  This format is
    /// compatible with Gnuplot.
    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), fmt::Error> {
        for [x, y] in self.iter() {
            if x.is_nan() {
                writeln!(f)?
            } else {
                writeln!(f, "{x:e} {y:e}")?
            }
        }
        Ok(())
    }
}

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

#[cfg(test)]
mod tests {
    use std::{error::Error,
              fs::File,
              io::Write, path::Path};
    use crate::{BoundingBox, Data, MergeSortedIterator, NoData, Point, Sampling};

    #[test]
    fn merge_sorted_iterator() {
        let i1 = [1, 3, 5, 7].into_iter();
        let i2 = [1, 2, 5, 8].into_iter();
        let i = MergeSortedIterator::new(i1, i2, |i1, i2| i1.cmp(i2));
        assert_eq!(i.collect::<Vec<_>>(), [1, 2, 3, 5, 7, 8])
    }

    type R<T> = Result<T, Box<dyn Error>>;

    fn xy_of_sampling<D>(s: &Sampling<D>) -> Vec<Option<(f64, f64)>> {
        s.path.iter().map(|p| {
            if p.is_valid() { Some((p.xy[0], p.xy[1])) } else { None }
        })
            .collect()
    }

    #[test]
    fn clone_sampling() {
        let s = Sampling::from_iter([[0.,0.], [1.,1.]]);
        let xy0: Vec<_> = s.iter().collect();
        let xy1: Vec<_> = s.clone().iter().collect();
        assert_eq!(xy0, xy1)
    }

    #[test]
    fn io() {
        let s = Sampling::from_iter([[0.,0.], [1.,1.]]);
        assert_eq!(xy_of_sampling(&s), vec![Some((0.,0.)), Some((1.,1.))]);
        let s = Sampling::from_iter([[0.,0.], [1.,1.], [f64::NAN, 1.],
                                     [2.,2.]]);
        assert_eq!(xy_of_sampling(&s),
                   vec![Some((0.,0.)), Some((1.,1.)), None, Some((2.,2.))]);
    }

    #[test]
    fn bounding_box_singleton() {
        let s = Sampling::singleton(
            Point::new_unchecked(0., [1., 2.], NoData::new()));
        let bb = BoundingBox {xmin: 1., xmax: 1., ymin: 2., ymax: 2.};
        assert_eq!(s.bounding_box(), bb);
    }

        #[test]
    fn singleton_has_empty_bounding_box() {
        let s = Sampling::singleton(
            Point::new_unchecked(0., [1., 2.], NoData::new()));
        assert!(s.bounding_box().is_empty());
    }

    fn test_clip(bb: BoundingBox,
                 path: Vec<[f64;2]>,
                 expected: Vec<Option<(f64,f64)>>) {
        let s = Sampling::from_iter(path).clip(bb);
        assert_eq!(xy_of_sampling(&s), expected);
    }

    #[test]
    fn clip_base () {
        let bb = BoundingBox { xmin: 0.,  xmax: 3., ymin: 0.,  ymax: 2.};
        test_clip(bb, vec![[0.,1.], [2.,3.]],
                  vec![Some((0.,1.)), Some((1.,2.))]);
        test_clip(bb,
                  vec![[-1.,0.], [2.,3.], [4.,1.]],
                  vec![Some((0.,1.)), Some((1.,2.)), None, Some((3., 2.))]);
        test_clip(bb,
                  vec![[0.,1.], [2.,3.], [4.,1.], [2.,1.], [2., -1.],
                       [0., -1.], [0., 1.] ],
                  vec![Some((0.,1.)), Some((1.,2.)), None,
                       Some((3., 2.)), None,
                       Some((3., 1.)), Some((2., 1.)), Some((2., 0.)), None,
                       Some((0., 0.)), Some((0., 1.))]);
    }

    #[test]
    fn clip_empty() {
        let bb = BoundingBox {xmin: 0., xmax: 1., ymin: 0., ymax: 1.};
        let path = xy_of_sampling(&Sampling::<()>::empty().clip(bb));
        assert_eq!(path, vec![]);
    }

    #[test]
    fn clip_double_cut() {
        let bb = BoundingBox { xmin: 0.,  xmax: 4., ymin: 0.,  ymax: 2.};
        test_clip(bb,
                  vec![[1., 2.], [2.,3.], [5.,0.], [-1., 0.] ],
                  vec![Some((1., 2.)), None,
                       Some((3., 2.)), Some((4., 1.)), None,
                       Some((4., 0.)), Some((0., 0.)) ]);
    }

    #[test]
    fn clip_almost_horiz() {
        let bb = BoundingBox { xmin: 0.,  xmax: 2., ymin: -1.,  ymax: 1.};
        test_clip(bb,
                  vec![[1., 1e-100], [3., -1e-100] ],
                  vec![Some((1., 1e-100)), Some((2., 0.))]);
    }

    #[test]
    fn clip_touches_1pt_cut() {
        let bb = BoundingBox {xmin: 0., xmax: 2., ymax: 0., ymin: -1.};
        let cut = [f64::NAN, f64::NAN];
        test_clip(bb,
                  vec![[0.,1.], cut, [1., 0.], cut, cut, [2., 1.]],
                  vec![Some((1., 0.))])
    }

    #[test]
    fn clip_final_cut() {
        let bb = BoundingBox {xmin: 0., xmax: 1., ymin: 0., ymax: 2.};
        test_clip(bb,
                  vec![[0., 0.], [2., 2.]],
                  vec![Some((0., 0.)), Some((1., 1.))])
    }

    #[test]
    fn small_n() {
        // The initial exploration uses at least 10 points, even if
        // fewer evaluations were requested.  The number of refinement
        // steps `n - n0` used to underflow in that case: a panic in
        // debug mode, an endless loop in release mode.
        for n in 2 ..= 12 {
            let mut evals = 0;
            let s = Sampling::fun(|x: f64| { evals += 1; x.sin() }, 0., 3.)
                .n(n).build();
            assert_eq!(evals, n.max(10), "n = {n}");
            assert_eq!(s.len(), evals);
        }
    }

    #[test]
    fn uniform1() {
        let s = Sampling::uniform(|x| x, 0., 4.).n(3)
            .init(&[1.]).init_pt([(3., 0.)]).build();
        assert_eq!(xy_of_sampling(&s),
                   vec![Some((0.,0.)), Some((1.,1.)), Some((2.,2.)),
                        Some((3., 0.)), Some((4.,4.))]);
    }

    #[test]
    fn uniform2() {
        let s = Sampling::uniform(|x| (4. - x).sqrt(), 0., 6.).n(4).build();
        assert_eq!(xy_of_sampling(&s),
                   vec![Some((0.,2.)), Some((2., 2f64.sqrt())),
                        Some((4., 0.)), None]);
    }

    #[test]
    fn iter_data() {
        let s = Sampling::uniform(|x| (x, Data(x as i32)), 0., 4.).n(3)
            .init(&[1.]).init_pt([(3., (0., Data(-1)))]).build();
        let expected = vec![
            ([0.,0.], &0), ([1.,1.], &1), ([2.,2.], &2), ([3., 0.], &-1),
            ([4.,4.], &4)];
        let (_, size_hint) = s.iter_data().size_hint();
        assert_eq!(size_hint, Some(expected.len()));
        for (i, d) in s.iter_data().enumerate() {
            assert_eq!(d, expected[i]);
        }
    }

    #[test]
    fn into_iter_data() {
        let s = Sampling::uniform(|x| (x, Data(x as i32)), 0., 4.).n(3)
            .init(&[1.]).init_pt([(3., (0., Data(-1)))]).build();
        let expected = vec![
            ([0.,0.], 0), ([1.,1.], 1), ([2.,2.], 2), ([3., 0.], -1),
            ([4.,4.], 4)];
        let (_, size_hint) = s.iter_data().size_hint();
        assert_eq!(size_hint, Some(expected.len()));
        for (i, d) in s.into_iter_data().enumerate() {
            assert_eq!(d, expected[i]);
        }
    }

    /// In order the judge the quality of the sampling, we save it
    /// with the internal cost data.
    fn write_with_point_costs<D>(
        s: &Sampling<D>,
        fname: impl AsRef<Path>
    ) -> R<()> {
        let mut fh = File::create(fname)?;
        for p in s.path.iter() {
            if p.is_valid() {
                let [x, y] = p.xy;
                writeln!(fh, "{x} {y} {}", p.cost)?;
            } else {
                writeln!(fh)?;
            }
        }
        Ok(())
    }

    fn write_segments<D>(mut s: Sampling<D>, fname: impl AsRef<Path>) -> R<()> {
        let mut fh = File::create(fname)?;
        let mut seg = vec![];
        loop {
            let priority = s.pq.max_priority();
            if let Some(p0) = s.pq.pop() {
                let p1 = unsafe { s.path.next(&p0).unwrap() };
                let p1 = unsafe { s.path.get(&p1) };
                let p0 = unsafe { s.path.get(&p0) };
                let tm = (p0.t + p1.t) / 2.;
                seg.push((tm, p0.coord(), p1.coord(), priority))
            } else {
                break;
            }
        }
        seg.sort_by(|(t1,_,_,_), (t2,_,_,_)| t1.partial_cmp(t2).unwrap());
        for (tm, p0, p1, priority) in seg {
            writeln!(fh, "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}",
                tm, p0.t, p0.xy[0], p0.xy[1],
                p1.t, p1.xy[0], p1.xy[1], priority)?;
        }
        Ok(())
    }

    #[derive(Clone)]
    struct Plot {
        xmin: f64,  xmax: f64,
        ymin: f64,  ymax: f64,
        n: usize,
        init: Vec<f64>,
    }

    impl Plot {
        /// Return the Gnuplot instructions to plot the data.
        fn plot<F>(&self, f: F, title: &str) -> R<String>
        where F: FnMut(f64) -> f64 {
            let vp = BoundingBox { xmin: self.xmin, xmax: self.xmax,
                                   ymin: self.ymin, ymax: self.ymax };
            let s = Sampling::fun(f, self.xmin, self.xmax)
                .n(self.n).init(self.init.iter()).viewport(vp).build();
            static mut NDAT: usize = 0;
            let ndat = unsafe { NDAT += 1;  NDAT };
            println!("Horror {ndat}");
            s.pq.print_stats();
            let dir = Path::new("target");
            let fname = format!("horror{}.dat", ndat);
            s.write(&mut File::create(dir.join(&fname))?)?;
            let fname_p = format!("horror{}_p.dat", ndat);
            write_with_point_costs(&s, dir.join(&fname_p))?;
            let fname_s = format!("horror{}_s.dat", ndat);
            write_segments(s, dir.join(&fname_s))?;
        Ok(format!(
            "unset title\n\
             unset y2tics\n\
             plot [{}:{}] \"{}\" with l lt 5 title \"{}\", \
             \"{}\" with p lt 3 pt 6 ps 0.2 title \"n={}\"\n\
             set title \"#{ndat} Restricted to viewport [{}:{}]×[{}:{}]\"\n\
             set y2tics\n\
             set y2range [-1e-6:]\n\
             plot [{}:{}] [{}:{}] \"{}\" with l lt 5 title \"{}\", \
             \"{}\" with p lt 3 pt 7 ps 0.2 title \"n={}\", \
             \"{}\" using 1:3  with lp ps 0.2 lt rgb \"#760b0b\" \
             title \"cost points\", \
             \"{}\" using 1:8 with lp ps 0.2 lt rgb \"#909090\" \
             axes x1y2 title \"cost segments\"\n",
            self.xmin, self.xmax, &fname, title, &fname, self.n,
            self.xmin, self.xmax, self.ymin, self.ymax,
            self.xmin, self.xmax, self.ymin, self.ymax, &fname, title,
            &fname, self.n, &fname_p, &fname_s))
        }
    }

    #[test]
    #[cfg_attr(miri, ignore)] // very slow under Miri
    fn horror() -> R<()> {
        let d = Plot {
            xmin: -5., xmax: 5., ymin: -5., ymax: 5., n: 100, init: vec![] };
        macro_rules! p {
            ($($id:ident $l:tt $e:expr),*) => {
                Plot{ $($id: $e,)* ..d.clone() } };
        }
        let s = [
            p!(n: 10).plot(|_| 2., "x ↦ 2")?,
            p!().plot(|x| x, "x ↦ x")?,
            p!().plot(|x| 5. * x, "x ↦ 5x")?,
            p!().plot(|x| 1e6 * x, "x ↦ 10⁶ x")?, // high slope
            p!().plot(|x| 1e50 * x, "x ↦ 10⁵⁰ x")?, // high slope
            p!().plot(|x| 1. / x, "x ↦ 1/x")?, // check singularity
            p!(xmin: 0., xmax: 5., ymax: 100.).plot(
                |x| 1. / x, "x ↦ 1/x")?, // singularity at starting point
            p!(xmin: -0.4, xmax: 2., ymin: 0., ymax: 1.6, n: 50).plot(
                |x| x.sqrt(), "x ↦ √x")?,
            // Test cuts also to the right:
            p!(xmin: -2., xmax: 1., ymin: 0., ymax: 1.6, n: 50).plot(
                |x| (-x).sqrt(), "x ↦ √(-x)")?,
            p!(n: 200).plot(|x| x.tan(), "tan")?,
            p!().plot(|x| 1. / x.abs(), "x ↦ 1/|x|")?,
            p!(xmin: -6., xmax: 6., ymin: -2., ymax: 2.).plot(
                |x| (1. + x.cos().sin()).ln(), "1 + sin(cos x)")?,
            p!(xmin: 0., xmax: 6.28, ymin: -1.5, ymax: 1.5, n: 400).plot(
                |x| x.powi(3).sin() + x.powi(3).cos(), "sin x³ + cos x³")?,
            p!(xmin: -5., xmax:200., ymin: -1., ymax: 1., n: 400).plot(
                |x| x.sin(), "sin")?,
            // Examples from R. Avitzur, O. Bachmann, N. Kajler, "From
            // Honest to Intelligent Plotting", proceedings of ISSAC'
            // 95, pages 32-41, July 1995.
            p!(xmin: -4., xmax: 4., ymin: -1., ymax: 1.).plot(
                |x| (300. * x).sin(), "sin(300 x)")?,
            p!(xmin: -4., xmax: 4., ymin: -1., ymax: 1., n: 1000).plot(
                |x| (300. * x).sin(), "sin(300 x)")?,
            p!(xmin: -2., xmax: 2., ymin: 0., ymax: 3.).plot(
                |x| 1. + x * x + 0.0125 * (1. - 3. * (x - 1.)).abs().ln(),
                "1 + x² + 0.0125 ln|1 - 3(x-1)|")?,
            p!(xmin: -2., xmax: 2., ymin: 0., ymax: 3., n: 300,
                init:vec![4. / 3.]).plot(
                    |x| 1. + x * x + 0.0125 * (1. - 3. * (x - 1.)).abs().ln(),
                   "1 + x² + 0.0125 ln|1 - 3(x-1)| (specifying x:4/3)")?,
            p!(xmin: -0.5, xmax: 0.5, ymin: -1., ymax: 1.).plot(
                |x| x * (1. / x).sin(), "x sin(1/x)")?,
            p!(xmin: -0.5, xmax: 0.5, ymin: -1., ymax: 1., n:200).plot(
                |x| x * (1. / x).sin(), "x sin(1/x)")?,
            p!(xmin: -2., xmax: 2., ymin: -1., ymax: 1.).plot(
                |x| (1. / x).sin(), "sin(1/x)")?,
            p!(xmin: -2., xmax: 2., ymin: -1., ymax: 1., n: 400).plot(
                |x| (1. / x).sin(), "sin(1/x)")?,
            p!(xmin: -4., xmax: 4., ymin: -1., ymax: 1.).plot(
                |x| x.powi(4).sin(), "sin(x⁴)")?,
            p!(xmin: -4., xmax: 4., ymin: -1., ymax: 1., n: 600).plot(
                |x| x.powi(4).sin(), "sin(x⁴)")?,
            p!(xmin: -6., xmax: 6., ymin: -1., ymax: 1.).plot(
                |x| x.exp().sin(), "sin(exp x)")?,
            p!(xmin: -6., xmax: 6., ymin: -1., ymax: 1., n: 500).plot(
                |x| x.exp().sin(), "sin(exp x)")?,
            p!(xmin: -10., xmax: 10., ymin: 0., ymax: 10.).plot(
                |x| 1. / x.sin(), "1 / sin x")?,
            p!(xmin: -6., xmax: 6., ymin: -1., ymax: 1.).plot(
                |x| x.sin() / x, "(sin x)/x")?,
            p!(xmin: -2., xmax: 2., ymin: -15., ymax: 15.).plot(
                |x| (x.powi(3) - x + 1.).tan() + 1. / (x + 3. * x.exp()),
                "tan(x³ - x + 1) + 1/(x + 3 eˣ)")?,
            p!(xmin: 0., xmax: 17., ymin: 0., ymax: 2.).plot(
                |x| (1. + x.cos()) * (-0.1 * x).exp(),
                "(1 + cos x) exp(-x/10)")?,
            p!(xmin: -2., xmax: 17., ymin: 0., ymax: 2.).plot(
                |x| (1. + x.cos()) * (-0.1 * x).exp(),
                "(1 + cos x) exp(-x/10)")?,
            p!(xmin: 0., xmax: 17., ymin: 0., ymax: 2.).plot(
                |x| (1. + x.cos()) * (-0.01 * x * x).exp(),
                "(1 + cos x) exp(-x²/100)")?,
        ].join("");
        let mut fh = File::create("target/horror.gp")?;
        write!(fh, "set terminal pdfcairo\n\
                    set output \"horror.pdf\"\n\
                    set grid\n\
                    {s}\n")?;
        Ok(())
    }
}