xad-rs 0.2.0

//! `LabeledHessian` + `compute_full_hessian` — string-keyed dense full
//! Hessian helper, gated on the `labeled-hessian` feature.
//!
//! This module is the only code in `xad-rs` v0.2.0 with compile-time
//! dependencies on both the labeled variable layer ([`VarRegistry`])
//! and [`Dual2Vec`](crate::dual2vec::Dual2Vec). It offers a labeled
//! convenience wrapper around a single `Dual2Vec` forward pass: build a
//! [`VarRegistry`] from the input name slice, seed one
//! [`Dual2Vec::variable`](crate::dual2vec::Dual2Vec::variable) per
//! position, invoke the user closure, and package the output
//! `(value, gradient, hessian)` alongside the registry.
//!
//! # Complexity and mode-selection guidance
//!
//! Each arithmetic operation on `Dual2Vec` is `O(n²)` in the number of
//! input variables (the cost of updating the `n × n` Hessian
//! accumulator). For problems with `n` larger than roughly `50-100`,
//! prefer seeded [`Dual2<T>`](crate::dual2::Dual2) in a loop over input
//! directions — see the [`dual2vec`](crate::dual2vec) module docs for
//! the full crossover discussion. `compute_full_hessian` inherits that
//! cost model verbatim; this module cross-references the `Dual2Vec`
//! docs rather than duplicating the guidance.
//!
//! # Feature gate
//!
//! This entire module is gated on `labeled-hessian`, which transitively
//! enables both `labeled` (for `VarRegistry`) and `dual2-vec` (for
//! `Dual2Vec`). The slim default build (`--no-default-features`) does
//! not link `indexmap` or `ndarray`.
//!
//! # No lifetime parameters
//!
//! Neither [`LabeledHessian`] nor [`compute_full_hessian`] carries any
//! lifetime parameter. The returned struct is `'static`, so it can be
//! freely stored, returned from functions, and passed between threads.
//! This matches the PyO3-future constraint in `PROJECT.md` line 107:
//! no lifetime parameters leaking into user-facing labeled types.
//!
//! # Example
//!
//! ```
//! # #[cfg(feature = "labeled-hessian")]
//! # {
//! use xad_rs::dual2vec::Dual2Vec;
//! use xad_rs::labeled::hessian::compute_full_hessian;
//!
//! // f(x, y) = x²·y + y³ at (1, 2) — Hessian [[4, 2], [2, 12]]
//! let inputs = vec![("x".to_string(), 1.0), ("y".to_string(), 2.0)];
//! let lh = compute_full_hessian(&inputs, |v: &[Dual2Vec]| {
//!     let x = &v[0];
//!     let y = &v[1];
//!     &(&(x * x) * y) + &(&(y * y) * y)
//! });
//!
//! assert_eq!(lh.vars.index_of("x"), Some(0));
//! assert_eq!(lh.vars.index_of("y"), Some(1));
//! assert!((lh.value - 10.0).abs() < 1e-13);
//! assert!((lh.hessian[[0, 0]] - 4.0).abs() < 1e-13);
//! assert!((lh.hessian[[1, 1]] - 12.0).abs() < 1e-13);
//! assert!((lh.hessian[[0, 1]] - 2.0).abs() < 1e-13);
//! # }
//! ```

use ndarray::Array2;

use crate::dual2vec::Dual2Vec;
use crate::labeled::VarRegistry;

/// Labeled dense full Hessian produced by [`compute_full_hessian`].
///
/// Field semantics:
///
/// - `vars` — frozen [`VarRegistry`] whose insertion order matches the
///   `inputs` slice passed to [`compute_full_hessian`]. Use
///   [`VarRegistry::index_of`] to map a name to the positional slot
///   used by `gradient` and `hessian`.
/// - `value` — `f(inputs)` evaluated at the seed point.
/// - `gradient` — `∂f/∂inputs[i]` as a `Vec<f64>` of length `n`.
/// - `hessian` — `∂²f/∂inputs[i]∂inputs[j]` as an `n × n`
///   [`ndarray::Array2<f64>`]. Structurally symmetric straight out of
///   the `Dual2Vec` arithmetic pipeline (see `Dual2Vec` docs D2VS-01);
///   no `symmetrize()` call is needed or performed.
///
/// No lifetime parameters — `LabeledHessian` is `'static` and can be
/// freely stored, returned from functions, and passed between threads.
pub struct LabeledHessian {
    pub vars: VarRegistry,
    pub value: f64,
    pub gradient: Vec<f64>,
    pub hessian: Array2<f64>,
}

/// Compute a labeled dense full Hessian via a single `Dual2Vec` forward
/// pass.
///
/// # Arguments
///
/// * `inputs` — `(name, value)` pairs in the order they should appear
///   as rows and columns of the Hessian. Names become the
///   `LabeledHessian.vars` registry; duplicate names are idempotent
///   (first insertion wins — see [`VarRegistry::from_names`]).
/// * `f` — closure operating on positional `&[Dual2Vec]` and returning
///   the output `Dual2Vec`. Each input is seeded via
///   [`Dual2Vec::variable(value, i, n)`](Dual2Vec::variable).
///
/// # Complexity
///
/// Each operation inside `f` costs `O(n²)` where `n = inputs.len()`.
/// For `n` larger than roughly `50-100`, prefer seeding a
/// [`Dual2<T>`](crate::dual2::Dual2) per input direction in a loop —
/// see the [`dual2vec`](crate::dual2vec) module docs for the full
/// discussion.
pub fn compute_full_hessian<F>(inputs: &[(String, f64)], f: F) -> LabeledHessian
where
    F: Fn(&[Dual2Vec]) -> Dual2Vec,
{
    let n = inputs.len();

    // Build the frozen VarRegistry from input names in insertion order.
    let names: Vec<String> = inputs.iter().map(|(name, _)| name.clone()).collect();
    let vars = VarRegistry::from_names(names);

    // Seed the Dual2Vec inputs at positions 0..n.
    let d2v_inputs: Vec<Dual2Vec> = inputs
        .iter()
        .enumerate()
        .map(|(i, (_, value))| Dual2Vec::variable(*value, i, n))
        .collect();

    // Run the user closure.
    let output = f(&d2v_inputs);

    // Extract value, gradient, hessian into owned types.
    let value = output.value();
    let gradient: Vec<f64> = output.gradient().to_vec();
    let hessian: Array2<f64> = output.hessian().clone();

    LabeledHessian {
        vars,
        value,
        gradient,
        hessian,
    }
}