basin 3.0.0

// Index-based loops mirror the OrthoMADS exposition (Abramson et al. 2009) and
// the flat poll/direction arithmetic; the lint is blanket-allowed here.
#![allow(clippy::needless_range_loop)]

//! MADS (Audet & Dennis 2006): mesh adaptive direct search, the deterministic
//! OrthoMADS instance (Abramson, Audet, Dennis & Le Digabel 2009).
//!
//! [`Mads`] is a directional direct-search method for **nonsmooth or
//! non-continuous** objectives, where Nelder-Mead has no convergence theory and
//! is known to stagnate. Each iteration polls a positive spanning set of `2n`
//! directions on a shrinking mesh around the incumbent; an improving mesh point
//! coarsens the mesh, a failed poll refines it. The set of poll directions
//! becomes asymptotically dense, which gives MADS its Clarke-stationarity
//! convergence guarantee on nonsmooth functions. It binds [`CostFunction`] only;
//! no derivatives.
//!
//! # Architecture
//!
//! This implements the **OrthoMADS** instance: the poll directions are generated
//! deterministically (no RNG) from the quasi-random Halton sequence (`halton`)
//! and a scaled Householder reflection (`geometry`), producing an orthogonal
//! integer basis whose `±` columns form the maximal positive basis
//! `D_k = [H −H]`. Two parameters drive the mesh: a fine mesh size `Δᵐ` and a
//! coarser poll size `Δᵖ`, linked through the integer mesh index `ℓ` by
//! `Δᵖ = 2^{-ℓ}`, `Δᵐ = min{1, 4^{-ℓ}}` (OrthoMADS eq. (1)); `Δᵐ` shrinks faster
//! than `Δᵖ`, so the number of available poll directions grows without bound.
//! Because the directions are integer and trial points are `x_k + Δᵐ·d`, every
//! poll point lies exactly on the mesh. The poll loop and bookkeeping live in
//! the `driver` submodule.
//!
//! # State and solver split
//!
//! Shared progress retains the evaluated incumbent and its publication metadata.
//! Unbounded and box-bounded modes use [`PointState`](crate::PointState);
//! progressive-barrier mode uses [`SelectedState`](crate::SelectedState) with
//! the matching constraint violation. The solver owns its mesh, Halton schedule,
//! poll size, and feasible/infeasible search history. Read diagnostics through
//! [`Mads::poll_size`] and [`Mads::mesh_index`].
//!
//! # Termination
//!
//! MADS's natural convergence is the poll size `Δᵖ` reaching the configured floor
//! `poll_size_min`; the solver signals it via
//! [`crate::TerminationCode::SolverConverged`]. Add [`max_cost_evals`](crate::Executor::max_cost_evals)
//! to cap the evaluation budget, or [`with_absolute_poll_size_tolerance`](crate::Mads::with_absolute_poll_size_tolerance) to
//! stop early at a coarser poll size.
//!
//! # Backends
//!
//! `Vec<F>`, nalgebra `DVector<F>`, ndarray `Array1<F>`, and faer `Col<F>`
//! work with `F = f64` or `f32` (the parameter type needs only [`Clone`], [`VectorLen`], and
//! `Index`/`IndexMut`). The poll geometry is pure-Rust integer/`f64` scratch with
//! no RNG, no `linalg`-tier ops, and no time: fully deterministic and wasm-clean.
//!
//! [`CostFunction`]: crate::core::problem::CostFunction
//! [`VectorLen`]: crate::core::math::VectorLen
//! [`crate::TerminationCode::SolverConverged`]: crate::TerminationCode::SolverConverged

pub(crate) mod driver;
pub(crate) mod geometry;
pub(crate) mod halton;
pub(crate) mod progressive_barrier;

use std::marker::PhantomData;

use crate::core::constraint::evaluation::{BoxProblem, invalid};
use crate::core::inner::InitialState;
use crate::core::math::{Scalar, VectorLen};
use crate::core::problem::{CostFunction, Problem};
use crate::core::solver::Solver;
use crate::core::state::{PointState, SelectedState, State};
use crate::core::termination::Termination;
use crate::solver::nelder_mead::Unbounded;
use crate::{ConstraintSolveError, NonlinearInequalityValues};
use crate::{ViolationMeasure, ViolationReport};

