basin 3.0.0

use crate::core::constraint::evaluation::BoxProblem;
use crate::core::convergence::optional_tolerance;
use crate::core::math::{Scalar, VectorIndex, VectorLen};
use crate::core::parallel::{MaybeSend, MaybeSync};
use crate::core::problem::{CostFunction, Problem};
use crate::core::solver::Solver;
use crate::core::state::{PointState, State};
use crate::core::termination::Termination;

mod rectangle;
use rectangle::{Rectangle, Work};

/// Original DIRECT (DIviding RECTangles) for deterministic global optimization
/// over finite box bounds, following Jones, Perttunen, and Stuckman (1993).
///
/// # Algorithm
///
/// The free coordinates are normalized to the unit hypercube. Initialization
/// evaluates its midpoint once. Each iteration selects all potentially optimal
/// rectangles using their center costs `f` and Euclidean half-diagonals `d`:
/// a rectangle `j` qualifies when some `K > 0` satisfies both
/// `f_j - K d_j <= f_i - K d_i` for every `i` and
/// `f_j - K d_j <= f_min - epsilon * abs(f_min)`.
/// The implementation groups equal sizes and constructs their lower convex
/// hull, retaining collinear eligible points and all exact cost ties.
///
/// The selected set is frozen for the iteration and processed from largest
/// to smallest, with creation order breaking ties. Each rectangle is sampled
/// at `c ± delta e_i` along every longest axis, where `delta` is one-third of
/// that side's length. Probes are ordered by axis, negative then positive.
/// Trisection proceeds along axes in ascending order of their best probe
/// cost, with axis index breaking ties. Each split divides the remaining
/// central rectangle, so the better samples occupy larger children. Integer
/// trisection depths identify equal sizes without floating-point comparisons.
///
/// This is original DIRECT, not the locally biased DIRECT-L variant. SciPy's
/// corresponding setting is `locally_biased=False`. Basin uses exact cost
/// ties; SciPy 1.16.2 groups costs within an absolute `1e-13`, so complete
/// trajectories and evaluation counts need not coincide.
///
/// # Initialization and continuation
///
/// Supply a [`PointState`] with the same dimension as the bounds. Its parameter
/// provides a backend template; its values do not select the starting point.
/// Every fresh run starts at the box midpoint and clears the partition. This
/// solver does not implement point warm starts. An entirely fixed box requires
/// just one evaluation. Exact continuation uses a solver-and-state
/// [`ExactCheckpoint`](crate::ExactCheckpoint), preserving the partition and
/// evaluation counts. With `serde`, the solver's evolving machinery is also
/// serializable. Bounds and the objective must remain unchanged during a run
/// and across exact continuation.
///
/// # Termination
///
/// By default the solver stops when the best rectangle's half-diagonal in
/// normalized coordinates is at most `1e-6`. Relative volume stopping is
/// disabled by default because volumes shrink quickly with dimension. The
/// enabled geometry checks compose with OR and return
/// [`SolverConverged`](crate::TerminationCode::SolverConverged). They measure
/// search resolution and do not certify global optimality.
///
/// Use [`Executor`](crate::Executor) for budgets, objective targets,
/// cancellation, and stagnation. Budgets and cancellation are checked between
/// complete subdivision sweeps, so the last sweep can exceed an evaluation
/// budget. Storage grows as `O(n * evaluations)`; cost ties can make sweeps
/// large. Independent probes use [`Problem::cost_batch`], preserving numerical
/// results with the optional `parallel` feature. Actual callback scheduling
/// may differ under parallel execution.
///
/// # Numerical safeguards
///
/// Bounds must be finite and ordered. Fixed coordinates are removed from the
/// normalized geometry and restored in every callback. Interpolation avoids
/// overflowing the width of a box spanning extreme finite values. If a split
/// cannot produce distinct probes in both normalized and original coordinates,
/// the solver returns [`NumericalNoProgress`](crate::TerminationCode::NumericalNoProgress)
/// without evaluating those probes.
///
/// NaN and positive infinity reject a sample. As an extension for rejected
/// samples, each sweep additionally explores the oldest largest rejected
/// rectangle; its cost is excluded from the finite selection diagram. This
/// permits recovery from a rejected midpoint. Geometry stopping requires a
/// usable incumbent. An all-rejected run stops through executor controls with
/// no incumbent, or reports [`SolverFailed`](crate::TerminationCode::SolverFailed)
/// for an entirely fixed box. Negative infinity is retained and stops with
/// `SolverConverged`. Typed objective errors propagate unchanged.
///
/// # Backends
///
/// `Vec<F>`, `nalgebra::DVector<F>`, `ndarray::Array1<F>`, and `faer::Col<F>`
/// support both `F = f64` (default) and `f32`, with their corresponding backend
/// features. Only [`VectorLen`], [`VectorIndex<F>`], and `Clone` are required.
/// With `parallel`, the problem and parameter must be `Sync`, and the output
/// and problem error must be `Send`.
///
/// # Example
///
/// ```
/// use basin::{CostFunction, Direct, Executor, PointState};
/// struct BoundedQuartic { lower: Vec<f64>, upper: Vec<f64> }
/// impl basin::Constraints for BoundedQuartic  {
///     type Set = basin::ConstraintSet<basin::BoxBounds<Vec<f64>>>;
///     fn constraints(&self) -> <Self::Set as basin::ConstraintDescription>::Borrowed<'_> { basin::ConstraintSet::new().with_bounds(basin::BoxBounds::new(Some(&self.lower), Some(&self.upper))) }
/// }
/// impl CostFunction for BoundedQuartic {
///     type Param = Vec<f64>;
///     type Output = f64;
///     type Error = std::convert::Infallible;
///     fn cost(&self, x: &Vec<f64>) -> Result<f64, Self::Error> {
///         Ok(x.iter().map(|x| 0.5 * (x.powi(4) - 16.0*x*x + 5.0*x)).sum())
///     }
/// }
///
/// let result = Executor::new(
///     BoundedQuartic { lower: vec![-5.0; 2], upper: vec![5.0; 2] },
///     Direct::new(),
///     PointState::new(vec![0.0; 2]),
/// ).target_objective(-78.33).max_cost_evals(2000).run().unwrap();
/// assert!(result.cost() < -78.33);
/// ```
///
/// # References
///
/// Donald R. Jones, Cary D. Perttunen, and Bruce E. Stuckman,
/// “Lipschitzian optimization without the Lipschitz constant,”
/// *Journal of Optimization Theory and Applications* 79, 157–181 (1993),
/// especially pp. 169–170, the subdivision procedure and Definition 4.1.
/// <https://doi.org/10.1007/BF00941892>.
///
/// Solution quality and evaluation counts are cross-checked with
/// [SciPy 1.16.2](https://github.com/scipy/scipy/tree/v1.16.2/scipy/optimize/_direct)
/// (BSD-3-Clause), whose DIRECT kernel derives from Gablonsky's DIRECT 2.0.4
/// (MIT). Basin implements the published algorithm independently.
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[derive(Clone, Debug)]
pub struct Direct<F: Scalar = f64> {
    epsilon: F,
    radius_tolerance: Option<F>,
    volume_tolerance: Option<F>,
    work: Option<Work<F>>,
}

