basin 3.0.0

// Index-based loops mirror Powell's exposition (NEWUOA, Powell 2006) and the
// dense index arithmetic of the H-factorization algebra; the lint is
// blanket-allowed for this module.
#![allow(clippy::needless_range_loop)]
// The standalone `minimize` driver, `NewuoaConfig`/`NewuoaOutcome`, and the
// model-core read surface have only `#[cfg(test)]` callers (the in-module tests
// and the PRIMA parity fixtures); the public `Newuoa` solver drives the model
// through `NewuoaWork` instead. Blanket-allow so that test-only `pub(crate)`
// surface does not trip dead-code analysis in the non-test build.
#![allow(dead_code)]

//! NEWUOA (Powell 2006): model-based derivative-free optimization.
//!
//! [`Newuoa`] is Powell's NEWUOA: an unconstrained trust-region method for
//! smooth objectives whose derivatives are unavailable. It maintains a quadratic
//! surrogate `Q` interpolating the objective on a set of `npt` points and updates
//! it by the least-Frobenius-norm rule, so each iteration needs only **one** new
//! objective value. It binds [`CostFunction`] only; no gradient.
//!
//! # Architecture
//!
//! The module is the spine of Powell's model-based family (NEWUOA → BOBYQA →
//! LINCOA): an internal `QuadraticModel`, maintained through a factored
//! inverse-KKT matrix `H = W⁻¹`, plus a swappable trust-region subproblem (here
//! TRSAPP). The two operations that maintain the model are the closed-form
//! initialization (§3) and the least-Frobenius-norm rank-2 update (§4; derivation
//! in Powell 2004a, `references/frobenius-update/`). The Figure-1 driver loop adds
//! the ρ/Δ schedule (§7), the BIGLAG/BIGDEN geometry steps (§6), origin shifts
//! (§7), and the Qint robustness modification (§8). These pieces are
//! crate-internal; only [`Newuoa`] and [`PointState`](crate::PointState) are
//! public.
//!
//! # State and solver split
//!
//! [`PointState`](crate::PointState) publishes the evaluated iterate and its
//! historical incumbent. The solver owns the quadratic model, factored `H`,
//! and trust-region radii. Read the current resolution through [`Newuoa::rho`].
//! The model uses internal `Vec<F>` storage and needs no matrix capability from
//! the parameter backend. Fresh initialization rebuilds it; exact continuation
//! retains the solver and progress together.
//!
//! # Termination
//!
//! NEWUOA's natural convergence is `ρ` reaching `ρ_end` (configured on the
//! solver, where it also drives the eq-7.6 schedule); 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_radius_tolerance`](crate::Newuoa::with_absolute_radius_tolerance) to stop early at a coarser `ρ`.
//!
//! # Backends
//!
//! `Vec<F>`, nalgebra `DVector<F>`, ndarray `Array1<F>`, and faer `Col<F>`
//! all work with `F = f64` or `f32` (the parameter type needs only element access and length:
//! [`Clone`], [`VectorLen`](crate::core::math::VectorLen), and indexing). wasm-
//! clean: the model algebra is pure-Rust `Vec<F>` with no BLAS/LAPACK.
//!
//! [`CostFunction`]: crate::core::problem::CostFunction
//! [`crate::TerminationCode::SolverConverged`]: crate::TerminationCode::SolverConverged

pub(crate) mod bigden;
pub(crate) mod biglag;
pub(crate) mod driver;
pub(crate) mod init;
pub(crate) mod trsapp;

#[cfg(test)]
mod parity;

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, State};
use crate::core::termination::Termination;

use driver::{NewuoaWork, Transition};