use driver::{MadsWork, Transition};
use progressive_barrier::PbWork;

/// Type-state marker: box-bound-constrained MADS (`Mads<Bounded>`), reached via
/// [`Mads::bounded`]. Bounds are enforced by the **extreme barrier**: infeasible
/// poll points are assigned `f = +∞` and rejected without evaluating the
/// objective (Audet & Dennis 2006, §2). Pair with a problem implementing
/// [`BoxBounds`](crate::BoxBounds).
pub struct Bounded;

/// Type-state marker: nonlinearly-constrained MADS (`Mads<Constrained>`),
/// reached via [`Mads::constrained`]. General nonlinear inequality constraints
/// `c(x) ≤ 0` are handled by the **progressive barrier** (Audet & Dennis 2009):
/// the solver tracks the aggregate violation `h(x) = Σⱼ max(cⱼ, 0)²` and walks a
/// barrier threshold down to zero, so an **infeasible start is allowed**. Pair
/// with a problem implementing [`NonlinearInequalityValues`].
pub struct Constrained;

/// MADS (Audet & Dennis 2006): mesh adaptive direct search, deterministic
/// OrthoMADS instance: derivative-free local optimization for nonsmooth or
/// non-continuous objectives.
///
/// Configure the poll-size schedule, then drive it with an
/// [`Executor`](crate::Executor) over a [`PointState`]:
///
/// ```
/// use basin::{CostFunction, Executor, Mads, PointState};
///
/// // A nonsmooth objective (an L1 "valley") with minimizer (1, 2).
/// struct AbsValley;
/// impl CostFunction for AbsValley {
///     type Param = Vec<f64>;
///     type Output = f64;
///     type Error = std::convert::Infallible;
///     fn cost(&self, x: &Vec<f64>) -> Result<f64, std::convert::Infallible> {
///         Ok((x[0] - 1.0).abs() + 2.0 * (x[1] - 2.0).abs())
///     }
/// }
///
/// let solver = Mads::new()
///     .with_initial_poll_size(1.0)
///     .with_minimum_poll_size(1e-9);
/// let state = PointState::new(vec![0.0, 0.0]);
/// let result = Executor::new(AbsValley, solver, state)
///     .max_cost_evals(5_000)
///     .run()
///     .unwrap();
/// assert!(result.best_cost() < 1e-6);
/// ```
///
/// # Configuration
///
/// - [`with_initial_poll_size`](Self::with_initial_poll_size): initial poll
///   size `Δ₀` (a reasonable initial change to the variables; default `1.0`). It
///   uniformly scales the mesh, so trial points stay on the scaled integer mesh.
/// - [`with_minimum_poll_size`](Self::with_minimum_poll_size): convergence floor on the
///   poll size (default `1e-6`). The run converges once `Δᵖ` reaches it; must be
///   `> 0` and `< Δ₀`.
///
/// # ConstraintSet
///
/// Two opt-in constrained modes specialize the default unconstrained solver:
///
/// - [`bounded`](Self::bounded) → `Mads<Bounded>` handles box bounds by the
///   **extreme barrier** (infeasible poll points get `f = +∞`); needs a problem
///   declaring a box-only set and a feasible (or clamped) start.
/// - [`constrained`](Self::constrained) → `Mads<Constrained>` handles general
///   nonlinear inequality constraints `c(x) ≤ 0` by the **progressive barrier**
///   (Audet & Dennis 2009): it tracks the aggregate violation
///   `h(x) = Σⱼ max(cⱼ, 0)²` and drives a barrier threshold down to zero,
///   exploring around both a feasible and an infeasible incumbent. Needs a
///   problem implementing [`NonlinearInequalityValues`] and a
///   [`SelectedState`]; an **infeasible start is allowed**.
///
/// # Progress and continuation
///
/// Use [`PointState::new`] for unbounded and box-bounded runs, and
/// [`SelectedState::new`] for progressive-barrier runs. [`crate::Executor::from_start`]
/// chooses the matching state. Shared progress preserves every raw evaluation
/// category. Constraint callbacks count as residual evaluations, including
/// failures; box-rejected probes make no objective call. The constrained
/// state's `current()` and `best()` return the selected point, objective, and
/// violation report with [`ViolationMeasure::SquaredL2`], unscaled and with
/// identity transformation. Its `value` is the squared-sum violation. NaN
/// residuals do not satisfy the positive-part comparison and contribute zero;
/// reporting preserves this native behavior. [`crate::ConstraintViolationState`]
/// borrows current and incumbent reports without evaluations.
/// Selection can increase the objective as feasibility
/// improves; unchanged records retain their publication metadata.
///
/// The solver owns every mesh and Halton counter and the progressive barrier's
/// history. [`poll_size`](Self::poll_size) and [`mesh_index`](Self::mesh_index)
/// return `None` before initialization; read them via
/// [`crate::Executor::run_with_solver`] or [`crate::Stepper::solver`]. Fresh
/// runs reset progress and rebuild those histories, reevaluating the seed.
/// Exact continuation retains a [`crate::ExactCheckpoint`] and skips init.
/// With `serde`, every mode serializes its model and matching shared state.
/// Legacy `MadsState` and `ConstrainedMadsState` payloads are incompatible.
/// Native and configured poll-size stopping still read the solver's schedule;
/// shared progress states do not expose [`crate::MeshState`].
///
/// # Backends
///
/// `Vec<F>`, nalgebra `DVector<F>`, ndarray `Array1<F>`, and faer `Col<F>`,
/// for every supported release and `F = f64` or `f32`, in every mode. Vectors
/// need only [`Clone`], [`VectorLen`], and `Index`/`IndexMut`. Poll geometry
/// uses internal integer and floating-point buffers, with no matrix capability
/// or BLAS/LAPACK dependency. The method is deterministic and works on WASM.
/// Choose poll-size limits appropriate to the scalar precision.
///
/// # References
///
/// C. Audet & J. E. Dennis, Jr., *Mesh adaptive direct search algorithms for
/// constrained optimization*, SIAM J. Optim. 17 (2006), pp. 188–217 (with the
/// 2008 erratum). M. A. Abramson, C. Audet, J. E. Dennis, Jr. & S. Le Digabel,
/// *OrthoMADS: A deterministic MADS instance with orthogonal directions*, SIAM
/// J. Optim. 20 (2009), pp. 948–966. C. Audet & J. E. Dennis, Jr., *A
/// progressive barrier for derivative-free nonlinear programming*, SIAM J.
/// Optim. 20 (2009), pp. 445–472 (the constrained mode).
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[cfg_attr(
    feature = "serde",
    serde(bound(
        serialize = "F: serde::Serialize",
        deserialize = "F: serde::Deserialize<'de>"
    ))
)]
pub struct Mads<Mode = Unbounded, F: Scalar = f64> {
    poll_tolerance: Option<F>,
    poll_size_init: F,
    poll_size_min: F,
    /// Built in [`Solver::init`]; the resumable poll loop + mesh schedule
    /// (unconstrained or box-bounded modes).
    work: Option<MadsWork<F>>,
    /// Built in [`Solver::init`] for the constrained mode; the resumable
    /// progressive-barrier poll loop. `None` in the other modes.
    pb_work: Option<PbWork<F>>,
    constraint_count: Option<usize>,
    _mode: PhantomData<fn() -> Mode>,
}