/// Borrowed geometry of the live incumbent rectangle after a DIRECT sweep.
///
/// Available after initialization when a usable incumbent exists. Reading this
/// view performs no evaluations and changes no partition. The center and depths
/// cover only free coordinates in the unit box. A hard error can interrupt a
/// sweep and does not promise a coherent checkpoint. Fresh initialization resets
/// the partition; exact solver-aware continuation retains it.
#[derive(Clone, Copy, Debug)]
#[non_exhaustive]
pub struct RectangleRef<'a, F: Scalar = f64> {
    /// Original coordinate indices participating in normalized geometry.
    pub free_coordinates: &'a [usize],
    /// Center in normalized free coordinates.
    pub center: &'a [F],
    /// Per-coordinate trisection depths.
    pub depths: &'a [usize],
    /// Euclidean half-diagonal in normalized free coordinates.
    pub radius: F,
    /// Sum of the trisection depths, defining relative volume as `3^-level`.
    pub level: u64,
    /// Number of rectangles retained in the partition.
    pub rectangle_count: usize,
}

impl<F: Scalar> Default for Direct<F> {
    fn default() -> Self {
        Self {
            epsilon: F::from_f64(1e-4).unwrap(),
            radius_tolerance: Some(F::from_f64(1e-6).unwrap()),
            volume_tolerance: None,
            work: None,
        }
    }
}

