numra_optim/
problem.rs

1//! Declarative optimization problem builder.
2//!
3//! Author: Moussa Leblouba
4//! Date: 5 March 2026
5//! Modified: 2 May 2026
6
7use numra_core::Scalar;
8
9use crate::augmented_lagrangian::AugLagOptions;
10use crate::error::OptimError;
11use crate::global::DEOptions;
12use crate::types::{OptimOptions, OptimResult};
13
14/// Boxed scalar function `&[S] -> S`.
15pub type ScalarFn<S> = Box<dyn Fn(&[S]) -> S>;
16/// Boxed vector function `(&[S], &mut [S])`.
17pub type VectorFn<S> = Box<dyn Fn(&[S], &mut [S])>;
18
19/// Variable type for mixed-integer problems.
20#[derive(Clone, Copy, Debug, PartialEq, Eq)]
21pub enum VarType {
22    Continuous,
23    Integer,
24    Binary,
25}
26
27/// Hint about problem structure for closure-based objectives.
28///
29/// When the objective is provided as a closure (`ObjectiveKind::Minimize`),
30/// the auto-selector cannot determine if the underlying function is linear
31/// or quadratic. Setting a hint informs solver selection without requiring
32/// coefficient extraction.
33#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
34pub enum ProblemHint {
35    /// No structural information (default).
36    #[default]
37    None,
38    /// The closure implements a linear function.
39    Linear,
40    /// The closure implements a quadratic function.
41    Quadratic,
42}
43
44/// Kind of constraint: equality h(x)=0 or inequality g(x)<=0.
45#[derive(Clone, Copy, Debug, PartialEq, Eq)]
46pub enum ConstraintKind {
47    Equality,
48    Inequality,
49}
50
51/// A single constraint with optional gradient.
52pub struct Constraint<S: Scalar> {
53    pub(crate) func: ScalarFn<S>,
54    pub(crate) grad: Option<VectorFn<S>>,
55    pub(crate) kind: ConstraintKind,
56}
57
58/// A linear constraint: a^T x {<=, =} b.
59#[derive(Clone, Debug)]
60pub struct LinearConstraint<S: Scalar> {
61    pub(crate) a: Vec<S>,
62    pub(crate) b: S,
63    pub(crate) kind: ConstraintKind,
64}
65
66/// Whether the objective is a scalar function or a least-squares residual.
67pub enum ObjectiveKind<S: Scalar> {
68    /// Linear objective: min c^T x
69    Linear { c: Vec<S> },
70    /// Quadratic objective: min 1/2 x^T H x + c^T x
71    /// H stored row-major as `Vec<S>` (n*n), must be positive semi-definite.
72    Quadratic {
73        h_row_major: Vec<S>,
74        c: Vec<S>,
75        n: usize,
76    },
77    Minimize {
78        func: ScalarFn<S>,
79        grad: Option<VectorFn<S>>,
80    },
81    LeastSquares {
82        residual: VectorFn<S>,
83        jacobian: Option<VectorFn<S>>,
84        n_residuals: usize,
85    },
86    /// Multi-objective: minimize [f1(x), f2(x), ...] simultaneously.
87    MultiObjective { funcs: Vec<ScalarFn<S>> },
88}
89
90/// Declarative optimization problem.
91///
92/// Build a problem with the fluent API, then call `.solve()` to automatically
93/// dispatch to the best solver, or `.solve_with(choice)` for explicit control.
94pub struct OptimProblem<S: Scalar> {
95    pub(crate) n: usize,
96    pub(crate) x0: Option<Vec<S>>,
97    pub(crate) bounds: Vec<Option<(S, S)>>,
98    pub(crate) var_types: Vec<VarType>,
99    pub(crate) objective: Option<ObjectiveKind<S>>,
100    pub(crate) constraints: Vec<Constraint<S>>,
101    pub(crate) linear_constraints: Vec<LinearConstraint<S>>,
102    pub(crate) options: OptimOptions<S>,
103    pub(crate) aug_lag_options: Option<AugLagOptions<S>>,
104    pub(crate) global: bool,
105    pub(crate) de_options: Option<DEOptions<S>>,
106    pub(crate) negate_result: bool,
107    pub(crate) hint: ProblemHint,
108}
109
110impl<S: Scalar> OptimProblem<S> {
111    /// Create a new problem with `n` decision variables.
112    pub fn new(n: usize) -> Self {
113        Self {
114            n,
115            x0: None,
116            bounds: vec![None; n],
117            var_types: vec![VarType::Continuous; n],
118            objective: None,
119            constraints: Vec::new(),
120            linear_constraints: Vec::new(),
121            options: OptimOptions::default(),
122            aug_lag_options: None,
123            global: false,
124            de_options: None,
125            negate_result: false,
126            hint: ProblemHint::None,
127        }
128    }
129
130    /// Set the initial point.
131    pub fn x0(mut self, x0: &[S]) -> Self {
132        self.x0 = Some(x0.to_vec());
133        self
134    }
135
136    /// Set bounds for variable `i`.
137    pub fn bounds(mut self, i: usize, lo_hi: (S, S)) -> Self {
138        self.bounds[i] = Some(lo_hi);
139        self
140    }
141
142    /// Set bounds for all variables at once.
143    pub fn all_bounds(mut self, bounds: &[(S, S)]) -> Self {
144        for (i, &b) in bounds.iter().enumerate() {
145            self.bounds[i] = Some(b);
146        }
147        self
148    }
149
150    /// Mark variable `i` as integer-valued.
151    pub fn integer_var(mut self, i: usize) -> Self {
152        self.var_types[i] = VarType::Integer;
153        self
154    }
155
156    /// Mark variable `i` as binary (0 or 1).
157    pub fn binary_var(mut self, i: usize) -> Self {
158        self.var_types[i] = VarType::Binary;
159        self.bounds[i] = Some((S::ZERO, S::ONE));
160        self
161    }
162
163    /// Set a scalar objective function to minimize.
164    pub fn objective<F: Fn(&[S]) -> S + 'static>(mut self, f: F) -> Self {
165        match &mut self.objective {
166            Some(ObjectiveKind::Minimize { func, .. }) => {
167                *func = Box::new(f);
168            }
169            _ => {
170                self.objective = Some(ObjectiveKind::Minimize {
171                    func: Box::new(f),
172                    grad: None,
173                });
174            }
175        }
176        self
177    }
178
179    /// Set the gradient for the scalar objective.
180    pub fn gradient<G: Fn(&[S], &mut [S]) + 'static>(mut self, g: G) -> Self {
181        match &mut self.objective {
182            Some(ObjectiveKind::Minimize { grad, .. }) => {
183                *grad = Some(Box::new(g));
184            }
185            _ => {
186                self.objective = Some(ObjectiveKind::Minimize {
187                    func: Box::new(|_| S::ZERO),
188                    grad: Some(Box::new(g)),
189                });
190            }
191        }
192        self
193    }
194
195    /// Set a scalar objective function to maximize.
196    ///
197    /// Internally converts to minimization by negating the objective and gradient.
198    pub fn maximize<F: Fn(&[S]) -> S + 'static>(mut self, f: F) -> Self {
199        self.objective = Some(ObjectiveKind::Minimize {
200            func: Box::new(move |x: &[S]| -f(x)),
201            grad: None,
202        });
203        self.negate_result = true;
204        self
205    }
206
207    /// Set the gradient for a maximization objective.
208    ///
209    /// Note: provide the gradient of the ORIGINAL function (not negated).
210    /// The builder handles negation internally.
211    pub fn maximize_gradient<G: Fn(&[S], &mut [S]) + 'static>(mut self, g: G) -> Self {
212        if let Some(ObjectiveKind::Minimize { grad, .. }) = &mut self.objective {
213            if self.negate_result {
214                *grad = Some(Box::new(move |x: &[S], gout: &mut [S]| {
215                    g(x, gout);
216                    for gi in gout.iter_mut() {
217                        *gi = -*gi;
218                    }
219                }));
220            }
221        }
222        self
223    }
224
225    /// Set a least-squares objective: minimize ||r(x)||^2.
226    pub fn least_squares<R: Fn(&[S], &mut [S]) + 'static>(
227        mut self,
228        n_residuals: usize,
229        residual: R,
230    ) -> Self {
231        self.objective = Some(ObjectiveKind::LeastSquares {
232            residual: Box::new(residual),
233            jacobian: None,
234            n_residuals,
235        });
236        self
237    }
238
239    /// Set the Jacobian for a least-squares objective.
240    pub fn jacobian<J: Fn(&[S], &mut [S]) + 'static>(mut self, j: J) -> Self {
241        if let Some(ObjectiveKind::LeastSquares { jacobian, .. }) = &mut self.objective {
242            *jacobian = Some(Box::new(j));
243        }
244        self
245    }
246
247    /// Set a linear objective: minimize c^T x.
248    pub fn linear_objective(mut self, c: &[S]) -> Self {
249        assert_eq!(c.len(), self.n, "linear objective dimension mismatch");
250        self.objective = Some(ObjectiveKind::Linear { c: c.to_vec() });
251        self
252    }
253
254    /// Set a quadratic objective: minimize 1/2 x^T H x + c^T x.
255    /// `h_row_major` is the n x n Hessian in row-major order.
256    pub fn quadratic_objective_dense(mut self, h_row_major: &[S], c: &[S]) -> Self {
257        let n = self.n;
258        assert_eq!(h_row_major.len(), n * n, "Hessian dimension mismatch");
259        assert_eq!(c.len(), n, "linear term dimension mismatch");
260        self.objective = Some(ObjectiveKind::Quadratic {
261            h_row_major: h_row_major.to_vec(),
262            c: c.to_vec(),
263            n,
264        });
265        self
266    }
267
268    /// Add a linear inequality constraint: a^T x <= b.
269    pub fn linear_constraint_ineq(mut self, a: &[S], b: S) -> Self {
270        assert_eq!(a.len(), self.n, "constraint dimension mismatch");
271        self.linear_constraints.push(LinearConstraint {
272            a: a.to_vec(),
273            b,
274            kind: ConstraintKind::Inequality,
275        });
276        self
277    }
278
279    /// Add a linear equality constraint: a^T x = b.
280    pub fn linear_constraint_eq(mut self, a: &[S], b: S) -> Self {
281        assert_eq!(a.len(), self.n, "constraint dimension mismatch");
282        self.linear_constraints.push(LinearConstraint {
283            a: a.to_vec(),
284            b,
285            kind: ConstraintKind::Equality,
286        });
287        self
288    }
289
290    /// Add an inequality constraint g(x) <= 0.
291    pub fn constraint_ineq<F: Fn(&[S]) -> S + 'static>(mut self, f: F) -> Self {
292        self.constraints.push(Constraint {
293            func: Box::new(f),
294            grad: None,
295            kind: ConstraintKind::Inequality,
296        });
297        self
298    }
299
300    /// Add an equality constraint h(x) = 0.
301    pub fn constraint_eq<F: Fn(&[S]) -> S + 'static>(mut self, f: F) -> Self {
302        self.constraints.push(Constraint {
303            func: Box::new(f),
304            grad: None,
305            kind: ConstraintKind::Equality,
306        });
307        self
308    }
309
310    /// Add an inequality constraint with gradient.
311    pub fn constraint_ineq_with_grad<F, G>(mut self, f: F, g: G) -> Self
312    where
313        F: Fn(&[S]) -> S + 'static,
314        G: Fn(&[S], &mut [S]) + 'static,
315    {
316        self.constraints.push(Constraint {
317            func: Box::new(f),
318            grad: Some(Box::new(g)),
319            kind: ConstraintKind::Inequality,
320        });
321        self
322    }
323
324    /// Add an equality constraint with gradient.
325    pub fn constraint_eq_with_grad<F, G>(mut self, f: F, g: G) -> Self
326    where
327        F: Fn(&[S]) -> S + 'static,
328        G: Fn(&[S], &mut [S]) + 'static,
329    {
330        self.constraints.push(Constraint {
331            func: Box::new(f),
332            grad: Some(Box::new(g)),
333            kind: ConstraintKind::Equality,
334        });
335        self
336    }
337
338    /// Set maximum iterations.
339    pub fn max_iter(mut self, n: usize) -> Self {
340        self.options.max_iter = n;
341        self
342    }
343
344    /// Set gradient tolerance.
345    pub fn gtol(mut self, tol: S) -> Self {
346        self.options.gtol = tol;
347        self
348    }
349
350    /// Override all options.
351    pub fn options(mut self, opts: OptimOptions<S>) -> Self {
352        self.options = opts;
353        self
354    }
355
356    /// Set augmented Lagrangian options for constrained problems.
357    pub fn aug_lag_options(mut self, opts: AugLagOptions<S>) -> Self {
358        self.aug_lag_options = Some(opts);
359        self
360    }
361
362    /// Set the constraint tolerance for constrained problems.
363    pub fn ctol(mut self, tol: S) -> Self {
364        let opts = self
365            .aug_lag_options
366            .get_or_insert_with(AugLagOptions::default);
367        opts.ctol = tol;
368        self
369    }
370
371    /// Enable global optimization (Differential Evolution).
372    pub fn global(mut self, enabled: bool) -> Self {
373        self.global = enabled;
374        self
375    }
376
377    /// Set Differential Evolution options.
378    pub fn de_options(mut self, opts: DEOptions<S>) -> Self {
379        self.de_options = Some(opts);
380        self
381    }
382
383    /// Set multiple objectives for multi-objective optimization (all minimized).
384    pub fn multi_objective(mut self, funcs: Vec<ScalarFn<S>>) -> Self {
385        self.objective = Some(ObjectiveKind::MultiObjective { funcs });
386        self
387    }
388
389    /// Provide a structural hint about the closure-based objective.
390    ///
391    /// This informs auto solver selection when the objective is a closure
392    /// that implements a linear or quadratic function but cannot be decomposed
393    /// into coefficient vectors.
394    pub fn hint(mut self, hint: ProblemHint) -> Self {
395        self.hint = hint;
396        self
397    }
398
399    // ── query methods ──
400
401    /// Returns `true` if any variable has bounds.
402    pub fn has_bounds(&self) -> bool {
403        self.bounds.iter().any(|b| b.is_some())
404    }
405
406    /// Returns `true` if constraints have been added.
407    pub fn has_constraints(&self) -> bool {
408        !self.constraints.is_empty()
409    }
410
411    /// Returns `true` if the objective is a least-squares residual.
412    pub fn is_least_squares(&self) -> bool {
413        matches!(self.objective, Some(ObjectiveKind::LeastSquares { .. }))
414    }
415
416    /// Returns `true` if the objective is structurally linear (coefficient vector).
417    pub fn is_linear(&self) -> bool {
418        matches!(self.objective, Some(ObjectiveKind::Linear { .. }))
419    }
420
421    /// Returns `true` if the objective is structurally quadratic (Hessian + linear).
422    pub fn is_quadratic(&self) -> bool {
423        matches!(self.objective, Some(ObjectiveKind::Quadratic { .. }))
424    }
425
426    /// Returns `true` if a linear hint is set (closure known to be linear).
427    pub fn is_hint_linear(&self) -> bool {
428        self.hint == ProblemHint::Linear
429    }
430
431    /// Returns `true` if a quadratic hint is set (closure known to be quadratic).
432    pub fn is_hint_quadratic(&self) -> bool {
433        self.hint == ProblemHint::Quadratic
434    }
435
436    /// Returns `true` if the objective is multi-objective.
437    pub fn is_multi_objective(&self) -> bool {
438        matches!(self.objective, Some(ObjectiveKind::MultiObjective { .. }))
439    }
440
441    /// Returns `true` if linear constraints have been added.
442    pub fn has_linear_constraints(&self) -> bool {
443        !self.linear_constraints.is_empty()
444    }
445
446    /// Returns `true` if any variable is integer or binary.
447    pub fn has_integer_vars(&self) -> bool {
448        self.var_types.iter().any(|v| *v != VarType::Continuous)
449    }
450
451    /// Number of equality constraints.
452    pub fn n_eq_constraints(&self) -> usize {
453        self.constraints
454            .iter()
455            .filter(|c| c.kind == ConstraintKind::Equality)
456            .count()
457    }
458
459    /// Number of inequality constraints.
460    pub fn n_ineq_constraints(&self) -> usize {
461        self.constraints
462            .iter()
463            .filter(|c| c.kind == ConstraintKind::Inequality)
464            .count()
465    }
466}
467
468impl<S: Scalar + faer::SimpleEntity + faer::Conjugate<Canonical = S> + faer::ComplexField>
469    OptimProblem<S>
470{
471    /// Solve using automatic solver selection.
472    pub fn solve(self) -> Result<OptimResult<S>, OptimError> {
473        crate::auto::auto_minimize(self)
474    }
475
476    /// Solve using an explicitly chosen solver.
477    pub fn solve_with(
478        self,
479        choice: crate::auto::SolverChoice,
480    ) -> Result<OptimResult<S>, OptimError> {
481        crate::auto::dispatch(self, choice)
482    }
483}
484
485/// Compute a finite-difference gradient of `f` at `x` using central differences.
486pub fn finite_diff_gradient<S: Scalar>(f: &dyn Fn(&[S]) -> S, x: &[S], g: &mut [S]) {
487    let n = x.len();
488    let eps = S::from_f64(1e-8);
489    let mut xp = x.to_vec();
490    for i in 0..n {
491        let xi_orig = xp[i];
492        xp[i] = xi_orig + eps;
493        let fp = f(&xp);
494        xp[i] = xi_orig - eps;
495        let fm = f(&xp);
496        g[i] = (fp - fm) / (S::TWO * eps);
497        xp[i] = xi_orig;
498    }
499}
500
501/// Compute a finite-difference Jacobian of residual `r` at `x`.
502/// `jac` is row-major m*n.
503pub fn finite_diff_jacobian<S: Scalar>(
504    r: &dyn Fn(&[S], &mut [S]),
505    x: &[S],
506    m: usize,
507    jac: &mut [S],
508) {
509    let n = x.len();
510    let eps = S::from_f64(1e-8);
511    let mut xp = x.to_vec();
512    let mut rp = vec![S::ZERO; m];
513    let mut rm = vec![S::ZERO; m];
514    for j in 0..n {
515        let xj_orig = xp[j];
516        xp[j] = xj_orig + eps;
517        r(&xp, &mut rp);
518        xp[j] = xj_orig - eps;
519        r(&xp, &mut rm);
520        for i in 0..m {
521            jac[i * n + j] = (rp[i] - rm[i]) / (S::TWO * eps);
522        }
523        xp[j] = xj_orig;
524    }
525}
526
527#[cfg(test)]
528mod tests {
529    use super::*;
530
531    #[test]
532    fn test_builder_construction() {
533        let problem = OptimProblem::<f64>::new(2)
534            .x0(&[1.0, 2.0])
535            .objective(|x: &[f64]| x[0] * x[0] + x[1] * x[1])
536            .gradient(|x: &[f64], g: &mut [f64]| {
537                g[0] = 2.0 * x[0];
538                g[1] = 2.0 * x[1];
539            })
540            .gtol(1e-10);
541
542        assert!(!problem.has_bounds());
543        assert!(!problem.has_constraints());
544        assert!(!problem.is_least_squares());
545        assert_eq!(problem.n_eq_constraints(), 0);
546        assert_eq!(problem.n_ineq_constraints(), 0);
547    }
548
549    #[test]
550    fn test_builder_with_bounds() {
551        let problem = OptimProblem::<f64>::new(3)
552            .x0(&[0.0, 0.0, 0.0])
553            .objective(|x: &[f64]| x.iter().copied().map(|xi| xi * xi).sum::<f64>())
554            .bounds(0, (-1.0, 1.0))
555            .bounds(2, (0.0, 10.0));
556
557        assert!(problem.has_bounds());
558        assert!(problem.bounds[0].is_some());
559        assert!(problem.bounds[1].is_none());
560        assert!(problem.bounds[2].is_some());
561    }
562
563    #[test]
564    fn test_builder_with_constraints() {
565        let problem = OptimProblem::<f64>::new(2)
566            .x0(&[1.0, 1.0])
567            .objective(|x: &[f64]| x[0] + x[1])
568            .constraint_eq(|x: &[f64]| x[0] * x[0] + x[1] * x[1] - 1.0)
569            .constraint_ineq(|x: &[f64]| x[0] - 0.5);
570
571        assert!(problem.has_constraints());
572        assert_eq!(problem.n_eq_constraints(), 1);
573        assert_eq!(problem.n_ineq_constraints(), 1);
574    }
575
576    #[test]
577    fn test_builder_least_squares() {
578        let problem = OptimProblem::<f64>::new(2).x0(&[0.0, 0.0]).least_squares(
579            3,
580            |x: &[f64], r: &mut [f64]| {
581                r[0] = x[0] - 1.0;
582                r[1] = x[1] - 2.0;
583                r[2] = x[0] + x[1] - 3.0;
584            },
585        );
586
587        assert!(problem.is_least_squares());
588    }
589
590    #[test]
591    fn test_builder_linear_objective() {
592        let p = OptimProblem::<f64>::new(3).linear_objective(&[2.0, 3.0, 1.0]);
593        assert!(p.is_linear());
594        assert!(!p.is_quadratic());
595        assert!(!p.is_least_squares());
596    }
597
598    #[test]
599    fn test_builder_quadratic_objective() {
600        let p = OptimProblem::<f64>::new(2)
601            .quadratic_objective_dense(&[2.0, 0.0, 0.0, 4.0], &[1.0, 1.0]);
602        assert!(p.is_quadratic());
603        assert!(!p.is_linear());
604    }
605
606    #[test]
607    fn test_builder_linear_constraints() {
608        let p = OptimProblem::<f64>::new(2)
609            .linear_objective(&[1.0, 2.0])
610            .linear_constraint_ineq(&[1.0, 1.0], 10.0)
611            .linear_constraint_eq(&[1.0, -1.0], 0.0);
612        assert!(p.has_linear_constraints());
613        assert_eq!(p.linear_constraints.len(), 2);
614    }
615
616    #[test]
617    fn test_builder_integer_vars() {
618        let p = OptimProblem::<f64>::new(3)
619            .linear_objective(&[1.0, 2.0, 3.0])
620            .integer_var(0)
621            .binary_var(2);
622        assert!(p.has_integer_vars());
623        assert_eq!(p.var_types[0], VarType::Integer);
624        assert_eq!(p.var_types[1], VarType::Continuous);
625        assert_eq!(p.var_types[2], VarType::Binary);
626        assert_eq!(p.bounds[2], Some((0.0, 1.0)));
627    }
628
629    #[test]
630    fn test_finite_diff_gradient() {
631        let f = |x: &[f64]| x[0] * x[0] + 4.0 * x[1] * x[1];
632        let x = [3.0, 2.0];
633        let mut g = [0.0; 2];
634        finite_diff_gradient(&f, &x, &mut g);
635        // df/dx0 = 2*x0 = 6.0, df/dx1 = 8*x1 = 16.0
636        assert!((g[0] - 6.0).abs() < 1e-5);
637        assert!((g[1] - 16.0).abs() < 1e-5);
638    }
639}
numra_optim/problem.rs

numra_optim/
problem.rs