impl<Mode, F: Scalar> Mads<Mode, F> {
    /// Current poll size, or `None` before initialization.
    pub fn poll_size(&self) -> Option<F> {
        self.work
            .as_ref()
            .map(MadsWork::poll_size)
            .or_else(|| self.pb_work.as_ref().map(PbWork::poll_size))
    }

    /// Current mesh index, or `None` before initialization.
    pub fn mesh_index(&self) -> Option<i32> {
        self.work
            .as_ref()
            .map(MadsWork::mesh_index)
            .or_else(|| self.pb_work.as_ref().map(PbWork::mesh_index))
    }

    /// Stop when the observed radius or step size is <= the tolerance.
    ///
    /// Disabled by default. `None` disables the test and zero requests an
    /// exact-zero threshold. The tolerance must be finite and nonnegative.
    /// Checked at initialized iteration boundaries. This observation does not
    /// change the algorithm's radius or step-size update schedule.
    pub fn with_absolute_poll_size_tolerance(
        mut self,
        value: impl Into<Option<F>>,
    ) -> Self {
        self.poll_tolerance =
            crate::core::convergence::optional_tolerance(value);
        self
    }
}

impl<F: Scalar> Mads<Unbounded, F> {
    /// A MADS solver with the default schedule (`Δ₀ = 1`, `poll_size_min =
    /// 1e-6`). Tune with the `with_*` builders.
    pub fn new() -> Self {
        Self {
            poll_size_init: F::from_f64(1.0).expect("1.0 representable"),
            poll_tolerance: None,
            poll_size_min: F::from_f64(1e-6).expect("1e-6 representable"),
            work: None,
            pb_work: None,
            constraint_count: None,
            _mode: PhantomData,
        }
    }
}