impl<F: Scalar> Direct<F> {
    /// Borrow the current incumbent's normalized rectangle geometry.
    pub fn best_rectangle(&self) -> Option<RectangleRef<'_, F>> {
        let work = self.work.as_ref()?;
        let rectangle = &work.rectangles[work.best?];
        Some(RectangleRef {
            free_coordinates: &work.active,
            center: &rectangle.center,
            depths: &rectangle.depths,
            radius: work.radius(rectangle),
            level: rectangle.level(),
            rectangle_count: work.rectangles.len(),
        })
    }
    /// Construct original DIRECT with epsilon `1e-4`, normalized radius
    /// tolerance `1e-6`, and no volume tolerance.
    pub fn new() -> Self {
        Self::default()
    }

    /// Set the potential-improvement factor (default `1e-4`).
    ///
    /// This is a rectangle-selection control, not a stopping tolerance.
    /// Smaller values allow more local refinement. Zero removes the required
    /// improvement margin. Panics unless `epsilon` is finite and nonnegative.
    pub fn with_epsilon(mut self, epsilon: F) -> Self {
        assert!(
            epsilon.is_finite() && epsilon >= F::zero(),
            "DIRECT epsilon must be finite and nonnegative"
        );
        self.epsilon = epsilon;
        self
    }

    /// Stop at a best-rectangle half-diagonal at most `value` in normalized
    /// unit-box coordinates (default `Some(1e-6)`).
    ///
    /// This is an absolute Euclidean radius after coordinate normalization,
    /// not a radius in the original units. Fixed coordinates do not contribute.
    /// `None` disables the check and zero requests exact zero. Panics for a
    /// non-finite or negative tolerance. Enabled geometry checks compose with OR.
    pub fn with_absolute_radius_tolerance(
        mut self,
        value: impl Into<Option<F>>,
    ) -> Self {
        self.radius_tolerance = optional_tolerance(value);
        self
    }

    /// Stop when the best rectangle's volume, relative to the original box
    /// over free coordinates, is at most `value` (default `None`).
    ///
    /// `None` disables the check and zero requests exact zero. The comparison
    /// uses logarithmic volumes to avoid treating underflow as convergence.
    /// Panics for a non-finite or negative tolerance. Enabled geometry checks
    /// compose with OR.
    pub fn with_relative_volume_tolerance(
        mut self,
        value: impl Into<Option<F>>,
    ) -> Self {
        self.volume_tolerance = optional_tolerance(value);
        self
    }
}