/// NEWUOA (Powell 2006): model-based derivative-free trust-region optimization.
///
/// Configure the trust-region radii and interpolation-set size, then drive it
/// with an [`Executor`](crate::Executor) over a [`PointState`]:
///
/// ```
/// use basin::{CostFunction, Executor, Newuoa, PointState};
///
/// struct Quadratic;
/// impl CostFunction for Quadratic {
///     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).powi(2) + 2.0 * (x[1] + 2.0).powi(2))
///     }
/// }
///
/// let solver = Newuoa::new()
///     .with_initial_radius(0.5)
///     .with_final_radius(1e-8);
/// let state = PointState::new(vec![0.0, 0.0]);
/// let result = Executor::new(Quadratic, solver, state)
///     .max_cost_evals(500)
///     .run()
///     .unwrap();
/// assert!(result.best_cost() < 1e-10);
/// ```
///
/// # Configuration
///
/// - [`with_initial_radius`](Self::with_initial_radius): initial trust-region radius `ρ_beg`
///   (a reasonable initial change to the variables; default `1.0`). Also the
///   initial `Δ`.
/// - [`with_final_radius`](Self::with_final_radius): final radius `ρ_end`, the required
///   accuracy in the variables (default `1e-6`). The run converges once `ρ`
///   reaches it; `ρ_end` also drives the eq-7.6 schedule, so it must satisfy
///   `ρ_beg > ρ_end > 0`.
/// - [`with_npt`](Self::with_npt): interpolation-set size `npt`, in
///   `[n+2, ½(n+1)(n+2)]` (default `2n+1`, Powell's recommendation).
///
/// # Progress and continuation
///
/// Use [`PointState::new`] or [`crate::Executor::from_start`] with an unevaluated
/// starting point. Fresh initialization clears progress and rebuilds the
/// interpolation model, even when reusing a populated state or solver. The
/// shared state retains the current evaluated point, its historical objective
/// incumbent, and all six raw evaluation categories. It supplies no gradient
/// or radius capability. The solver owns every model buffer and the radius
/// schedule; [`rho`](Self::rho) returns `None` before initialization. Inspect
/// it through [`crate::Executor::run_with_solver`] or [`crate::Stepper::solver`].
///
/// Preserve [`crate::ExactCheckpoint`] for continuation without initialization
/// or reevaluation. With `serde`, the solver and shared state support
/// serialization, including the complete numerical model. A progress-only
/// snapshot is a fresh point warm start. Legacy algorithm-specific state
/// payloads cannot be read as `PointState` or as a new exact checkpoint.
/// Radius-based convergence remains solver-owned; objective targets and stalls
/// use the shared state's historical incumbent.
///
/// # Backends
///
/// `Vec<F>`, nalgebra `DVector<F>`, ndarray `Array1<F>`, and faer `Col<F>`,
/// for every supported backend release and `F = f64` or `f32`. The parameter
/// type needs only [`Clone`], [`VectorLen`], and `Index`/`IndexMut` element access.
/// The model uses pure-Rust `Vec<F>` storage, with no BLAS/LAPACK dependency.
/// Choose radii appropriate to the scalar precision.
///
/// # References
///
/// M. J. D. Powell, *The NEWUOA software for unconstrained optimization without
/// derivatives*, in Large-Scale Nonlinear Optimization (2006), pp. 255–297.
/// Cross-validated against [PRIMA](https://github.com/libprima/prima) v0.7.2,
/// the authoritative source for the exact formulas. PRIMA is BSD 3-Clause
/// licensed; its required notice is retained in the crate's `COPYRIGHT` file.
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Newuoa<F: Scalar = f64> {
    radius_tolerance: Option<F>,
    rho_beg: F,
    rho_end: F,
    npt: Option<usize>,
    /// Built in [`Solver::init`]; the resumable model + ρ/Δ schedule.
    work: Option<NewuoaWork<F>>,
}

impl<F: Scalar> Newuoa<F> {
    /// Current trust-region resolution, or `None` before initialization.
    /// Read it through `run_with_solver()` or `Stepper::solver()`.
    pub fn rho(&self) -> Option<F> {
        self.work.as_ref().map(|work| work.rho())
    }

    /// Copy the current interpolation geometry and native model gradient.
    ///
    /// Returns the trust radius, effective final resolution, absolute interpolation
    /// points, their sampled costs, and the model gradient at its best point.
    /// This opt-in snapshot performs no callbacks and changes no model or schedule.
    /// The gradient belongs to the interpolation model, not the true objective.
    /// Returns `None` before successful initialization.
    #[allow(clippy::type_complexity)]
    pub fn geometry_diagnostics(
        &self,
    ) -> Option<(F, F, Vec<Vec<F>>, Vec<F>, Vec<F>)> {
        self.work.as_ref().map(|work| work.geometry_diagnostics())
    }

    /// 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_radius_tolerance(
        mut self,
        value: impl Into<Option<F>>,
    ) -> Self {
        self.radius_tolerance =
            crate::core::convergence::optional_tolerance(value);
        self
    }

    /// A NEWUOA solver with the default schedule (`ρ_beg = 1`, `ρ_end = 1e-6`,
    /// `npt = 2n+1`). Tune with the `with_*` builders.
    pub fn new() -> Self {
        Self {
            rho_beg: F::from_f64(1.0).expect("1.0 representable"),
            radius_tolerance: None,
            rho_end: F::from_f64(1e-6).expect("1e-6 representable"),
            npt: None,
            work: None,
        }
    }

    /// Configure the initial radius.
    /// Retains the algorithm's existing formula, validation, and default.
    pub fn with_initial_radius(mut self, rho_beg: F) -> Self {
        self.rho_beg = rho_beg;
        self
    }

    /// Configure the final radius.
    /// Retains the algorithm's existing formula, validation, and default.
    pub fn with_final_radius(mut self, rho_end: F) -> Self {
        self.rho_end = rho_end;
        self
    }

    /// Set the interpolation-set size `npt`, in `[n+2, ½(n+1)(n+2)]`. Defaults to
    /// `2n+1` (Powell's recommendation) when unset.
    pub fn with_npt(mut self, npt: usize) -> Self {
        self.npt = Some(npt);
        self
    }
}

impl<F: Scalar> Default for Newuoa<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
}

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

impl<P, V, F> Solver<P, PointState<V, F>> for Newuoa<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;
        let n = state.param().vec_len();
        assert!(n >= 1, "Newuoa requires a non-empty start point");
        let npt = self.npt.unwrap_or(2 * n + 1);

        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))
            };
            NewuoaWork::try_init(
                x0,
                self.rho_beg,
                self.rho_end,
                npt,
                &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("Newuoa::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 evaluations into the running best (the reported
        // iterate). NEWUOA's reported point is the least-F seen, 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);
            }
        }

        // ρ reaching ρ_end (with the work there complete) is NEWUOA's natural
        // convergence; signal it as a clean stop.
        let reason = match out.transition {
            Transition::Converged => Some(Termination::custom(
                "newuoa.final_radius",
                "The algorithm completed its final trust-radius stage; this does not certify feasibility or global optimality.",
                vec![crate::Measurement {
                    name: "final_radius".into(),
                    value: self.rho_end,
                }],
            )),
            Transition::Continue | Transition::RhoReduced => None,
        };
        Ok(crate::SolverStep::from((state, reason)))
    }

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