impl<Mode, F: Scalar> Mads<Mode, F> {
    /// Set the initial poll size `Δ₀`, a reasonable initial change to the
    /// variables. Uniformly scales the mesh.
    pub fn with_initial_poll_size(mut self, poll_size_init: F) -> Self {
        self.poll_size_init = poll_size_init;
        self
    }

    /// Configure the minimum poll size.
    /// Retains the algorithm's existing formula, validation, and default.
    pub fn with_minimum_poll_size(mut self, poll_size_min: F) -> Self {
        self.poll_size_min = poll_size_min;
        self
    }
}

impl<F: Scalar> Mads<Unbounded, F> {
    /// Switch to the box-bound-constrained variant ([`Bounded`]). The resulting
    /// `Mads<Bounded>` requires a problem declaring a box-only set and
    /// enforces the bounds by the extreme barrier (infeasible poll points get
    /// `f = +∞`).
    /// An infeasible start is clamped into the box.
    pub fn bounded(self) -> Mads<Bounded, F> {
        Mads {
            poll_size_init: self.poll_size_init,
            poll_tolerance: self.poll_tolerance,
            poll_size_min: self.poll_size_min,
            work: None,
            pb_work: None,
            constraint_count: None,
            _mode: PhantomData,
        }
    }

    /// Switch to the nonlinearly-constrained variant ([`Constrained`]). The
    /// resulting `Mads<Constrained>` requires a problem implementing
    /// [`NonlinearInequalityValues`] (`c(x) ≤ 0`) and handles the
    /// constraints by the **progressive barrier**, so an **infeasible start is
    /// allowed**. Drive it with a [`SelectedState`].
    pub fn constrained(self) -> Mads<Constrained, F> {
        Mads {
            poll_size_init: self.poll_size_init,
            poll_tolerance: self.poll_tolerance,
            poll_size_min: self.poll_size_min,
            work: None,
            pb_work: None,
            constraint_count: None,
            _mode: PhantomData,
        }
    }
}

impl<F: Scalar> Default for Mads<Unbounded, F> {
    fn default() -> Self {
        Self::new()
    }
}

/// Build a `V` from a flat `&[F]`, reusing `template` for type and length.
fn fill_from<V, F>(template: &V, slice: &[F]) -> V
where
    V: Clone + std::ops::IndexMut<usize, Output = F>,
    F: Copy,
{
    let mut v = template.clone();
    for (i, &x) in slice.iter().enumerate() {
        v[i] = x;
    }
    v
}

/// Copy a backend vector `V` into a flat `Vec<F>` of length `n`.
fn to_vec<V, F>(v: &V, n: usize) -> Vec<F>
where
    V: std::ops::Index<usize, Output = F>,
    F: Copy,
{
    (0..n).map(|i| v[i]).collect()
}

/// Whether any component of `x` lies outside `[lower, upper]`.
fn out_of_box<F: Scalar>(x: &[F], lower: &[F], upper: &[F]) -> bool {
    x.iter()
        .zip(lower)
        .zip(upper)
        .any(|((&xi, &lo), &up)| xi < lo || xi > up)
}

impl<V, F> InitialState<V> for Mads<Unbounded, F>
where
    F: Scalar,
    V: Clone,
{
    type State = PointState<V, F>;
    fn seed(&self, x: &V) -> Self::State {
        PointState::new(x.clone())
    }
}

impl<V, F> InitialState<V> for Mads<Bounded, F>
where
    F: Scalar,
    V: Clone,
{
    type State = PointState<V, F>;
    fn seed(&self, x: &V) -> Self::State {
        PointState::new(x.clone())
    }
}

impl<V, F> InitialState<V> for Mads<Constrained, F>
where
    F: Scalar,
    V: Clone,
{
    type State = SelectedState<V, F>;
    fn seed(&self, x: &V) -> Self::State {
        SelectedState::new(x.clone())
    }
}