impl<P, V, F> Solver<P, PointState<V, F>> for Direct<F>
where
    F: Scalar + MaybeSend,
    P: CostFunction<Param = V, Output = F> + BoxProblem<Param = V> + MaybeSync,
    P::Error: MaybeSend,
    V: Clone + VectorLen + VectorIndex<F> + MaybeSync,
{
    type Error = crate::ConstraintSolveError<P::Error>;

    fn init(
        &mut self,
        problem: &mut Problem<P>,
        mut state: PointState<V, F>,
    ) -> Result<PointState<V, F>, Self::Error> {
        crate::core::constraint::evaluation::validate_coordinate_box_problem(
            problem.inner(),
            state.param().vec_len(),
            true,
        )
        .map_err(crate::core::constraint::evaluation::invalid)?;
        self.work = None;
        let lower = problem.inner().lower();
        let upper = problem.inner().upper();
        let n = lower.vec_len();
        assert!(n > 0, "DIRECT requires a nonempty box");
        assert_eq!(upper.vec_len(), n, "DIRECT bound dimensions differ");
        assert_eq!(
            state.param().vec_len(),
            n,
            "DIRECT state and bound dimensions differ"
        );
        let mut bounds = Vec::with_capacity(n);
        let mut active = Vec::new();
        for i in 0..n {
            let lo = lower.get_scalar(i);
            let hi = upper.get_scalar(i);
            assert!(
                lo.is_finite() && hi.is_finite() && lo <= hi,
                "DIRECT requires finite ordered bounds at coordinate {i}"
            );
            bounds.push((lo, hi));
            if lo < hi {
                active.push(i);
            }
        }
        let root = Rectangle::new(
            vec![F::from_f64(0.5).unwrap(); active.len()],
            F::infinity(),
        );
        let mut work = Work::new(bounds, active, root);
        let point = work.point(state.param(), &work.rectangles[0].center);
        let cost = problem.cost(&point)?;
        work.rectangles[0].cost = cost;
        work.consider(0);
        state.reset();
        state.replace(point, cost);
        self.work = Some(work);
        Ok(state)
    }

    fn next_iter(
        &mut self,
        problem: &mut Problem<P>,
        mut state: PointState<V, F>,
    ) -> Result<crate::SolverStep<PointState<V, F>>, Self::Error> {
        let work = self.work.as_mut().expect("DIRECT must be initialized");
        let selected = work.select(self.epsilon);
        if selected.is_empty() {
            return Ok(crate::SolverStep::from((
                state,
                Some(Termination::numerical_stall(
                    "direct: no numerical progress",
                )),
            )));
        }
        for id in selected {
            let Some(probes) = work.probes(state.param(), id) else {
                return Ok(crate::SolverStep::from((
                    state,
                    Some(Termination::numerical_stall(
                        "direct: no numerical progress",
                    )),
                )));
            };
            let costs = problem.cost_batch(&probes.params)?;
            work.divide(id, probes.axes, probes.centers, costs);
            if let Some(best) = work.best {
                let rectangle = &work.rectangles[best];
                state.replace(
                    work.point(state.param(), &rectangle.center),
                    rectangle.cost,
                );
                if rectangle.cost == F::neg_infinity() {
                    return Ok(crate::SolverStep::from((
                        state,
                        Some(Termination::custom(
                            "direct.negative_infinity",
                            "The selected rectangle has objective negative infinity.",
                            vec![crate::Measurement {
                                name: "cost".into(),
                                value: rectangle.cost,
                            }],
                        )),
                    )));
                }
            }
        }
        Ok(crate::SolverStep::from((state, None)))
    }

    fn terminate(&self, _state: &PointState<V, F>) -> Option<Termination<F>> {
        let work = self.work.as_ref()?;
        let Some(best) = work.best else {
            return work.active.is_empty().then(|| Termination::numerical_failure("DIRECT exhausted its rectangles without a valid incumbent."));
        };
        let rectangle = &work.rectangles[best];
        if rectangle.cost == F::neg_infinity() || work.active.is_empty() {
            return Some(Termination::custom(
                "direct.exhausted_or_unbounded",
                "The selected objective is negative infinity or every splittable rectangle is exhausted.",
                vec![
                    crate::Measurement {
                        name: "cost".into(),
                        value: rectangle.cost,
                    },
                    crate::Measurement {
                        name: "active_rectangles".into(),
                        value: F::from_usize(work.active.len()).unwrap(),
                    },
                ],
            ));
        }
        let radius = work.radius(rectangle);
        let radius_small = self.radius_tolerance.is_some_and(|t| radius <= t);
        let volume_small = self.volume_tolerance.is_some_and(|t| {
            F::from_u64(rectangle.level()).unwrap()
                * F::from_f64(3.0).unwrap().ln()
                >= -t.ln()
        });
        let radius_report = radius_small.then(|| {
            Termination::custom(
                "direct.radius",
                "The scaled best-rectangle radius is <= the configured tolerance, including exact zero.",
                vec![
                    crate::Measurement {
                        name: "radius".into(),
                        value: radius,
                    },
                    crate::Measurement {
                        name: "tolerance".into(),
                        value: self.radius_tolerance.unwrap(),
                    },
                ],
            )
        });
        let volume_report = volume_small.then(|| {
            Termination::custom(
                "direct.volume",
                "level * ln(3) >= -ln(volume_tolerance); a zero threshold requires a truly zero volume, not underflow.",
                vec![
                    crate::Measurement {
                        name: "level".into(),
                        value: F::from_u64(rectangle.level()).unwrap(),
                    },
                    crate::Measurement {
                        name: "tolerance".into(),
                        value: self.volume_tolerance.unwrap(),
                    },
                ],
            )
        });
        radius_report
            .into_iter()
            .chain(volume_report)
            .reduce(Termination::merge)
    }
}