math_audio_optimisation/
cobyla.rs

1//! Pure-Rust COBYLA local optimizer.
2//!
3//! Wraps the [`cobyla`](https://crates.io/crates/cobyla) crate (a faithful
4//! pure-Rust port of NLopt's `cobyla.c`, itself a port of M.J.D. Powell's
5//! 1994 FORTRAN implementation) behind an API consistent with this crate's
6//! [`differential_evolution`](crate::differential_evolution) and
7//! [`levenberg_marquardt`](crate::levenberg_marquardt) entry points.
8//!
9//! # Constraint convention
10//!
11//! Inequality constraints are passed in the **autoeq/NLopt convention**:
12//! `g_i(x) <= 0` is feasible. Internally this module flips the sign before
13//! handing them to the `cobyla` crate, which uses the opposite convention
14//! (`fc(x) >= 0` is feasible).
15//!
16//! # Example
17//!
18//! ```rust
19//! use math_audio_optimisation::cobyla::{cobyla, CobylaConfig, CobylaConstraint, CobylaRhoBegin};
20//! use ndarray::Array1;
21//! use std::sync::Arc;
22//!
23//! // Minimize 10*(x0+1)^2 + x1^2 subject to x0 >= 0
24//! // (i.e. constraint  -x0 <= 0  in our convention)
25//! let f = |x: &Array1<f64>| 10.0 * (x[0] + 1.0).powi(2) + x[1].powi(2);
26//! let constraints = vec![CobylaConstraint {
27//!     fun: Arc::new(|x: &Array1<f64>| -x[0]),
28//! }];
29//! let cfg = CobylaConfig {
30//!     x0: Array1::from(vec![1.0, 1.0]),
31//!     bounds: vec![(-10.0, 10.0), (-10.0, 10.0)],
32//!     maxeval: 200,
33//!     rho_begin: CobylaRhoBegin::All(0.5),
34//!     ..Default::default()
35//! };
36//! let report = cobyla(&f, &constraints, cfg).expect("cobyla should run");
37//! assert!(report.x[0] >= -1e-3, "x0 must be roughly >= 0");
38//! assert!(report.fun < 10.1, "constrained minimum is at (0, 0) with value 10");
39//! ```
40
41use ndarray::Array1;
42use std::sync::Arc;
43
44use crate::error::{DEError, Result};
45
46/// Erased inequality-constraint closure. Returns `<= 0` when feasible
47/// (NLopt convention; the wrapper flips the sign for the `cobyla` crate).
48pub type CobylaConstraintFn = Arc<dyn Fn(&Array1<f64>) -> f64 + Send + Sync>;
49
50/// A single inequality constraint `fun(x) <= 0`.
51#[derive(Clone)]
52pub struct CobylaConstraint {
53    /// Constraint function. Feasible when `<= 0`.
54    pub fun: CobylaConstraintFn,
55}
56
57/// Termination tolerances passed through to NLopt's `nlopt_stopping`.
58#[derive(Debug, Clone, Copy)]
59pub struct CobylaStopTols {
60    /// Stop when objective drops at or below `stopval`. Disabled if
61    /// non-positive (mapped to `-INFINITY` internally). The upstream
62    /// `cobyla` crate's `minimize()` hardcoded this to `-INFINITY` —
63    /// re-exposing it as a tuning knob.
64    pub stopval: f64,
65    /// Stop when |Δf| < `ftol_abs`. Disabled if non-positive.
66    pub ftol_abs: f64,
67    /// Stop when |Δf|/|f| < `ftol_rel`. Disabled if non-positive.
68    pub ftol_rel: f64,
69    /// Stop when ||Δx||_inf < `xtol_abs`. Disabled if non-positive.
70    pub xtol_abs: f64,
71    /// Stop when ||Δx||/||x|| < `xtol_rel`. Disabled if non-positive.
72    pub xtol_rel: f64,
73}
74
75impl Default for CobylaStopTols {
76    fn default() -> Self {
77        // Match NLopt-via-autoeq settings (the previous `optim/nlopt.rs`
78        // before the C-FFI dep was dropped):
79        //   set_stopval(1e-4); set_ftol_rel(1e-6); set_xtol_rel(1e-4);
80        Self {
81            stopval: 1e-4,
82            ftol_abs: 0.0,
83            ftol_rel: 1e-6,
84            xtol_abs: 0.0,
85            xtol_rel: 1e-4,
86        }
87    }
88}
89
90/// Initial trust-region radius spec.
91///
92/// Powell's COBYLA accepts a single `rhobeg` applied to all variables, or
93/// per-dimension initial steps. For autoeq's PEQ problem the parameters
94/// differ in scale by orders of magnitude (log10-frequency span ~0.5,
95/// gain span ~24 dB), so per-dimension steps are recommended.
96#[derive(Debug, Clone)]
97pub enum CobylaRhoBegin {
98    /// Same initial step for every dimension.
99    All(f64),
100    /// Per-dimension initial steps. Length must equal `x0.len()`.
101    PerDim(Vec<f64>),
102}
103
104impl Default for CobylaRhoBegin {
105    fn default() -> Self {
106        Self::All(0.5)
107    }
108}
109
110/// Configuration for [`cobyla`].
111#[derive(Clone)]
112pub struct CobylaConfig {
113    /// Initial guess. Must lie within `bounds`.
114    pub x0: Array1<f64>,
115    /// `(lower, upper)` per dimension.
116    pub bounds: Vec<(f64, f64)>,
117    /// Initial trust-region radius. See [`CobylaRhoBegin`].
118    pub rho_begin: CobylaRhoBegin,
119    /// Maximum number of objective evaluations.
120    pub maxeval: usize,
121    /// Termination tolerances.
122    pub stop_tol: CobylaStopTols,
123}
124
125impl Default for CobylaConfig {
126    fn default() -> Self {
127        Self {
128            x0: Array1::zeros(0),
129            bounds: Vec::new(),
130            rho_begin: CobylaRhoBegin::default(),
131            maxeval: 1000,
132            stop_tol: CobylaStopTols::default(),
133        }
134    }
135}
136
137/// Result of a [`cobyla`] run.
138#[derive(Clone)]
139pub struct CobylaReport {
140    /// Best parameter vector found.
141    pub x: Array1<f64>,
142    /// Objective value at `x`.
143    pub fun: f64,
144    /// Whether the underlying solver returned a success status.
145    pub success: bool,
146    /// Human-readable status (Debug-formatted underlying status).
147    pub message: String,
148    /// Number of function evaluations consumed (best-effort — the
149    /// underlying crate does not expose an exact count, so this is the
150    /// configured `maxeval` budget when unknown).
151    pub nfev: usize,
152}
153
154impl std::fmt::Debug for CobylaReport {
155    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
156        f.debug_struct("CobylaReport")
157            .field("x_len", &self.x.len())
158            .field("fun", &self.fun)
159            .field("success", &self.success)
160            .field("message", &self.message)
161            .field("nfev", &self.nfev)
162            .finish()
163    }
164}
165
166/// Minimize `f` subject to bounds and inequality constraints `g_i(x) <= 0`.
167///
168/// On success returns the best point found. On a non-success termination
169/// the report is still populated with the last point evaluated so callers
170/// can decide whether the result is usable.
171pub fn cobyla<F>(
172    f: &F,
173    constraints: &[CobylaConstraint],
174    config: CobylaConfig,
175) -> Result<CobylaReport>
176where
177    F: Fn(&Array1<f64>) -> f64 + Sync,
178{
179    let n = config.x0.len();
180    if n == 0 {
181        return Err(DEError::BoundsMismatch {
182            lower_len: 0,
183            upper_len: 0,
184        });
185    }
186    if config.bounds.len() != n {
187        return Err(DEError::BoundsMismatch {
188            lower_len: config.bounds.len(),
189            upper_len: n,
190        });
191    }
192    for (i, (lo, hi)) in config.bounds.iter().enumerate() {
193        if lo > hi {
194            return Err(DEError::InvalidBounds {
195                index: i,
196                lower: *lo,
197                upper: *hi,
198            });
199        }
200    }
201
202    // Translate config into the `cobyla_native` API (hand-translated
203    // Rust port of NLopt 2.7.1's cobyla.c, replacing the c2rust output).
204    use crate::cobyla_native as native;
205
206    // Build per-dim initial steps (rhobeg).
207    let dx: Vec<f64> = match &config.rho_begin {
208        CobylaRhoBegin::All(r) => vec![*r; n],
209        CobylaRhoBegin::PerDim(v) => {
210            if v.len() != n {
211                return Err(DEError::BoundsMismatch {
212                    lower_len: v.len(),
213                    upper_len: n,
214                });
215            }
216            v.clone()
217        }
218    };
219
220    // Stop criteria — hand-port honours these directly (no wrapper bugs to
221    // work around).
222    let stopval = if config.stop_tol.stopval > 0.0 {
223        config.stop_tol.stopval
224    } else {
225        f64::NEG_INFINITY
226    };
227    let stop = native::StopCriteria {
228        stopval,
229        ftol_abs: config.stop_tol.ftol_abs,
230        ftol_rel: config.stop_tol.ftol_rel,
231        xtol_abs: config.stop_tol.xtol_abs,
232        xtol_rel: config.stop_tol.xtol_rel,
233        maxeval: config.maxeval,
234    };
235
236    // Wrap user-facing closures (Array1) as &[f64] for the native API.
237    // Constraint convention is identical: feasible when `g(x) <= 0`.
238    let f_native = |x: &[f64]| -> f64 {
239        let arr = Array1::from(x.to_vec());
240        f(&arr)
241    };
242    type SliceObj = Box<dyn Fn(&[f64]) -> f64>;
243    let cons_native: Vec<SliceObj> = constraints
244        .iter()
245        .map(|c| {
246            let g = c.fun.clone();
247            let b: SliceObj = Box::new(move |x: &[f64]| {
248                let arr = Array1::from(x.to_vec());
249                g(&arr)
250            });
251            b
252        })
253        .collect();
254
255    let mut x_buf: Vec<f64> = config.x0.to_vec();
256    let report = native::cobyla_native(
257        n,
258        f_native,
259        &cons_native,
260        &config.bounds,
261        &mut x_buf,
262        &dx,
263        &stop,
264    )?;
265
266    Ok(CobylaReport {
267        x: Array1::from(report.x),
268        fun: report.fun,
269        success: report.success,
270        message: report.message.to_string(),
271        nfev: report.nfev,
272    })
273}
274
275#[cfg(test)]
276mod tests {
277    use super::*;
278
279    #[test]
280    fn paraboloid_unconstrained() {
281        // f(x) = (x0+1)^2 + x1^2, minimum at (-1, 0) with value 0.
282        let f = |x: &Array1<f64>| (x[0] + 1.0).powi(2) + x[1].powi(2);
283        let cfg = CobylaConfig {
284            x0: Array1::from(vec![1.0, 1.0]),
285            bounds: vec![(-10.0, 10.0), (-10.0, 10.0)],
286            rho_begin: CobylaRhoBegin::All(0.5),
287            maxeval: 500,
288            // Disable stopval so the test can drive convergence below 1e-4.
289            stop_tol: CobylaStopTols {
290                stopval: 0.0,
291                ..CobylaStopTols::default()
292            },
293        };
294        let report = cobyla(&f, &[], cfg).expect("cobyla failed");
295        assert!(report.fun < 1e-6, "fun = {} should be ~0", report.fun);
296        assert!((report.x[0] - (-1.0)).abs() < 1e-3);
297        assert!(report.x[1].abs() < 1e-3);
298    }
299
300    #[test]
301    fn paraboloid_with_inequality() {
302        // Minimize (x0+1)^2 + x1^2 subject to x0 >= 0 (i.e. -x0 <= 0).
303        // Constrained optimum at (0, 0) with value 1.
304        let f = |x: &Array1<f64>| (x[0] + 1.0).powi(2) + x[1].powi(2);
305        let constraints = vec![CobylaConstraint {
306            fun: Arc::new(|x: &Array1<f64>| -x[0]),
307        }];
308        let cfg = CobylaConfig {
309            x0: Array1::from(vec![1.0, 1.0]),
310            bounds: vec![(-10.0, 10.0), (-10.0, 10.0)],
311            rho_begin: CobylaRhoBegin::All(0.5),
312            maxeval: 500,
313            ..Default::default()
314        };
315        let report = cobyla(&f, &constraints, cfg).expect("cobyla failed");
316        // Active constraint at x0 = 0 ⇒ f = 1.
317        assert!(
318            (report.fun - 1.0).abs() < 1e-3,
319            "fun = {} should be ~1",
320            report.fun
321        );
322        assert!(
323            report.x[0] >= -1e-3,
324            "x0 = {} should respect x0 >= 0",
325            report.x[0]
326        );
327    }
328
329    #[test]
330    fn rejects_bad_bounds() {
331        let f = |x: &Array1<f64>| x[0] * x[0];
332        let cfg = CobylaConfig {
333            x0: Array1::from(vec![0.0]),
334            bounds: vec![(1.0, -1.0)], // lo > hi
335            rho_begin: CobylaRhoBegin::All(0.1),
336            maxeval: 100,
337            ..Default::default()
338        };
339        let err = cobyla(&f, &[], cfg).unwrap_err();
340        matches!(err, DEError::InvalidBounds { .. });
341    }
342}
math_audio_optimisation/cobyla.rs

math_audio_optimisation/
cobyla.rs
