problemreductions/models/algebraic/
feasible_basis_extension.rs

1//! Feasible Basis Extension problem implementation.
2//!
3//! Given an m x n integer matrix A (m < n), a column vector a_bar of length m,
4//! and a subset S of column indices with |S| < m, determine whether there exists
5//! a feasible basis B (a set of m column indices including S) such that the
6//! m x m submatrix A_B is nonsingular and A_B^{-1} a_bar >= 0.
7//!
8//! NP-complete (Murty, 1972).
9
10use crate::registry::{FieldInfo, ProblemSchemaEntry};
11use crate::traits::Problem;
12use serde::{Deserialize, Serialize};
13
14inventory::submit! {
15    ProblemSchemaEntry {
16        name: "FeasibleBasisExtension",
17        display_name: "Feasible Basis Extension",
18        aliases: &[],
19        dimensions: &[],
20        module_path: module_path!(),
21        description: "Given matrix A, vector a_bar, and required columns S, find a feasible basis extending S",
22        fields: &[
23            FieldInfo { name: "matrix", type_name: "Vec<Vec<i64>>", description: "m x n integer matrix A (row-major)" },
24            FieldInfo { name: "rhs", type_name: "Vec<i64>", description: "Column vector a_bar of length m" },
25            FieldInfo { name: "required_columns", type_name: "Vec<usize>", description: "Subset S of column indices that must be in the basis" },
26        ],
27    }
28}
29
30/// The Feasible Basis Extension problem.
31///
32/// Given an m x n integer matrix A with m < n, a column vector a_bar of length m,
33/// and a subset S of column indices with |S| < m, determine whether there exists
34/// a feasible basis B of m columns (including all of S) such that the submatrix
35/// A_B is nonsingular and A_B^{-1} a_bar >= 0.
36///
37/// # Representation
38///
39/// Each non-required column has a binary variable: `x_j = 1` if column j is
40/// selected. A valid config must select exactly m - |S| additional columns.
41/// The problem is satisfiable iff the resulting A_B is nonsingular and
42/// A_B^{-1} a_bar >= 0.
43///
44/// # Example
45///
46/// ```
47/// use problemreductions::models::algebraic::FeasibleBasisExtension;
48/// use problemreductions::{Problem, Solver, BruteForce};
49///
50/// let matrix = vec![
51///     vec![1, 0, 1, 2, -1, 0],
52///     vec![0, 1, 0, 1,  1, 2],
53///     vec![0, 0, 1, 1,  0, 1],
54/// ];
55/// let rhs = vec![7, 5, 3];
56/// let required = vec![0, 1];
57/// let problem = FeasibleBasisExtension::new(matrix, rhs, required);
58/// let solver = BruteForce::new();
59/// let solution = solver.find_witness(&problem);
60/// assert!(solution.is_some());
61/// ```
62#[derive(Debug, Clone, Serialize, Deserialize)]
63pub struct FeasibleBasisExtension {
64    matrix: Vec<Vec<i64>>,
65    rhs: Vec<i64>,
66    required_columns: Vec<usize>,
67}
68
69impl FeasibleBasisExtension {
70    /// Create a new FeasibleBasisExtension instance.
71    ///
72    /// # Panics
73    ///
74    /// Panics if:
75    /// - The matrix is empty or has inconsistent row lengths
76    /// - m >= n (must have more columns than rows)
77    /// - rhs length does not equal m
78    /// - |S| >= m (must have room for at least one additional column)
79    /// - Any required column index is out of bounds
80    /// - Required columns contain duplicates
81    pub fn new(matrix: Vec<Vec<i64>>, rhs: Vec<i64>, required_columns: Vec<usize>) -> Self {
82        let m = matrix.len();
83        assert!(m > 0, "Matrix must have at least one row");
84        let n = matrix[0].len();
85        for row in &matrix {
86            assert_eq!(row.len(), n, "All rows must have the same length");
87        }
88        assert!(
89            m < n,
90            "Number of rows ({m}) must be less than number of columns ({n})"
91        );
92        assert_eq!(
93            rhs.len(),
94            m,
95            "rhs length ({}) must equal number of rows ({m})",
96            rhs.len()
97        );
98        assert!(
99            required_columns.len() < m,
100            "|S| ({}) must be less than m ({m})",
101            required_columns.len()
102        );
103        for &col in &required_columns {
104            assert!(col < n, "Required column index {col} out of bounds (n={n})");
105        }
106        // Check for duplicates
107        let mut sorted = required_columns.clone();
108        sorted.sort_unstable();
109        for i in 1..sorted.len() {
110            assert_ne!(
111                sorted[i - 1],
112                sorted[i],
113                "Duplicate required column index {}",
114                sorted[i]
115            );
116        }
117        Self {
118            matrix,
119            rhs,
120            required_columns,
121        }
122    }
123
124    /// Returns the matrix A.
125    pub fn matrix(&self) -> &[Vec<i64>] {
126        &self.matrix
127    }
128
129    /// Returns the right-hand side vector a_bar.
130    pub fn rhs(&self) -> &[i64] {
131        &self.rhs
132    }
133
134    /// Returns the required column indices S.
135    pub fn required_columns(&self) -> &[usize] {
136        &self.required_columns
137    }
138
139    /// Returns the number of rows (m).
140    pub fn num_rows(&self) -> usize {
141        self.matrix.len()
142    }
143
144    /// Returns the number of columns (n).
145    pub fn num_columns(&self) -> usize {
146        self.matrix[0].len()
147    }
148
149    /// Returns the number of required columns (|S|).
150    pub fn num_required(&self) -> usize {
151        self.required_columns.len()
152    }
153
154    /// Returns the indices of non-required columns (the "free" columns).
155    fn free_columns(&self) -> Vec<usize> {
156        let required_set: std::collections::HashSet<usize> =
157            self.required_columns.iter().copied().collect();
158        (0..self.num_columns())
159            .filter(|c| !required_set.contains(c))
160            .collect()
161    }
162
163    /// Check if basis columns form a nonsingular system and the solution is non-negative.
164    ///
165    /// Uses exact rational arithmetic via integer Gaussian elimination with
166    /// numerator/denominator tracking to avoid floating-point errors.
167    #[allow(clippy::needless_range_loop)]
168    fn check_feasible_basis(&self, basis_cols: &[usize]) -> bool {
169        let m = self.num_rows();
170        assert_eq!(basis_cols.len(), m);
171
172        // Build augmented matrix [A_B | a_bar] in i128 to avoid overflow
173        // during Bareiss fraction-free Gaussian elimination.
174        let mut aug128: Vec<Vec<i128>> = Vec::with_capacity(m);
175        for i in 0..m {
176            let mut row = Vec::with_capacity(m + 1);
177            for &col in basis_cols {
178                row.push(self.matrix[i][col] as i128);
179            }
180            row.push(self.rhs[i] as i128);
181            aug128.push(row);
182        }
183
184        // Bareiss algorithm: fraction-free Gaussian elimination.
185        // After elimination, the system is upper-triangular.
186        let mut prev_pivot = 1i128;
187
188        for k in 0..m {
189            // Partial pivoting
190            let mut max_row = k;
191            let mut max_val = aug128[k][k].abs();
192            for i in (k + 1)..m {
193                if aug128[i][k].abs() > max_val {
194                    max_val = aug128[i][k].abs();
195                    max_row = i;
196                }
197            }
198            if max_val == 0 {
199                return false; // singular
200            }
201            if max_row != k {
202                aug128.swap(k, max_row);
203            }
204
205            for i in (k + 1)..m {
206                for j in (k + 1)..=m {
207                    aug128[i][j] =
208                        (aug128[k][k] * aug128[i][j] - aug128[i][k] * aug128[k][j]) / prev_pivot;
209                }
210                aug128[i][k] = 0;
211            }
212            prev_pivot = aug128[k][k];
213        }
214
215        // Back-substitution to solve. We solve in rational form: x_i = num_i / det.
216        // The solution x = A_B^{-1} a_bar must satisfy x >= 0, which means
217        // num_i / det >= 0 for all i, i.e., num_i and det have the same sign (or num_i = 0).
218        // Back-substitution using rational arithmetic to check x >= 0.
219        // Simple rational back-substitution:
220        // x[i] = (aug128[i][m] - sum_{j>i} aug128[i][j] * x[j]) / aug128[i][i]
221        // We track x[i] as (numerator, denominator) pairs.
222
223        let mut x_nums = vec![0i128; m];
224        let mut x_dens = vec![1i128; m];
225
226        for i in (0..m).rev() {
227            // numerator of (aug128[i][m] - sum_{j>i} aug128[i][j] * x[j])
228            let mut num = aug128[i][m];
229            let mut den = 1i128;
230
231            for j in (i + 1)..m {
232                // subtract aug128[i][j] * (x_nums[j] / x_dens[j])
233                // num/den - aug128[i][j] * x_nums[j] / x_dens[j]
234                // = (num * x_dens[j] - den * aug128[i][j] * x_nums[j]) / (den * x_dens[j])
235                let a = aug128[i][j];
236                num = num * x_dens[j] - den * a * x_nums[j];
237                den *= x_dens[j];
238                // Simplify to avoid overflow
239                let g = gcd_i128(num.abs(), den.abs());
240                if g > 1 {
241                    num /= g;
242                    den /= g;
243                }
244            }
245            // x[i] = (num/den) / aug128[i][i] = num / (den * aug128[i][i])
246            let diag = aug128[i][i];
247            x_nums[i] = num;
248            x_dens[i] = den * diag;
249            // Normalize sign: make denominator positive
250            if x_dens[i] < 0 {
251                x_nums[i] = -x_nums[i];
252                x_dens[i] = -x_dens[i];
253            }
254            let g = gcd_i128(x_nums[i].abs(), x_dens[i].abs());
255            if g > 1 {
256                x_nums[i] /= g;
257                x_dens[i] /= g;
258            }
259        }
260
261        // Check x >= 0: each x_nums[i] / x_dens[i] >= 0
262        // Since x_dens[i] > 0 (normalized), we need x_nums[i] >= 0
263        x_nums.iter().take(m).all(|&num| num >= 0)
264    }
265}
266
267/// Compute GCD of two i128 values.
268fn gcd_i128(mut a: i128, mut b: i128) -> i128 {
269    while b != 0 {
270        let t = b;
271        b = a % b;
272        a = t;
273    }
274    a
275}
276
277impl Problem for FeasibleBasisExtension {
278    const NAME: &'static str = "FeasibleBasisExtension";
279    type Value = crate::types::Or;
280
281    fn variant() -> Vec<(&'static str, &'static str)> {
282        crate::variant_params![]
283    }
284
285    fn dims(&self) -> Vec<usize> {
286        vec![2; self.num_columns() - self.num_required()]
287    }
288
289    fn evaluate(&self, config: &[usize]) -> crate::types::Or {
290        let free_cols = self.free_columns();
291        let num_free = free_cols.len();
292
293        if config.len() != num_free {
294            return crate::types::Or(false);
295        }
296        if config.iter().any(|&v| v >= 2) {
297            return crate::types::Or(false);
298        }
299
300        let m = self.num_rows();
301        let s = self.num_required();
302        let needed = m - s;
303
304        // Count selected free columns
305        let selected_free: Vec<usize> = config
306            .iter()
307            .enumerate()
308            .filter(|(_, &v)| v == 1)
309            .map(|(i, _)| free_cols[i])
310            .collect();
311
312        if selected_free.len() != needed {
313            return crate::types::Or(false);
314        }
315
316        // Form basis: required columns + selected free columns
317        let mut basis_cols: Vec<usize> = self.required_columns.clone();
318        basis_cols.extend_from_slice(&selected_free);
319
320        crate::types::Or(self.check_feasible_basis(&basis_cols))
321    }
322}
323
324crate::declare_variants! {
325    default FeasibleBasisExtension => "2^num_columns * num_rows^3",
326}
327
328#[cfg(feature = "example-db")]
329pub(crate) fn canonical_model_example_specs() -> Vec<crate::example_db::specs::ModelExampleSpec> {
330    vec![crate::example_db::specs::ModelExampleSpec {
331        id: "feasible_basis_extension",
332        // 3x6 matrix, rhs=[7,5,3], required={0,1}, select col 2 -> B={0,1,2}, x=(4,5,3)>=0
333        instance: Box::new(FeasibleBasisExtension::new(
334            vec![
335                vec![1, 0, 1, 2, -1, 0],
336                vec![0, 1, 0, 1, 1, 2],
337                vec![0, 0, 1, 1, 0, 1],
338            ],
339            vec![7, 5, 3],
340            vec![0, 1],
341        )),
342        optimal_config: vec![1, 0, 0, 0], // select col 2 (first free column)
343        optimal_value: serde_json::json!(true),
344    }]
345}
346
347#[cfg(test)]
348#[path = "../../unit_tests/models/algebraic/feasible_basis_extension.rs"]
349mod tests;
problemreductions/models/algebraic/feasible_basis_extension.rs

problemreductions/models/algebraic/
feasible_basis_extension.rs