impl<P, V, F> Solver<P, PointState<V, F>> for Mads<Unbounded, F>
where
    F: Scalar,
    P: CostFunction<Param = V, Output = F>,
    V: Clone
        + VectorLen
        + std::ops::Index<usize, Output = F>
        + std::ops::IndexMut<usize, Output = F>,
{
    type Error = P::Error;

    fn init(
        &mut self,
        problem: &mut Problem<P>,
        mut state: PointState<V, F>,
    ) -> Result<PointState<V, F>, Self::Error> {
        state.reset();
        self.work = None;
        self.pb_work = None;
        let n = state.param().vec_len();
        assert!(n >= 1, "Mads requires a non-empty start point");

        let x0: Vec<F> = (0..n).map(|i| state.param()[i]).collect();
        let template = state.param().clone();

        let (work, best_x, best_f) = {
            let mut eval = |slice: &[F]| -> Result<F, P::Error> {
                problem.cost(&fill_from(&template, slice))
            };
            MadsWork::try_init(
                x0,
                self.poll_size_init,
                self.poll_size_min,
                &mut eval,
            )?
        };

        state.replace(fill_from(&template, &best_x), best_f);
        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 template = state.param().clone();
        let work = self
            .work
            .as_mut()
            .expect("Mads::init must run before next_iter");

        let out = {
            let mut eval = |slice: &[F]| -> Result<F, P::Error> {
                problem.cost(&fill_from(&template, slice))
            };
            work.step(&mut eval)?
        };

        // Fold the step's improving evaluation into the running best (the
        // reported iterate). MADS accepts only improving mesh points, so
        // param/cost are monotone and coincide with best_param/best_cost.
        let mut best_f = state.cost();
        for (xabs, f_new) in &out.evaluated {
            if *f_new < best_f {
                best_f = *f_new;
                state.replace(fill_from(&template, xabs), best_f);
            }
        }

        // The poll size reaching its floor is MADS's natural convergence.
        let reason = match out.transition {
            Transition::Converged => Some(Termination::custom(
                "mads.final_poll",
                "The poll schedule reached its configured minimum size.",
                vec![crate::Measurement {
                    name: "minimum_poll_size".into(),
                    value: self.poll_size_min,
                }],
            )),
            Transition::Continue => None,
        };
        Ok(crate::SolverStep::from((state, reason)))
    }

    fn terminate(&self, _state: &PointState<V, F>) -> Option<Termination<F>> {
        let tolerance = self.poll_tolerance?;
        let metric = self.poll_size()?;
        (metric.is_finite() && metric <= tolerance).then(|| {
            Termination::upper_bound(
                crate::ConvergenceTest::Mesh,
                metric,
                tolerance,
                tolerance,
                None,
            )
        })
    }
}

/// Aggregate constraint violation `h(x) = Σⱼ max(cⱼ(x), 0)²` from the constraint
/// vector `c(x)` (length `m`). `h = 0` exactly on the feasible set `c(x) ≤ 0`.
fn violation<V, F>(c: &V, m: usize) -> F
where
    V: std::ops::Index<usize, Output = F>,
    F: Scalar,
{
    let mut h = F::zero();
    for j in 0..m {
        let cj = c[j];
        if cj > F::zero() {
            h = h + cj * cj;
        }
    }
    h
}

