Skip to main content

gam_problem/
linear_constraints.rs

1use ndarray::{Array1, Array2};
2
3#[derive(Clone, Debug)]
4pub struct LinearInequalityConstraints {
5    pub a: Array2<f64>,
6    pub b: Array1<f64>,
7}
8
9impl LinearInequalityConstraints {
10    /// Construct with the equal-row-count invariant enforced. The dimensions
11    /// `a.nrows() == b.len()` are required by every downstream KKT / active-set
12    /// routine; routing every construction site through this constructor
13    /// eliminates a class of "rows out of sync" bugs at the type boundary.
14    #[inline]
15    pub fn new(a: Array2<f64>, b: Array1<f64>) -> Result<Self, String> {
16        if a.nrows() != b.len() {
17            return Err(format!(
18                "LinearInequalityConstraints: row count mismatch (A has {} rows, b has length {})",
19                a.nrows(),
20                b.len(),
21            ));
22        }
23        Ok(Self { a, b })
24    }
25
26    /// Build the per-coordinate `β_i >= lower_bounds[i]` inequality system.
27    /// Non-finite entries are treated as "no bound" and skipped; returns
28    /// `None` when every entry is non-finite so callers can short-circuit
29    /// the no-constraint case without allocating the empty A/b pair.
30    pub fn from_per_coordinate_lower_bounds(lower_bounds: &Array1<f64>) -> Option<Self> {
31        let active_rows: Vec<usize> = (0..lower_bounds.len())
32            .filter(|&i| lower_bounds[i].is_finite())
33            .collect();
34        if active_rows.is_empty() {
35            return None;
36        }
37        let p = lower_bounds.len();
38        let mut a = Array2::<f64>::zeros((active_rows.len(), p));
39        let mut b = Array1::<f64>::zeros(active_rows.len());
40        for (r, &idx) in active_rows.iter().enumerate() {
41            a[[r, idx]] = 1.0;
42            b[r] = lower_bounds[idx];
43        }
44        Some(Self { a, b })
45    }
46}