impl<P, V, F> Solver<P, SelectedState<V, F>> for Mads<Constrained, F>
where
    F: Scalar,
    P: CostFunction<Param = V, Output = F>
        + crate::Constraints<
            Set = crate::ConstraintSet<
                crate::NoBounds,
                crate::NoLinearConstraints,
                crate::NonlinearBlocks<
                    crate::NoNonlinearBlock,
                    crate::NonlinearBlock,
                >,
            >,
        > + NonlinearInequalityValues<Param = V, Error = <P as CostFunction>::Error>,
    <P as NonlinearInequalityValues>::Values:
        VectorLen + std::ops::Index<usize, Output = F>,
    V: Clone
        + VectorLen
        + std::ops::Index<usize, Output = F>
        + std::ops::IndexMut<usize, Output = F>,
{
    type Error = ConstraintSolveError<<P as CostFunction>::Error>;

    fn init(
        &mut self,
        problem: &mut Problem<P>,
        mut state: SelectedState<V, F>,
    ) -> Result<SelectedState<V, F>, Self::Error> {
        state.reset();
        self.work = None;
        self.pb_work = None;
        let n = state.param().vec_len();
        assert!(n >= 1, "Mads requires a non-empty start point");

        let m = crate::Constraints::constraints(problem.inner())
            .nonlinear()
            .inequalities()
            .count();
        self.constraint_count = Some(m);
        let x0: Vec<F> = (0..n).map(|i| state.param()[i]).collect();
        let template = state.param().clone();

        let (work, best_x, best_f, best_h) = {
            let mut eval = |slice: &[F]| -> Result<
                (F, F),
                ConstraintSolveError<<P as CostFunction>::Error>,
            > {
                let xv = fill_from(&template, slice);
                let f = problem.cost(&xv)?;
                let c = checked_inequalities(problem, &xv, m)?;
                Ok((f, violation(&c, m)))
            };
            PbWork::try_init(
                x0,
                self.poll_size_init,
                self.poll_size_min,
                &mut eval,
            )?
        };

        state.replace(
            fill_from(&template, &best_x),
            best_f,
            ViolationReport::new(best_h, ViolationMeasure::SquaredL2),
        );
        state.select_current();
        self.pb_work = Some(work);
        Ok(state)
    }

    fn next_iter(
        &mut self,
        problem: &mut Problem<P>,
        mut state: SelectedState<V, F>,
    ) -> Result<crate::SolverStep<SelectedState<V, F>>, Self::Error> {
        let m = self
            .constraint_count
            .expect("MADS initialization declares the constraint count");
        let template = state.param().clone();
        let work = self
            .pb_work
            .as_mut()
            .expect("Mads::init must run before next_iter");

        let out = {
            let mut eval = |slice: &[F]| -> Result<
                (F, F),
                ConstraintSolveError<<P as CostFunction>::Error>,
            > {
                let xv = fill_from(&template, slice);
                let f = problem.cost(&xv)?;
                let c = checked_inequalities(problem, &xv, m)?;
                Ok((f, violation(&c, m)))
            };
            work.step(&mut eval)?
        };

        // The progressive-barrier driver owns incumbent selection; adopt its
        // reported incumbent (feasible if found, else best infeasible).
        let (x_star, f_star, h_star) = out.incumbent;
        let changed = state.current().is_none_or(|(x, f, h)| {
            f != f_star
                || h.value != h_star
                || (0..x_star.len()).any(|i| x[i] != x_star[i])
        });
        state.replace(
            fill_from(&template, &x_star),
            f_star,
            ViolationReport::new(h_star, ViolationMeasure::SquaredL2),
        );
        if changed {
            state.select_current();
        }

        let reason = match out.transition {
            Transition::Converged => Some(Termination::custom(
                "mads.final_poll",
                "The poll schedule reached its configured minimum size.",
                vec![crate::Measurement {
                    name: "minimum_poll_size".into(),
                    value: self.poll_size_min,
                }],
            )),
            Transition::Continue => None,
        };
        Ok(crate::SolverStep::from((state, reason)))
    }

    fn terminate(
        &self,
        _state: &SelectedState<V, F>,
    ) -> Option<Termination<F>> {
        let tolerance = self.poll_tolerance?;
        let metric = self.poll_size()?;
        (metric.is_finite() && metric <= tolerance).then(|| {
            Termination::upper_bound(
                crate::ConvergenceTest::Mesh,
                metric,
                tolerance,
                tolerance,
                None,
            )
        })
    }
}

impl<P, V, F> Solver<P, PointState<V, F>> for Mads<Bounded, F>
where
    F: Scalar,
    P: CostFunction<Param = V, Output = F> + BoxProblem<Param = V>,
    V: Clone
        + VectorLen
        + std::ops::Index<usize, Output = F>
        + std::ops::IndexMut<usize, Output = F>,
    V: crate::VectorLen + std::ops::Index<usize, Output = F>,
{
    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_box_problem(
            problem.inner(),
            Some(state.param().vec_len()),
            false,
        )
        .map_err(crate::core::constraint::evaluation::invalid)?;
        state.reset();
        self.work = None;
        self.pb_work = None;
        let n = state.param().vec_len();
        assert!(n >= 1, "Mads requires a non-empty start point");

        let lower = to_vec(problem.inner().lower(), n);
        let upper = to_vec(problem.inner().upper(), n);
        // Clamp an infeasible start into the box so the extreme barrier has a
        // feasible incumbent to descend from.
        let x0: Vec<F> = (0..n)
            .map(|i| {
                let xi = state.param()[i];
                if xi < lower[i] {
                    lower[i]
                } else if xi > upper[i] {
                    upper[i]
                } else {
                    xi
                }
            })
            .collect();
        let template = state.param().clone();

        let (work, best_x, best_f) = {
            let mut eval = |slice: &[F]| -> Result<
                F,
                crate::ConstraintSolveError<P::Error>,
            > {
                if out_of_box(slice, &lower, &upper) {
                    Ok(F::infinity()) // extreme barrier: reject without evaluating
                } else {
                    Ok(problem.cost(&fill_from(&template, slice))?)
                }
            };
            MadsWork::try_init(
                x0,
                self.poll_size_init,
                self.poll_size_min,
                &mut eval,
            )?
        };

        state.replace(fill_from(&template, &best_x), best_f);
        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 n = state.param().vec_len();
        let lower = to_vec(problem.inner().lower(), n);
        let upper = to_vec(problem.inner().upper(), n);
        let template = state.param().clone();
        let work = self
            .work
            .as_mut()
            .expect("Mads::init must run before next_iter");

        let out = {
            let mut eval = |slice: &[F]| -> Result<
                F,
                crate::ConstraintSolveError<P::Error>,
            > {
                if out_of_box(slice, &lower, &upper) {
                    Ok(F::infinity())
                } else {
                    Ok(problem.cost(&fill_from(&template, slice))?)
                }
            };
            work.step(&mut eval)?
        };

        let mut best_f = state.cost();
        for (xabs, f_new) in &out.evaluated {
            if *f_new < best_f {
                best_f = *f_new;
                state.replace(fill_from(&template, xabs), best_f);
            }
        }

        let reason = match out.transition {
            Transition::Converged => Some(Termination::custom(
                "mads.final_poll",
                "The poll schedule reached its configured minimum size.",
                vec![crate::Measurement {
                    name: "minimum_poll_size".into(),
                    value: self.poll_size_min,
                }],
            )),
            Transition::Continue => None,
        };
        Ok(crate::SolverStep::from((state, reason)))
    }

    fn terminate(&self, _state: &PointState<V, F>) -> Option<Termination<F>> {
        let tolerance = self.poll_tolerance?;
        let metric = self.poll_size()?;
        (metric.is_finite() && metric <= tolerance).then(|| {
            Termination::upper_bound(
                crate::ConvergenceTest::Mesh,
                metric,
                tolerance,
                tolerance,
                None,
            )
        })
    }
}

fn checked_inequalities<P, V, F>(
    p: &mut Problem<P>,
    x: &V,
    m: usize,
) -> Result<Vec<F>, ConstraintSolveError<<P as CostFunction>::Error>>
where
    P: CostFunction<Param = V, Output = F>
        + crate::Constraints<
            Set = crate::ConstraintSet<
                crate::NoBounds,
                crate::NoLinearConstraints,
                crate::NonlinearBlocks<
                    crate::NoNonlinearBlock,
                    crate::NonlinearBlock,
                >,
            >,
        > + NonlinearInequalityValues<Param = V, Error = <P as CostFunction>::Error>,
    <P as NonlinearInequalityValues>::Values:
        VectorLen + std::ops::Index<usize, Output = F>,
    F: Scalar,
{
    if crate::Constraints::constraints(p.inner())
        .nonlinear()
        .inequalities()
        .count()
        != m
    {
        return Err(invalid(crate::ConstraintError::LayoutChanged));
    }
    if m == 0 {
        return Ok(vec![]);
    }
    p.counts_mut().residual_evals += 1;
    let values =
        <P as NonlinearInequalityValues>::nonlinear_inequalities(p.inner(), x)?;
    crate::core::constraint::model::dimension(
        crate::ConstraintBlock::NonlinearInequalities,
        crate::ConstraintAxis::Length,
        m,
        values.vec_len(),
    )
    .map_err(invalid)?;
    Ok((0..m).map(|i| values[i]).collect())
}